分散分析（多変量解析入門Ⅱ　河口至商先生）

分散分析のコードです。

```{r}

#分散分析(一元配置)
library(dplyr)
#工場における不良率の差異を調べる
data=data.frame(num=1:6,factory1=c(3.1,4.1,3.3,3.9,3.7,2.4),factory2=c(4.7,5.6,4.3,5.9,6.1,4.2),factory3=c(5.1,3.7,4.5,6,3.9,5.4))
ave_x=apply(data[,2:ncol(data)],2,mean)
X_bar=mean(ave_x)
matrix=array(0,dim=c(3,3))
S_p=nrow(data)*sum((ave_x-X_bar)^2)
S_E=0
for(j in 1:(ncol(data)-1)){
S_E=S_E+sum((data[,2:ncol(data)][,j]-ave_x[j])^2)
}

p=0.05
K=ncol(data)-1
L=nrow(data)-1
N=nrow(data)
V_p=S_p/(K-1)
V_E=S_E/((N-1)*K)
F_0=V_p/V_E
if(qf(1-p, K-1, K*L)>F_0){
print(paste0("棄却率",p*100,"パーセントで不良率に差はない") ) 
}else{
print(paste0("棄却率",p*100,"パーセントで不良率に差があるといえる")  )
}

#一元配置の必要なサンプルサイズを調べる
#サンプルサイズの決め方 朝倉書店 永田靖先生 p.156

data=data.frame(num=1:6,factory1=c(3.1,4.1,3.3,3.9,3.7,2.4),factory2=c(4.7,5.6,4.3,5.9,6.1,4.2),factory3=c(5.1,3.7,4.5,6,3.9,5.4))
ave_x=apply(data[,2:ncol(data)],2,mean)
X_bar=mean(ave_x)
k=ncol(data)-1　#3群
d=ave_x-X_bar
#分散が未知なのでとりあえず最尤推定量にする
sig=sum(apply((t(data[,-1])-ave_x)^2,1,sum))/(k*nrow(data))

#標本数5~20で確認
ns=seq(5,20,1)
alpha=0.01
pthiA=k-1
power=c()
for(l in 1:length(ns)){
n=ns[l]
pthiE=k*(n-1)
w=qf(1-alpha,pthiA,pthiE)
delta=sum(d^2)/(2*sig)
lam=n*delta
CA=(pthiA+2*lam)/(pthiA+lam)
pthiA_star=(pthiA+lam)^2/(pthiA+2*lam)
u=(sqrt(w/pthiE)*sqrt(2*pthiE-1)-sqrt(CA/pthiA)*sqrt(2*pthiA_star-1))/sqrt(CA/pthiA+w/pthiE)
power=c(power,1-pnorm(u))
}

power_dat=data.frame(n=ns,power)
plot(ns,power,xlab="n",ylab="power",type="o",main=paste0("alpha:",alpha))




```

```{r}

#エンドウ豆の肥料と品種における収量の比較

K=2

cultivars1=data.frame(fertilizer1=c(73,78,77),fertilizer2=c(74,62,65),fertilizer3=c(82,101,78),fertilizer4=c(66,82,70),fertilizer5=c(66,89,93))

cultivars2=data.frame(fertilizer1=c(78,65,80),fertilizer2=c(60,81,64),fertilizer3=c(79,96,83),fertilizer4=c(80,69,63),fertilizer5=c(85,83,64))

data=rbind(cultivars1,cultivars2)-70

cultivars1=cultivars1-70;cultivars2=cultivars2-70

summary=rbind(apply(cultivars1,2,sum),apply(cultivars2,2,sum))

summary=rbind(summary,apply(summary,2,sum),apply(summary,2,sum)^2)

L=ncol(data);N=nrow(cultivars1);

total_var_x=sum((data-sum(data)/(L*N*K))^2)

S_p=L*N*((sum(cultivars1)/(L*N)-sum(data)/(L*N*K))^2+(sum(cultivars2)/(L*N)-sum(data)/(L*N*K))^2)

S_Q=K*N*sum(((apply(cultivars1,2,sum)+apply(cultivars2,2,sum))/(K*N)-sum(data)/(L*N*K))^2)

x_12=rbind(apply(cultivars1,2,mean),apply(cultivars2,2,mean))

x_1=c(sum(cultivars1)/(L*N),sum(cultivars2)/(L*N))

x_2=(apply(cultivars1,2,sum)+apply(cultivars2,2,sum))/(K*N)

S_I=0

for(i in 1:nrow(x_12)){
for(j in 1:ncol(x_12)){
  
 S_I=S_I+N*(x_12[i,j]-x_1[i]-x_2[j]+sum(data)/(L*N*K))^2

}
}

S_E=0

for(j in 1:ncol(cultivars1)){
  
 S_E=S_E+sum((cultivars1[,j]-x_12[1,j])^2)
 
 S_E=S_E+sum((cultivars2[,j]-x_12[2,j])^2)

}

summary=array(0,dim=c(4,5))

summary[1,1]=S_p;summary[2,1]=S_Q;summary[3,1]=S_I;summary[4,1]=S_E
summary[1,2]=K-1;summary[2,2]=L-1;summary[3,2]=(K-1)*(L-1);summary[4,2]=K*L*(N-1)

summary[,3]=summary[,1]/summary[,2]

V_p=S_p/(K-1);V_E=S_E/(K*L*(N-1));V_Q=S_Q/(L-1);V_I=S_I/((K-1)*(L-1))

vec_F=c(V_p,V_Q,V_I)/V_E

p=0.05

for(j in 1:(nrow(summary)-1)){
  
summary[j,4]=vec_F[j]
  
summary[j,5]=qf(1-p,summary[j,2],summary[nrow(summary),2])  
  
}

summary=cbind(c("cultivars","fetilizers","interaction","level"),summary)

```