生物群集の多変量解析　小林四朗先生　蒼樹書房

```{r}

#主成分分析 p.120

#matの縦は地点、横は種

mat=cbind(c(2,1,3,6,4),c(1,3,4,3,2))

L=nrow(mat)

n_ave=apply(mat,2,mean)

A=(t(mat)-n_ave)/sqrt(L-1)

#分散共分散
R=A%*%t(A)

eig_values=eigen(R)$values

eig_vectors=eigen(R)$vector

factor=eig_vectors

for(j in 1:ncol(factor)){
  
factor[,j]=eig_vectors[,j]*sqrt(eig_values[j])/sqrt(diag(R))

}


```

```{r}

#主成分分析

#Frame 法

#matの縦は地点、横は種

mat=cbind(c(2,1,3,6,4),c(1,3,4,3,2))

L=nrow(mat)

n_ave=apply(mat,2,mean)

A=(t(mat)-n_ave)/sqrt(L-1)

#分散共分散
A=A%*%t(A)

#A=matrix(c(202,-176,228,-262,440,-462,-369,504,-561),ncol=3)

row=nrow(A)

a=c(-sum(diag(A)))

A2=A

for(j in 1:(row-1)){

As=A2+diag(a[length(a)],row)

A2=A%*%As

a=c(a,-sum(diag(A2))/(j+1))

}

lam=Re(sort(polyroot(c(a[length(a):1],1)),decreasing=T)) 

Q=function(x){

mat=diag(x^(row-1),row)

a=c(-sum(diag(A)))

A2=A

for(j in 1:(row-1)){

As=A2+diag(a[length(a)],row)

mat=mat+As*x^(row-1-j)

A2=A%*%As

a=c(a,-sum(diag(A2))/(j+1))

}

return(mat)

}


v=array(0,dim=dim(A))

for(j in 1:row){

v[,j]=Q(lam[j])[,1]

}



```

```{r}
#主成分分析
#QR法　
#固有値(A_matの対角成分)と固有ベクトル(Q)
#matの縦は地点、横は種
mat=cbind(c(2,1,3,6,4),c(1,3,4,3,2))
L=nrow(mat)
n_ave=apply(mat,2,mean)
A=(t(mat)-n_ave)/sqrt(L-1)
#分散共分散
A=A%*%t(A)
A_mat=A
ite=1000
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
q_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
r_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
q1=A_mat[,1]/sqrt(sum(A_mat[,1]^2))
q_mat[,1]=q1
r_mat[1,1]=sqrt(sum(A_mat[,1]^2))
for(j in 2:ncol(A_mat)){
u=A_mat[,j]-apply(t(t(q_mat[,1:(j-1)])*c(t(q_mat[,1:(j-1)])%*%A_mat[,j])),1,sum)
r_mat[j,j]=sqrt(sum(u^2))
q_mat[,j]=u/sqrt(sum(u^2))
}
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}
```


```{r}



#主座標分析 p.131

mat=cbind(c(2,1,3,6,4),c(1,3,4,3,2))

L=nrow(mat)

d=array(0,dim=c(nrow(mat),nrow(mat)))

for(i in 1:ncol(d)){
for(j in 1:nrow(d)){  
  
d[i,j]=sum(abs(mat[i,]-mat[j,]))  
  
}}

del=d

A=(t(t(-del^2+apply(del^2,2,sum)/L)+apply(del^2,1,sum)/L)-sum(del^2)/(L^2))/2

eig=eigen(A)$values

prop=eig/(sum(eig[1:(length(eig)-1)])+(L-1)*abs(eig[length(eig)]))

prop=cumsum(prop)

eig_vector=eigen(A)$vector[,1:ifelse(length(prop[prop<0.8])>2,length(prop[prop<0.8]),2)]

eig_vector=eig_vector*c(sqrt(eig[1:ifelse(length(prop[prop<0.8])>2,length(prop[prop<0.8]),2)]))



```

```{r}

#主座標分析

#Frame 法

mat=cbind(c(2,1,3,6,4),c(1,3,4,3,2))

L=nrow(mat)

d=array(0,dim=c(nrow(mat),nrow(mat)))

for(i in 1:ncol(d)){
for(j in 1:nrow(d)){  

d[i,j]=sum(abs(mat[i,]-mat[j,]))  

}}

del=d

A=(t(t(-del^2+apply(del^2,2,sum)/L)+apply(del^2,1,sum)/L)-sum(del^2)/(L^2))/2

row=nrow(A)

a=c(-sum(diag(A)))

A2=A

for(j in 1:(row-1)){

As=A2+diag(a[length(a)],row)

A2=A%*%As

a=c(a,-sum(diag(A2))/(j+1))

}

lam=Re(sort(polyroot(c(a[length(a):1],1)),decreasing=T)) 

Q=function(x){

mat=diag(x^(row-1),row)

a=c(-sum(diag(A)))

A2=A

for(j in 1:(row-1)){

As=A2+diag(a[length(a)],row)

mat=mat+As*x^(row-1-j)

A2=A%*%As

a=c(a,-sum(diag(A2))/(j+1))

}

return(mat)

}


v=array(0,dim=dim(A))

for(j in 1:row){

v[,j]=Q(lam[j])[,1]

}



```

```{r}

#交互平均法（対応分析）p.137

nij=matrix(c(1,0,0,2,4,0,0,1,0,1,2,0,0,1,0,2,2,1,0,0,0,2,3,0,1,2,5,3,1,0,0,3,5,1,0,1,0,0,0,1,3,0,0,0,2,0,0,0),ncol = 6)

I=nrow(nij);J=ncol(nij)

w=sample(1:9,I,replace=T)

w_pre=w+10

while(sum(abs(w_pre-w))>10^(-6))({
  
w_pre=w

y=t(nij)%*%w/apply(nij,2,sum)

x=nij%*%y/apply(nij,1,sum)  

w=100*(x-min(x))/(max(x)-min(x))
  
print(sum(abs(w-w_pre)))

})

E=(max(x)-min(x))/100

#第1軸
Y1=100*(y-min(y))/(max(y)-min(y))

R=apply(nij,1,sum)

z=sum(R*w)/sum(R)

V=R*w-R*z

w2=w-sample(1:9,length(w),replace=T)/10

u=sum(V*w2)/sum(V*w)

w2=w2-u*w

w2_pre=w2+10

while(sum(abs(w2_pre-w2))>10^(-4))({
  
w2_pre=w2

y2=t(nij)%*%w2/apply(nij,2,sum)

x2=nij%*%y2/apply(nij,1,sum)  

w2=100*(x2-min(x2))/(max(x2)-min(x2))
  
V=R*w2-R*z

#u=sum(V*w2)/sum(V*w)

#w2=w2-u*w

print(sum(abs(w2-w2_pre)))

})

E2=(max(x2)-min(x2))/100

#第2軸
Y2=100*(y2-min(y2))/(max(y2)-min(y2))

```

```{r}

#p.149 ニッチ分析


#p.157

#H:2種が完全に重なり合い、利用曲線が等しい

mat=matrix(c(29,12,2,9,26,5,12,21,1,15,15,10,3,16,34,9,4,18),ncol=6)

S=dim(mat)[1];L=dim(mat)[2]

Nj=apply(mat,1,sum)

pij=mat/c(Nj)

E=diag(0,nrow(mat))

for(h in 1:nrow(E)){
for(i in 1:ncol(E)){  
  
E[h,i]=sum(pij[h,]*log(pij[i,]))-sum(pij[h,]*log(pij[h,]))  
  
}}

SO=exp(E)

U=-2*(log(SO))*c(Nj)

qchisq(1-0.05,df=L-1)


#H:種の利用曲線が完全に重なり合う

t=sum(Nj);Mi=apply(mat,2,sum)
  
E_val=sum(mat*log(c(Nj/t))-mat*log(t(t(mat)/c(Mi))))/t

GO=exp(E_val)

V=-2*t*log(GO)

qchisq(1-0.05,df=(S-1)*(L-1))


```