経時データ解析　朝倉書店　船渡川伊久子・隆先生

```{r}

#p.39の例

group1=matrix(c(95,91,72,76,111,75,59,51,95,79,96,94,91,68,60,58,107,91,85,84,100,75,66,57,99,76,65,58,117,88,69,55),nrow=4)

group2=matrix(c(110,101,88,92,106,87,92,81,108,86,76,72,109,83,81,82,116,113,105,122,95,100,103,96,100,77,85,81,103,79,84,72),nrow=4)

#実薬群
K=prod(dim(group1))
p=4

X=diag(1,p)

loglik=function(param){
beta=param[1:p]
param=param[-c(1:p)]
sigG=param[1]
sigep=param[2]
V=array(sigG^2,dim=c(p,p))+diag(sigep^2,p)
value=-(K-p)*log(2*pi)/2-ncol(group1)*log(det(V))/2-ncol(group1)*log(det(X%*%solve(V)%*%X))/2
for(j in 1:ncol(group1)){
value=value-as.numeric(c(group1[,j]-X%*%beta)%*%solve(V)%*%c(group1[,j]-X%*%beta))}
return(value)
}

SIGG=1
SIGep=2
BETA=rep(1,p)

ite=10^3
h=0.001
eta=0.01

for(l in 1:ite){
  
Vmat=array(SIGG^2,dim=c(p,p))+diag(SIGep^2,p)  
XVX=array(0,dim=c(p,p))
XVY=rep(0,p)
for(j in 1:ncol(group1)){  
XVX=XVX+(X%*%Vmat%*%X)
XVY=XVY+X%*%Vmat%*%group1[,j]
}

BETA=BETA*(1-eta)+eta*solve(XVX)%*%XVY

parameter=c(BETA,SIGG,SIGep)
sub_parameter=parameter
sub_parameter[p+1]=sub_parameter[p+1]+h
SIGG=SIGG+eta*(loglik(sub_parameter)-loglik(parameter))/h

sub_parameter=parameter
sub_parameter[p+2]=sub_parameter[p+2]+h
SIGep=SIGep+eta*(loglik(sub_parameter)-loglik(parameter))/h

print(loglik(parameter))
}



b=c()
for(j in 1:ncol(group1)){
b=c(b,(SIGG^2)*rep(1,p)%*%solve(Vmat)%*%(group1[,j]-X%*%BETA))
}  

Yhat=array(0,dim=dim(group1))
for(i in 1:nrow(Yhat)){
for(j in 1:ncol(Yhat)){  
Yhat[i,j]=BETA[i]+b[j]  
}}

  
```

```{r}

#p.42の例1

group1=matrix(c(95,91,72,76,111,75,59,51,95,79,96,94,91,68,60,58,107,91,85,84,100,75,66,57,99,76,65,58,117,88,69,55),nrow=4)

group2=matrix(c(110,101,88,92,106,87,92,81,108,86,76,72,109,83,81,82,116,113,105,122,95,100,103,96,100,77,85,81,103,79,84,72),nrow=4)

#実薬群
K=prod(dim(group1))
p=4

X=diag(1,p)
Z=cbind(rep(1,p),c(0:3))



loglik=function(param){
beta=param[1:p]
param=param[-c(1:p)]
sigG=matrix(c(param[1:(ncol(Z)^2)]),nrow=ncol(Z))
sigep=param[-c(1:(ncol(Z)^2))]
V=Z%*%sigG%*%t(Z)+(sigep^2)*diag(rep(1,p))
value=-(K-p)*log(2*pi)/2-ncol(group1)*log(det(V))/2-ncol(group1)*log(det(X%*%solve(V)%*%X))/2
for(j in 1:ncol(group1)){
value=value-as.numeric(c(group1[,j]-X%*%beta)%*%solve(V)%*%c(group1[,j]-X%*%beta))}
return(value)
}

G=diag(10,ncol(Z))
SIGep=5
BETA=rep(1,p)

ite=10^4
h=0.001
eta=0.001

for(l in 1:ite){
  
Vmat=Z%*%G%*%t(Z)+(SIGep^2)*diag(1,p) 
XVX=array(0,dim=c(p,p))
XVY=rep(0,p)
dG=array(0,dim=c(ncol(Z),ncol(Z)))
for(j in 1:ncol(group1)){  
XVX=XVX+(X%*%Vmat%*%X)
XVY=XVY+X%*%Vmat%*%group1[,j]
dG=dG+t(c(group1[,j]-X%*%BETA)%*%solve(Vmat)%*%solve(Vmat)%*%Z)%*%t(t(Z)%*%(c(group1[,j]-X%*%BETA)))
}

BETA=BETA*(1-eta)+eta*solve(XVX)%*%XVY

dG=dG-ncol(group1)*t(t(solve(t(X)%*%(Vmat)%*%X))%*%t(X)%*%Z)%*%t(t(Z)%*%X)-ncol(group1)*t(solve(Vmat)%*%Z)%*%(Z)

G=G+eta*dG

sub_parameter=parameter
sub_parameter[length(sub_parameter)]=sub_parameter[length(sub_parameter)]+h
SIGep=SIGep+eta*(loglik(sub_parameter)-loglik(parameter))/h

parameter=(c(BETA,c(G),SIGep))
print(loglik(parameter))
}



b=c()
for(j in 1:ncol(group1)){
b=c(b,G%*%t(Z)%*%solve(Vmat)%*%(group1[,j]-X%*%BETA))
}  
b=matrix(b,nrow=ncol(Z))

Yhat=array(0,dim=dim(group1))
for(j in 1:ncol(Yhat)){
Yhat[,j]=X%*%BETA+Z%*%b[,j]
}

  
```

```{r}

#p.42の例2

group1=matrix(c(95,91,72,76,111,75,59,51,95,79,96,94,91,68,60,58,107,91,85,84,100,75,66,57,99,76,65,58,117,88,69,55),nrow=4)

group2=matrix(c(110,101,88,92,106,87,92,81,108,86,76,72,109,83,81,82,116,113,105,122,95,100,103,96,100,77,85,81,103,79,84,72),nrow=4)

#実薬群
K=prod(dim(group1))

p=2
Z=cbind(rep(1,p),c(0:3))

loglik=function(param){
beta=param[1:p]
param=param[-c(1:p)]
sigG=matrix(c(param[1:(ncol(Z)^2)]),nrow=ncol(Z))
sigep=param[-c(1:(ncol(Z)^2))]
V=Z%*%sigG%*%t(Z)+(sigep^2)*diag(rep(1,4))
value=-(K-p)*log(2*pi)/2-ncol(group1)*log(det(V))/2-ncol(group1)*log(det(t(Z)%*%solve(V)%*%Z))/2
for(j in 1:ncol(group1)){
value=value-as.numeric(c(group1[,j]-Z%*%beta)%*%solve(V)%*%c(group1[,j]-Z%*%beta))}
return(value)
}

G=diag(10,ncol(Z))
SIGep=5
BETA=rep(1,p)

ite=10^4
h=0.001
eta=0.001

for(l in 1:ite){
  
Vmat=Z%*%G%*%t(Z)+(SIGep^2)*diag(1,4) 
XVX=array(0,dim=c(p,p))
XVY=rep(0,p)
dG=array(0,dim=c(ncol(Z),ncol(Z)))
for(j in 1:ncol(group1)){  
XVX=XVX+(t(Z)%*%Vmat%*%Z)
XVY=XVY+t(Z)%*%Vmat%*%group1[,j]
dG=dG+t(c(group1[,j]-Z%*%BETA)%*%solve(Vmat)%*%solve(Vmat)%*%Z)%*%t(t(Z)%*%(c(group1[,j]-Z%*%BETA)))
}

BETA=BETA*(1-eta)+eta*solve(XVX)%*%XVY

dG=dG-ncol(group1)*t(t(solve(t(Z)%*%(Vmat)%*%Z))%*%t(Z)%*%Z)%*%t(t(Z)%*%Z)-ncol(group1)*t(solve(Vmat)%*%Z)%*%(Z)

G=G+eta*dG

parameter=(c(BETA,c(G),SIGep))
sub_parameter=parameter
sub_parameter[length(sub_parameter)]=sub_parameter[length(sub_parameter)]+h
SIGep=SIGep+eta*(loglik(sub_parameter)-loglik(parameter))/h


print(loglik(parameter))
}



b=c()
for(j in 1:ncol(group1)){
b=c(b,G%*%t(Z)%*%solve(Vmat)%*%(group1[,j]-Z%*%BETA))
}  
b=matrix(b,nrow=ncol(Z))

Yhat=array(0,dim=dim(group1))
for(j in 1:ncol(Yhat)){
Yhat[,j]=Z%*%BETA+Z%*%b[,j]
}

  
```


```{r}


#p.42の例(p.32参照)

#非線形混合変量効果モデル  

#Taylor展開1次までで展開し、非線形関数を線形の形に変更(薬物動態の統計数理　医学統計ライブラリー　後藤・大門先生 p.101)


group1=matrix(c(95,91,72,76,111,75,59,51,95,79,96,94,91,68,60,58,107,91,85,84,100,75,66,57,99,76,65,58,117,88,69,55),nrow=4)

group2=matrix(c(110,101,88,92,106,87,92,81,108,86,76,72,109,83,81,82,116,113,105,122,95,100,103,96,100,77,85,81,103,79,84,72),nrow=4)

#実薬群
K=prod(dim(group1))
p=4


X=diag(1,p)

#非線形関数の設定　
#exp(X%*%beta)がいい
f=function(beta){
return(exp(X%*%beta)+X%*%beta)  
}

loglik=function(param){
beta=param[1:p]
Z=array(0,dim=c(nrow(X),p))  
for(j in 1:p){
beta2=beta
beta2[j]=beta2[j]+h
Z[,j]=(f(beta2)-f(beta))/h
}
param=param[-c(1:p)]
sigG=matrix(c(param[1:(p^2)]),nrow=p)
sigep=param[-c(1:(p^2))]
V=Z%*%sigG%*%t(Z)+(sigep^2)*diag(rep(1,p))
value=-(K-p)*log(2*pi)/2-ncol(group1)*log(det(V))/2-ncol(group1)*log(det(X%*%solve(V)%*%X))/2
for(j in 1:ncol(group1)){
value=value-as.numeric(c(group1[,j]-f(beta))%*%solve(V)%*%c(group1[,j]-f(beta)))}
return(value)
}

G=diag(10,p)
SIGep=5
BETA=rep(0.0001,p)

ite=2*10^4
h=0.001
eta=0.0001

for(l in 1:ite){
Z=array(0,dim=c(nrow(X),p))  
for(j in 1:p){
BETA2=BETA
BETA2[j]=BETA2[j]+h
Z[,j]=(f(BETA2)-f(BETA))/h
}
Vmat=Z%*%G%*%t(Z)+(SIGep^2)*diag(1,p) 
dG=array(0,dim=c(p,p))
XVX=array(0,dim=c(p,p))
XVY=rep(0,p)
dBETA=rep(0,length(BETA))
for(j in 1:ncol(group1)){  
#dBETA=dBETA+2*(X*c(f(BETA)))%*%solve(Vmat)%*%(group1[,j]-f(BETA))
dBETA=dBETA+2*(Z)%*%solve(Vmat)%*%(group1[,j]-f(BETA))
dG=dG+t(c(group1[,j]-f(BETA))%*%solve(Vmat)%*%solve(Vmat)%*%Z)%*%t(t(Z)%*%(c(group1[,j]-f(BETA))))
}

BETA=BETA+eta*dBETA

dG=dG-ncol(group1)*t(t(solve(t(X)%*%(Vmat)%*%X))%*%t(X)%*%Z)%*%t(t(Z)%*%X)-ncol(group1)*t(solve(Vmat)%*%Z)%*%(Z)

G=G+eta*dG
parameter=(c(BETA,c(G),SIGep))
sub_parameter=parameter
sub_parameter[length(sub_parameter)]=sub_parameter[length(sub_parameter)]+h
SIGep=SIGep+eta*(loglik(sub_parameter)-loglik(parameter))/h

parameter=(c(BETA,c(G),SIGep))
print(loglik(parameter))
}



b=c()
for(j in 1:ncol(group1)){
b=c(b,G%*%t(Z)%*%solve(Vmat)%*%(group1[,j]-f(BETA)))
}  
b=matrix(b,nrow=ncol(Z))

#変量パラメータ=0の周りでの1次までのTaylor展開を用いた場合
Yhat1=array(0,dim=dim(group1))
for(j in 1:ncol(Yhat1)){
Yhat1[,j]=f(BETA)+Z%*%b[,j]
}

#もともとの非線形関数での予測
Yhat2=array(0,dim=dim(group1))
for(j in 1:ncol(Yhat2)){
Yhat2[,j]=f(BETA+b[,j])
}

  
```