経時データ分析　藤越康祝・菅民郎・土方裕子先生　オーム社

```{r}

#p.52

library(dplyr)

mat=matrix(c(117,134,133,133,126,117,119,123,126,128,117,114,122,122,117,120,115,108),ncol=3)

n=nrow(mat);p=ncol(mat)

sb=n*sum((apply(mat,2,mean)-mean(mat))^2)

sd=p*sum((apply(mat,1,mean)-mean(mat))^2)

se=sum((t(mat-apply(mat,1,mean))-apply(mat,2,mean)+mean(mat))^2)

st=sum((mat-mean(mat))^2)

#H:bは等しい

f1=(sb/(p-1))/(se/((p-1)*(n-1)))

qf(1-0.05,p-1,(p-1)*(n-1))

#H:sig_d=0

f2=(sd/(n-1))/(se/((p-1)*(n-1)))

qf(1-0.05,n-1,(p-1)*(n-1))

sig=se/((n-1)*(p-1))

sig_d=(sd/(n-1)-se/((n-1)*(p-1)))/p

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

var_b=sig*(1-1/p)/n

#共分散構造の検証

sig_s=array(sig_d,dim=c(p,p))+diag(rep(sig,p))

C=t(matrix(c(1,-1,1,1,-2,0),ncol=3))

C=t(t(C)/sqrt(apply(C^2,2,sum)))

S=cov(mat)

lam=(det(t(C)%*%S%*%C)/((sum(diag(t(C)%*%S%*%C))/(p-1))^(p-1)))^(n/2)

f=(p-2)*(p+1)/2

kai2=-2*(n-1)*(1-(2*p^2-3*p+3)/(6*(p-1)*(n-1)))*log(lam)/n

qchisq(1-0.05,df=f)

```

```{r}

#多群の場合 p.56

#交互作用がある場合の分散分析

mat=matrix(c(rep(1,6),rep(2,7),156,128,136,155,156,139,127,147,143,136,147,122,130,149,112,125,138,149,137,124,151,140,137,151,117,139,136,117,117,145,136,135,124,107,131,117,125,103,119),ncol=4)

mat=cbind(c(c(1:6),c(1:7)),mat)

q=2;p=3;n=nrow(mat)

colnames(mat)=c("被験","群",paste0("t",c(1:3)))

dat=data.frame(mat)

group=dat%>%group_by(群)%>%summarise(n=n())

group_ave=c()

for(j in 1:q){
  
sub=dat%>%filter(群==j)%>%select(-被験,-群)

group_ave=c(group_ave,mean(as.matrix(sub)))
  
}

yij_ave=apply(mat[,-c(1,2)],1,mean)

yjk=dat%>%group_by(群)%>%summarise(t1=mean(t1),t2=mean(t2),t3=mean(t3))

sd=p*sum((yij_ave-group_ave[c(mat[,2])])^2)

sa=sum(group$n*(group_ave-mean(mat[,-c(1:2)]))^2)*p

sb=n*sum((apply(mat[,-c(1:2)],2,mean)-mean(mat[,-c(1:2)]))^2)

sr=sum(t((as.matrix(yjk[,-1])-(t(t(matrix(rep(group_ave,(ncol(mat)-2)),ncol=ncol(mat)-2))+apply(mat[,-c(1:2)],2,mean)))+mean(mat[,-c(1:2)]))^2)%*%group$n)

st=sum((mat[,-c(1:2)]-mean(mat[,-c(1:2)]))^2)

se=st-(sa+sb+sr+sd)

#群間の有意性検定

fa=(sa/(q-1))/(sd/(n-q))

qf(1-0.05,q-1,n-q)

#群と時点の交互作用の検定

fr=(sr/((q-1)*(p-1)))/(se/((p-1)*(n-q)))

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

#時点間の有意性検定

fb=(sb/(p-1))/(se/((p-1)*(n-q)))

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

#個体差の有意性検定

fd=(sd/((n-q)))/(se/((p-1)*(n-q)))

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

```


```{r}

#p.71

library(dplyr)

mat_A=matrix(c(43,45,55,34,40,45,31,43,50,51,45,45,59,39,62,28,47,48,80,48,55,70,50,52,35,50,62),ncol=3)

mat_B=matrix(c(20,38,16,30,43,41,40,21,33,25,45,40,23,30,50,40,23,33,50,50,3,50,55,56),ncol=3)

n1=nrow(mat_A);n2=nrow(mat_B);n=n1+n2;m=n-2

data=data.frame(class=c(rep("A",9),rep("B",8)),rbind(mat_A,mat_B))

colnames(data)=c("class","1times","2times","3times")

summary=as.matrix(cbind(c("A","B"),rbind(apply(mat_A,2,mean),apply(mat_B,2,mean))))

colnames(summary)=c("class","1times","2times","3times")

D=array(0,dim=c(4,3))

for(j in 1:ncol(D)){

D[j:(j+1),j]=c(1,-1)

}

D=t(D)

#H0:平行性仮説

S=(cov(mat_A)+cov(mat_B))/m

b=sqrt(n1*n2/(n1+n2))

u=apply(mat_A,2,mean)-apply(mat_B,2,mean)

t=(b^2)*t(t(D)%*%u)%*%solve(t(D)%*%S%*%D)%*%t(D)%*%u


```


```{r}

#プロファイル分析モデル　p.71

#2群の場合

library(dplyr)

#がくやべんの特徴が種と交互作用があるかどうか（あると予想）

#帰無仮説：種と交互作用はない（平行）

data(iris)

Speciess=unique(iris$Species)

iris1=iris%>%filter(Species==Speciess[1])%>%select(-Species)

iris2=iris%>%filter(Species==Speciess[2])%>%select(-Species)

n1=nrow(iris1);n2=nrow(iris2);n=n1+n2

p=ncol(iris1)

D=array(0,dim=c(p,p-1))

for(j in 1:(p-1)){
  
D[j:(j+1),j]=c(1,-1)  
  
}

y1=apply(iris1,2,mean)

y2=apply(iris2,2,mean)

m=n1+n2-2

S=(cov(iris1)*(n1-1)+cov(iris2)*(n2-1))/m

b=sqrt((n1*n2)/(n1+n2))

u=y1-y2

t=(b^2)*t(t(D)%*%u)%*%solve(t(D)%*%S%*%D)%*%t(D)%*%u

test_statics=(m-p+2)*t/(m*(p-1))

#Trueなら有意差あり

test_statics>qf(1-0.05,p-1,m-p+2)


```


```{r}

#4.3 多群の場合 p.76

#平行性仮説と検定 

library(dplyr)

#がくやべんの特徴が種と交互作用があるかどうか（あると予想）

#帰無仮説H1：種と交互作用はない（平行）

data(iris)

ns=iris%>%group_by(Species)%>%summarise(n=n())

Speciess=ns$Species;k=length(Speciess)

p=ncol(iris)-1

D=array(0,dim=c(p,p-1))

for(j in 1:(p-1)){
  
D[j:(j+1),j]=c(1,-1)  
  
}

Sh=array(0,dim=c(p,p))

Se=array(0,dim=c(p,p))

n=nrow(iris);ns=ns$n;m=n-2

y_ave=apply(iris[,!(colnames(iris) %in% c("Species"))],2,mean)

for(j in 1:k){
  
iris_sub=iris%>%filter(Species==Speciess[j])%>%select(-Species)
 
Sh=Sh+ns[j]*t(t(apply(iris_sub,2,mean)-y_ave))%*%t(apply(iris_sub,2,mean)-y_ave)
  
Se=Se+cov(iris_sub)*(nrow(iris_sub)-1)
  
}

lam1=det(t(D)%*%Se%*%D)/det(t(D)%*%(Se+Sh)%*%D)

test=-(m+(k-1-(p-1)-1)/2)*log(lam1)

#Trueなら有意差あり

test>qchisq(1-0.05,df=(p-1)*(k-1))

```

```{r}

#プロファイル分析モデル　p.71

#2群の場合

library(dplyr)

#帰無仮説：教育方法とテストの平均点に対する交互作用はない

mat=matrix(c(43,45,55,34,40,45,31,43,50,20,38,16,30,43,41,40,21,51,45,45,59,39,62,28,47,48,33,25,45,40,23,30,50,40,80,48,55,70,50,52,35,50,62,23,33,50,50,35,50,55,56,78,60,56,73,55,60,41,60,65,35,40,62,60,50,55,45,35),ncol=4)

iris=data.frame(mat,Species=c(rep(1,9),rep(2,8)))

ns=iris%>%group_by(Species)%>%summarise(n=n())

Speciess=unique(iris$Species)

iris1=iris%>%filter(Species==Speciess[1])%>%select(-Species)

iris2=iris%>%filter(Species==Speciess[2])%>%select(-Species)

n1=nrow(iris1);n2=nrow(iris2);n=n1+n2

p=ncol(iris1)

D=array(0,dim=c(p,p-1))

for(j in 1:(p-1)){
  
D[j:(j+1),j]=c(1,-1)  
  
}

y1=apply(iris1,2,mean)

y2=apply(iris2,2,mean)

m=n1+n2-2

S=(cov(iris1)*(n1-1)+cov(iris2)*(n2-1))/m

b=sqrt((n1*n2)/(n1+n2))

u=y1-y2

t=(b^2)*t(t(D)%*%u)%*%solve(t(D)%*%S%*%D)%*%t(D)%*%u

test_statics=(m-p+2)*t/(m*(p-1))

#Trueなら有意差あり

test_statics>qf(1-0.05,p-1,m-p+2)


r=t(rep(1,p))%*%solve(S)%*%u/as.numeric(t(rep(1,p))%*%solve(S)%*%t(t(rep(1,p))))

```

```{r}

#教育方法とテストの点のデータについて、多群として扱ってみる


mat=matrix(c(43,45,55,34,40,45,31,43,50,20,38,16,30,43,41,40,21,51,45,45,59,39,62,28,47,48,33,25,45,40,23,30,50,40,80,48,55,70,50,52,35,50,62,23,33,50,50,35,50,55,56,78,60,56,73,55,60,41,60,65,35,40,62,60,50,55,45,35),ncol=4)

iris=data.frame(mat,Species=c(rep(1,9),rep(2,8)))

ns=iris%>%group_by(Species)%>%summarise(n=n())

Speciess=ns$Species;k=length(Speciess)

p=ncol(iris)-1

D=array(0,dim=c(p,p-1))

for(j in 1:(p-1)){
  
D[j:(j+1),j]=c(1,-1)  
  
}

Sh=array(0,dim=c(p,p))

Se=array(0,dim=c(p,p))

n=nrow(iris);ns=ns$n;m=n-2

y_ave=apply(iris[,!(colnames(iris) %in% c("Species"))],2,mean)

y_ave_vec=array(0,dim=c(k,p))

for(j in 1:k){
  
iris_sub=iris%>%filter(Species==Speciess[j])%>%select(-Species)
 
y_ave_vec[j,]=apply(iris_sub,2,mean)

Sh=Sh+ns[j]*t(t(apply(iris_sub,2,mean)-y_ave))%*%t(apply(iris_sub,2,mean)-y_ave)
  
Se=Se+cov(iris_sub)*(nrow(iris_sub)-1)
  
}

lam1=det(t(D)%*%Se%*%D)/det(t(D)%*%(Se+Sh)%*%D)

test=-(m+(k-1-(p-1)-1)/2)*log(lam1)

#Trueなら有意差あり

test>qchisq(1-0.05,df=(p-1)*(k-1))




#平行幅の有意性

lam2=(det(Se)/det(Se+Sh))/lam1

test=(1-lam2)*(m-p+1)/(lam2*(k-1))

#Trueなら有意差あり（2群しかないので有意差は出ないはず）

test_statics>qf(1-0.05,p-1,m-p+2)


#平行性モデルのもとでの推測

Z=array(0,dim=c(k-1,p))

for(j in 1:(k-1)){
  
Z[j,]=y_ave_vec[j,]-y_ave_vec[k,]  
  
}

r=Z%*%solve(Se)%*%t(t(rep(1,p)))/as.numeric(t(rep(1,p))%*%solve(Se)%*%t(t(rep(1,p))))

sig=Se/m



```

```{r}

#p.83

mat=matrix(c(144,171,160,157,171,137,159,132,149,140,125,133,120,133,120,130,125,115,120,103,119,118,128,120,116,118,110,110),ncol=4)

Y=mat

n=nrow(Y)

m=n-1

S=cov(Y)

q=2

t=c(1,3,5,7)

p=length(t)

X=array(0,dim=c(p,q-1))

for(l in 1:(q-1)){
  
X[,l]=c(t^l)  
  
}

X=cbind(rep(1,p),X)

sita=solve(t(X)%*%solve(S)%*%X)%*%t(X)%*%solve(S)%*%apply(Y,2,mean)

var_sita=(m/(m-(p-q)))*((m-1)/(m-(p-q)-1))*solve(t(X)%*%solve(S)%*%X)/n

#sita=0の検定

D=c(0,1)

Se=m*t(D)%*%solve(t(X)%*%solve(S)%*%X)%*%D

t2pq=n*apply(Y,2,mean)%*%(solve(S)-solve(S)%*%X%*%solve(t(X)%*%solve(S)%*%X)%*%t(X)%*%solve(S))%*%apply(Y,2,mean)

R=as.numeric((1+t2pq/m)/n)

Sh=(t(D)%*%sita)%*%t(t(D)%*%sita)/R

lam=as.numeric(Se/(Se+Sh))

-((lam-1)/lam)*(m-(p-q))

qf(1-0.05,1,m-(p-q)+1-1)

#多項式モデルの適合性の検定

Se=m*S

Sh=n*t(t(c(apply(Y,2,mean))))%*%t(c(apply(Y,2,mean)))





```



```{r}


#p.93

mat1=matrix(c(144,171,160,157,171,137,159,132,149,140,125,133,120,133,120,130,125,115,120,103,119,118,128,120,116,118,110,110),ncol=4)

mat2=matrix(c(168,144,150,132,155,120,140,150,135,140,120,145,110,135,140,118,138,118,130,108,130,135,115,130,110,120,108,120),ncol=4)

n1=nrow(mat1);n2=nrow(mat2);k=2;n=n1+n2;m=n-k

t=c(1,3,5,7)

p=length(t)

q=2

X=array(0,dim=c(p,q-1))

for(l in 1:(q-1)){

X[,l]=c(t^l)  

}

X=cbind(rep(1,p),X)

S=(cov(mat1)*(n1-1)+cov(mat2)*(n2-1))/(n-k)

sita1=solve(t(X)%*%solve(S)%*%X)%*%t(X)%*%solve(S)%*%apply(mat1,2,mean)

sita2=solve(t(X)%*%solve(S)%*%X)%*%t(X)%*%solve(S)%*%apply(mat2,2,mean)

var_sita1=(m-1)*(m)*solve(t(X)%*%solve(S)%*%X)/(n1*(m-(p-q)-1)*(m-(p-q)))

var_sita2=(m-1)*(m)*solve(t(X)%*%solve(S)%*%X)/(n2*(m-(p-q)-1)*(m-(p-q)))

D=diag(rep(1,q))

C=diag(rep(1,k-1))

C=cbind(C,rep(-1,k-1))

sita_mat=t(cbind(sita1,sita2))

Se=t(D)%*%solve(t(X)%*%solve(S)%*%X)%*%D

y_vec_mat=rbind(apply(mat1,2,mean),apply(mat2,2,mean))

(C%*%sita_mat)%*%(t(X)%*%solve(S)%*%X)%*%t(C%*%sita_mat)

R=diag(1/c(n1,n2))+ y_vec_mat%*%(solve(S)-solve(S)%*%X%*%solve(t(X)%*%solve(S)%*%X)%*%t(X)%*%solve(S))%*%t(y_vec_mat)/m

lam2_value=(C%*%sita_mat)%*%(t(X)%*%solve(S)%*%X)%*%t(C%*%sita_mat)/as.numeric(m*C%*%R%*%t(C))

(m-(p-q))*lam2_value

qf(1-0.05,q,m-(p-q)+1-q)


```