実験計画　生沢雅夫先生　新曜社

```{r}

#実験計画　生沢雅夫著

#3章　線形模型における仮説検定

library(dplyr)

v=c(1,2,3,3.7)

x=c(50,60,84,97)

dim=2

pthi=array(0,dim=c(length(v),dim+1))

pthi[,1]=rep(1,length(v))

for(j in 1:dim){

pthi[,j+1]=v^j

}

dim1=dim

params1=solve(t(pthi)%*%pthi)%*%t(pthi)%*%x

#誤差
sum(x^2)-as.numeric(t(params1)%*%t(pthi)%*%x)

SS1=as.numeric(t(params1)%*%t(pthi)%*%x)



v=c(1,2,3,3.7)

x=c(50,60,84,97)

dim=1

pthi=array(0,dim=c(length(v),dim+1))

pthi[,1]=rep(1,length(v))

for(j in 1:dim){

pthi[,j+1]=v^j

}

dim2=dim

params2=solve(t(pthi)%*%pthi)%*%t(pthi)%*%x

#誤差
sum(x^2)-as.numeric(t(params2)%*%t(pthi)%*%x)

SS2=as.numeric(t(params2)%*%t(pthi)%*%x)

SSR=SS2-SS1;MSR=SSR/(length(x)-(dim2+dim1))

MS2=(SS2/dim2)

f=MS2/MSR

qf(1-0.05,dim2,length(x)-(dim2+dim1))

#4章　完全無作為配置

dat=data.frame(group=c(rep(1,3),rep(2,4),rep(3,4)),value=c(5,8,8,9,11,8,8,12,10,5,13))

groups=sort(unique(dat$group))

CF=sum(dat$value)^2/nrow(dat)

SST=sum(dat$value^2)-CF

summary=dat%>%group_by(group)%>%summarise(value=sum(value),ave=mean(value),r=n())

SSA=sum(summary$value^2/summary$r)-CF

SSR=SST-SSA

MSA=SSA/(length(groups)-1)

MSR=SSR/(nrow(dat)-length(groups))

#H：各群の水準は等しい

f=MSA/MSR

qf(1-0.05,length(groups)-1,nrow(dat)-length(groups))

#Hが棄却された場合：

t=array(0,dim=c(length(groups),length(groups)))

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

sub1=summary%>%filter(group==i)

sub2=summary%>%filter(group==j)

t[i,j]=(sub1$ave-sub2$ave)/(MSR*sqrt(1/sub1$r+1/sub2$r))

}}

qf(1-0.05,1,nrow(dat)-length(groups))

#p.64 4-6 固定効果模型と変量効果模型

dat=data.frame(group=c(rep(1,4),rep(2,4),rep(3,4)),value=c(11,10,7,12,0,2,0,-2,7,1,3,-1))

CF=sum(dat$value)^2/nrow(dat)

y=dat$value

groups=sort(unique(dat$group))

SST=sum(y^2)-CF

summary=dat%>%group_by(group)%>%summarise(value=sum(value),r=n())

SSA=sum(summary$value^2/summary$r)-CF

SSR=SST-SSA

MSA=SSA/(length(groups)-1)

MSR=SSR/(nrow(dat)-length(groups))

f=MSA/MSR

qf(1-0.05,length(groups)-1,nrow(dat)-length(groups))

sig_hat=MSR

sig_a_hat=(MSA-MSR)/4

#5.乱塊配置　p.73

data=data.frame(factor=c(1:3),matrix(c(93.75,154.38,206.25,100.13,161.25,208.75,107.5,161.88,220,101.25,157.5,232.5),ncol=4))

mat=as.matrix(data[,-1])-100

As=apply(mat,1,mean)-mean(mat)

Bs=apply(mat,2,mean)-mean(mat)

mu=mean(mat)

CF=sum(mat)^2/prod(dim(mat))

SSA=sum(apply(mat,1,sum)^2/ncol(mat))-CF

SSB=sum(apply(mat,2,sum)^2/nrow(mat))-CF

SSR=sum(mat^2)-CF-SSA-SSB

MSA=SSA/(nrow(mat)-1)

MSB=SSB/(ncol(mat)-1)

MSR=SSR/((dim(mat)[1]-1)*(dim(mat)[2]-1))

SST=sum(mat^2)-CF

#要因Aの検定

f_A=MSA/MSR

qf(1-0.01,nrow(mat)-1,(dim(mat)[1]-1)*(dim(mat)[2]-1))

#要因Bの検定

f_B=MSB/MSR

qf(1-0.01,ncol(mat)-1,(dim(mat)[1]-1)*(dim(mat)[2]-1))


#6.2要因配置

r=4;p=2;q=3

mat1=matrix(c(3,5,5,6,8,7,9,8),ncol=r)

rownames(mat1)=c("A1","A2")

colnames(mat1)=rep("B1",r)

mat2=matrix(c(7,5,9,8,11,9,12,10),ncol=r)

rownames(mat2)=c("A1","A2")

colnames(mat2)=rep("B2",r)

mat3=matrix(c(7,7,8,9,11,10,11,11),ncol=r)

rownames(mat3)=c("A1","A2")

colnames(mat3)=rep("B3",r)

mat=cbind(mat1,mat2,mat3)

#合計

total=cbind(apply(mat1,1,sum),apply(mat2,1,sum),apply(mat3,1,sum))

#平均

average=cbind(apply(mat1,1,mean),apply(mat2,1,mean),apply(mat3,1,mean))

mu_hat=mean(mat)

alpha_A=apply(total,1,sum)/(q*r)-mu_hat

alpha_B=apply(total,2,sum)/(p*r)-mu_hat

alpha_AB=total/r

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

alpha_AB[i,j]=alpha_AB[i,j]-alpha_A[i]-alpha_B[j]-mu_hat

}}

CF=sum(mat)^2/prod(dim(mat))

SST=sum(mat^2)-CF

SSA=sum(apply(total,1,sum)^2)/(q*r)-CF

SSB=sum(apply(total,2,sum)^2)/(p*r)-CF

SSAB=sum(total^2)/r-sum(apply(total,1,sum)^2)/(q*r)-sum(apply(total,2,sum)^2)/(p*r)+CF

SSR=sum(mat^2)-sum(total^2)/r

MSA=SSA/(p-1);MSB=SSB/(q-1);MSAB=SSAB/((p-1)*(q-1));MSR=SSR/(p*q*(r-1))

n=prod(dim(mat))

fA=MSA/MSR;

qf(1-0.01,p-1,n-p*q)

fB=MSB/MSR;

qf(1-0.01,q-1,n-p*q)

fAB=MSAB/MSR;

qf(1-0.01,(p-1)*(q-1),n-p*q)

#2要因配置で要因の一方が変量模型の時

r=4;p=2;q=3

mat1=matrix(c(19,46,20,49,27,43,22,46),ncol=r)

rownames(mat1)=c("A1","A2")

colnames(mat1)=rep("B1",r)

mat2=matrix(c(25,41,20,48,19,39,20,40),ncol=r)

rownames(mat2)=c("A1","A2")

colnames(mat2)=rep("B2",r)

mat3=matrix(c(16,33,20,41,26,42,20,34),ncol=r)

rownames(mat3)=c("A1","A2")

colnames(mat3)=rep("B3",r)

mat=cbind(mat1,mat2,mat3)

#合計

total=cbind(apply(mat1,1,sum),apply(mat2,1,sum),apply(mat3,1,sum))

#平均

average=cbind(apply(mat1,1,mean),apply(mat2,1,mean),apply(mat3,1,mean))

CF=sum(mat)^2/prod(dim(mat))

SST=sum(mat^2)-CF

SSA=sum(apply(total,1,sum)^2)/(q*r)-CF

SSB=sum(apply(total,2,sum)^2)/(p*r)-CF

SSAB=sum(total^2)/r-sum(apply(total,1,sum)^2)/(q*r)-sum(apply(total,2,sum)^2)/(p*r)+CF

SSR=sum(mat^2)-sum(total^2)/r

MSA=SSA/(p-1);MSB=SSB/(q-1);MSAB=SSAB/((p-1)*(q-1));MSR=SSR/(p*q*(r-1))

n=prod(dim(mat))

fA=MSA/MSAB;

qf(1-0.01,p-1,(p-1)*(q-1))

fB=MSB/MSR;

qf(1-0.01,q-1,n-p*q)

fAB=MSAB/MSR;

qf(1-0.01,(p-1)*(q-1),n-p*q)







```

```{r}

#p.100 2要因配置　繰り返し数の異なる場合

library(dummies)

A=c(rep(1,8),rep(2,9))
B=c(1,1,1,2,2,3,3,3,1,1,2,2,2,2,3,3,3)
AB=c(1,1,1,2,2,3,3,3,4,4,5,5,5,5,6,6,6)
mat=cbind(dummy(A),dummy(B),dummy(AB))
X=c(3,8,9,9,11,7,8,11,5,8,5,8,9,10,7,10,11)

r=apply(mat,2,sum)
r=r[-c(1:(length(unique(A))+length(unique(B))))]
r=t(matrix(r,ncol=length(unique(A))))

pthix=t(mat)%*%X

S=pthix[1:length(unique(A))]-r%*%c(t(dummy(B))%*%X/apply(dummy(B),2,sum))

r1=apply(r,1,sum)
r2=apply(r,2,sum)

Mat=array(0,dim=c(length(unique(A)),length(unique(A))))
for(i in 1:length(unique(A))){
for(j in 1:length(unique(A))){  
Mat[i,j]=-sum(r[i,]*r[j,]/r2)
if(i==j){Mat[i,j]=Mat[i,j]+r1[i]}
}}

#修正ガウス・ザイデル法
AA=Mat　#行列のサンプル
C=S  #ベクトルのサンプル
rs=C/diag(AA)　#r_iを計算
d=AA/c(diag(AA))　#d_{i,j}を計算
diag(d)=0 #d_{i,j}を計算
mat=cbind(rs,-d) #x_iの更新に必要な行列を作成
ite=1000　#反復計算回数
x=rep(1,ncol(AA))　#x_iの初期値
eta=0.01
for(l in 1:ite){
xpre=x  
z=c(1,x)
x=x*(1-eta)+eta*mat%*%z
x=x-mean(x)
print(sum(abs(xpre-x)))
}　

alpha=x
mu_B=c(t(dummy(B))%*%X/apply(dummy(B),2,sum))-t(r)%*%alpha/r2

Mat2=array(0,dim=c(length(unique(A)),length(unique(B))))
for(i in 1:length(unique(A))){
for(j in 1:length(unique(B))){  
  Mat2[i,j]=alpha[i]+mu_B[j]
}}

alpha_beta=t(matrix(t(dummy(AB))%*%X/apply(dummy(AB),2,sum),ncol=length(unique(A))))-Mat2

PB=dummy(B)%*%solve(t(dummy(B))%*%dummy(B))%*%t(dummy(B))
Ppthi=cbind(dummy(A),dummy(B))%*%solve(t(cbind(dummy(A),dummy(B)))%*%cbind(dummy(A),dummy(B)))%*%t(cbind(dummy(A),dummy(B)))


```