回帰モデルの種類と検定

```{r}
#回帰モデルの種類と検定

library(dplyr)


data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

cor(cbind(data$y,data$x2))

beta=cov(data$y,data$x2)/cov(data$x2,data$x2)

beta0=mean(data$y)-beta*mean(data$x2)

e_vec=(data$y-beta0-beta*data$x2)^2



value=data$x2
dim=2
vec=ans=rep(1,nrow(data))

for(j in 1:dim){
  vec=vec*value
  ans=cbind(ans,vec)

}
X=t(matrix(ans,ncol=ncol(ans),nrow=nrow(ans)))
a=(solve(X%*%t(X)))%*%X%*%data$y
y=t(X)%*%a
result=data.frame(num=1:nrow(data),y=data$y,predict=y)
V=nrow(result)
result=result%>% dplyr::mutate(error=(y-predict)^2/V)

plot(result$num,result$y,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=2)
par(new=T)
plot(result$num,result$predict,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=3)

print(paste0("sum of errors is ",sum(result$error)))

R2=(sum((result$predict-mean(result$y))^2)/sum((result$y-mean(result$y))^2))

R=sqrt(R2)

n=nrow(data)

kai2=n*R2

alpha=0.01

#white test

#H:r1=r2=0

if(kai2>qchisq(1-alpha,df=2)){
  
print("分散は不均一")  
  
}else{
  
print("分散は均一")   
  
}

HCSE=sqrt(sum(e_vec*(data$x2-mean(data$x2))^2)/sum((data$x2-mean(data$x2))^2)^2)

V_b=HCSE^2

```

```{r}

#H:beta1+beta2=1

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=as.numeric(data$y)

alpha=cov(data$y-data$x2,data$x1-data$x2)/cov(data$x1-data$x2,data$x1-data$x2)

alpha0=mean(data$y-data$x2)-alpha*mean(data$x1-data$x2)

RSSr=sum((Y-data$x2-alpha0-(data$x1-data$x2)*alpha)^2)

RSSu=sum((Y-data$x2-mean(Y-data$x2))^2)

f=(RSSr-RSSu)/(RSSu/(nrow(data)-3))

qf(1-0.01,1,nrow(data)-3)



```

```{r}

#serial correlations(cochrane-0rcutt method)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=as.numeric(data$y)

alpha=solve(t(X)%*%X)%*%t(X)%*%Y

alpha0=mean(Y)-sum(alpha*apply(X,2,mean))

e_vec=Y-alpha0-X%*%alpha

p=sum(e_vec[1:(length(e_vec)-1)]*e_vec[2:length(e_vec)])/sum(e_vec^2)

serial_error=data.frame(e_value=e_vec,predict=0)

for(j in 2:nrow(serial_error)){
  
serial_error$predict[j]=p*e_vec[j-1]

}




plot(c(1:nrow(serial_error)),serial_error$e_value,type="o",col=2,xlab="num",ylab="values",xlim=c(1,nrow(serial_error)),ylim=c(min(serial_error$e_value,serial_error$predict),max(serial_error$e_value,serial_error$predict)))

par(new=T)

plot(c(1:nrow(serial_error)),serial_error$predict,type="o",col=3,xlab="num",ylab="values",xlim=c(1,nrow(serial_error)),ylim=c(min(serial_error$e_value,serial_error$predict),max(serial_error$e_value,serial_error$predict)))



data=data %>% mutate(predict=alpha0+X%*%alpha,time_series=alpha0+X%*%alpha+serial_error$predict)

cor(cbind(data$y,data$predict,data$time_series))



DW=sum((e_vec[2:length(e_vec)]-e_vec[1:(length(e_vec)-1)])^2)/sum(e_vec^2)


```





```{r}

#一般化線形モデルを使用した損害保険料率の算定　日本アクチュアリー会

#Exponential dispersion model(EDM)


#relative pois model

library(dplyr)

data=data.frame(cell=c(1:17),x1=c(rep(1,9),rep(2,8)),x2=c(rep(1,7),2,2,1,2,2,2,2,2,2,2),x3=c(1,2,3,4,5,6,7,1,2,7,1,2,3,4,5,6,7),w=c(62.9,112.9,133.1,376.6,9.4,70.8,4.4,352.1,840.1,14.5,844.8,1296,1214.9,3740.7,109.4,404.7,66.3),y=c(270,62,68,19,0,14,228,148,82,0,111,76,30,15,37,12,15))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=c(data$y);W=c(data$w);

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

ite=50

beta=rep(1,ncol(X))

for(j in 1:ite){
  
mu=exp(X%*%beta)  
  
W_mat=diag(c(W*mu))

U=(W*t(Y*log(mu)-mu))%*%X

XtWX <- t(X) %*% W_mat %*% X

beta <- beta + solve(XtWX) %*% t(U)
    
#print(beta)

C=c()

for(k in 1:nrow(X)){
  
if(Y[k]>0){  
  
C=c(C,W[k]*Y[k]*log(W[k])-sum(log(c(1:(W[k]*Y[k])))))  
  
}else{
  
C=c(C,W[k]*Y[k]*log(W[k]))  
  
}
  
}  

print(sum(W*(Y*log(mu)-mu))+sum(C))

}

plot(c(1:length(Y)),Y*log(mu),col=2,type="o",xlab="index",ylab="values",ylim=c(min(Y*log(mu)),max(Y*log(mu))))

par(new=T)

plot(c(1:length(Y)),mu,col=3,type="o",xlab="index",ylab="values",ylim=c(min(Y*log(mu)),max(Y*log(mu))))

```

```{r}

#Exponential dispersion model(EDM) p.１５１


#乗法型ガンマモデル

library(dplyr)

data=data.frame(cell=c(1:17),x1=c(rep(1,9),rep(2,8)),x2=c(rep(1,7),2,2,1,2,2,2,2,2,2,2),x3=c(1,2,3,4,5,6,7,1,2,7,1,2,3,4,5,6,7),w=c(62.9,112.9,133.1,376.6,9.4,70.8,4.4,352.1,840.1,14.5,844.8,1296,1214.9,3740.7,109.4,404.7,66.3),y=c(270,62,68,19,10,14,228,148,82,10,111,76,30,15,37,12,15))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

#B-spline transform

X=exp(-X)

Y=c(data$y);W=c(data$w);

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

ite=50

lam=1

beta=rep(1,ncol(X))

for(j in 1:ite){
  
beta_pre=beta  

mu=exp(X%*%beta)  

W_mat=diag(c(W*Y/mu))

U=1-mu/Y+X%*%beta

XtWX <- t(X) %*% W_mat %*% X

beta<-solve(XtWX+diag(rep(lam,ncol(X))))%*%t(X)%*%W_mat%*%(U)

print(sum(abs(beta_pre-beta)))

}

plot(c(1:length(Y)),log(Y/mu),col=2,type="o",xlab="index",ylab="values",ylim=c(min(log(Y/mu)),max(log(Y/mu))))

par(new=T)

plot(c(1:length(Y)),Y/mu-1,col=3,type="o",xlab="index",ylab="values",ylim=c(min(log(Y/mu)),max(log(Y/mu))))


```

```{r}

#tweedie model

library(dplyr)

data=data.frame(cell=c(1:17),x1=c(rep(1,9),rep(2,8)),x2=c(rep(1,7),2,2,1,2,2,2,2,2,2,2),x3=c(1,2,3,4,5,6,7,1,2,7,1,2,3,4,5,6,7),w=c(62.9,112.9,133.1,376.6,9.4,70.8,4.4,352.1,840.1,14.5,844.8,1296,1214.9,3740.7,109.4,404.7,66.3),y=c(270,62,68,19,0,14,228,148,82,0,111,76,30,15,37,12,15))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=c(data$y);W=c(data$w);W=abs(W-mean(W))/sd(W)

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

ite=50

beta=rep(1,ncol(X));pthi=3;p=2

for(j in 1:ite){
  
mu=exp(X%*%beta)  
  
H=-t(X)%*%diag(c(W*mu^p)/pthi)%*%X  

W_mat=diag(c(W*mu))
  
beta=beta-solve(H)%*%(t(X)%*%W_mat%*%Y-t(X)%*%W_mat%*%mu)  
  
}

plot(c(1:length(Y)),Y,col=2,type="o",xlab="index",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

plot(c(1:length(Y)),mu,col=3,type="o",xlab="index",ylab="values",ylim=c(min(Y),max(Y)))

```

```{r}

#normal model

library(dplyr)

data=data.frame(cell=c(1:17),x1=c(rep(1,9),rep(2,8)),x2=c(rep(1,7),2,2,1,2,2,2,2,2,2,2),x3=c(1,2,3,4,5,6,7,1,2,7,1,2,3,4,5,6,7),w=c(62.9,112.9,133.1,376.6,9.4,70.8,4.4,352.1,840.1,14.5,844.8,1296,1214.9,3740.7,109.4,404.7,66.3),y=c(270,62,68,19,0,14,228,148,82,0,111,76,30,15,37,12,15))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=c(data$y);W=c(data$w);W=abs(W-mean(W))/sd(W)

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

ite=500

beta=rep(1,ncol(X))

var_Y=1

for(j in 1:ite){

mu=(X%*%beta)  

W_mat=diag(c(var_Y/W))

U=((W/var_Y)*t((Y-mu)))%*%X

XtWX <- t(X) %*% W_mat %*% X

beta <- beta + solve(XtWX) %*% t(U)

#print(beta)

print((-1/2)*sum((W/var_Y)*(Y-mu)^2)-sum(log(2*pi*var_Y/W))/2)

}

plot(c(1:length(Y)),Y,col=2,type="o",xlab="index",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

plot(c(1:length(Y)),mu,col=3,type="o",xlab="index",ylab="values",ylim=c(min(Y),max(Y)))

```
```{r}

#p.42

library(dplyr)

m1=10;m2=20

Y=matrix(sample(0,m1*m2,replace = T),ncol=m1)

W=matrix(sample(0,m1*m2,replace = T),ncol=m1)


for(j in 1:m1){
  
Y[,j]=rpois(m2,sample(10:100,1))

W[,j]=rpois(m2,sample(30:50,1))
  
}

for(i in 1:m2){
  
Y[i,]=Y[i,]+rpois(m1,sample(10:100,1))

W[i,]=W[i,]+rpois(m1,sample(30:50,1))
  
}

r0=1;r_1=rep(1,m1);r_2=rep(1,m2)

ite=10

for(k in 1:ite){
  
r_1_pre=r_1  
  
r0=sum(W*Y)/sum(t((W)*r_2)*r_1) 
  
r_1=c(apply((W*Y),2,sum)/(r0*t(W)%*%r_2))

r_2=c(apply((W*Y),1,sum)/(r0*W%*%r_1))

log_lik=sum(W*(Y*log(r0*t(t(r_2))%*%t(r_1))-r0*t(t(r_2))%*%t(r_1)))

print(paste0("log_lik:",log_lik))


}


```

```{r}

#線形混合モデル

n=100;p1=15;p2=19

Y=sample(1:100,n,replace=T)

X=array(0,dim=c(n,p1))

Z=array(0,dim=c(n,p2))

for(j in 1:n){

  X[j,]=rpois(p1,sample(1:20,1))  
  
  Z[j,]=rpois(p2,sample(20:80,1))
  
}

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

sigma_s=15

sigma_c=c(11,13,16,18)

D1=diag(sigma_s,15);D2=diag(sigma_c)

Z1=Z[,1:15];Z2=Z[,16:19]

V=Z1%*%D1%*%t(Z1)+Z2%*%D2%*%t(Z2)+diag(1,n)

beta=solve(t(X)%*%solve(V)%*%X)%*%t(X)%*%solve(V)%*%Y

D=diag(c(diag(D1),diag(D2)))

u=D%*%t(Z)%*%solve(V)%*%(Y-X%*%beta)

predict=X%*%beta+Z%*%u


plot(c(1:n),Y,type="o",col=2,xlim=c(1,n),ylim=c(min(c(predict,Y)),max(c(Y,predict))),xlab="index",ylab="values")

par(new=T)

plot(c(1:n),predict,type="o",col=3,xlim=c(1,n),ylim=c(min(c(predict,Y)),max(c(predict,Y))),xlab="index",ylab="values")

res=data.frame(no=c(1:n),Y,predict) %>% mutate(error=abs(Y-predict))

cor(predict,Y)

```


```{r}
#nonparametric regression methods for logitudal data analysis

#Penalized spline method p.234

library(dplyr)

data=data.frame(longley)/1000



x=as.matrix(data[,colnames(data) %in% c("GNP.deflator","Unemployed","Employed")]);Y=data[,colnames(data) %in% "GNP"]

c=as.matrix(data[,colnames(data) %in% c("Armed.Forces","Population")])


train_no=c(1:floor(nrow(x)))



c_train=c[train_no,];c_test=c[-train_no,]

x_train=x[train_no,];x_test=x[-train_no,]

y_train=Y[train_no];y_test=Y[-train_no]


lambda_dat=data.frame(lambda=seq(0.1,10000,0.1),GCV=0)

for(j in 1:nrow(lambda_dat)){
  
lambda=lambda_dat$lambda[j]

P=diag(1,nrow(x_train))-(x_train%*%solve(t(x_train)%*%x_train+diag(lambda,ncol(x_train))))%*%t(x_train)
  
Q=diag(1,nrow(x_train))-c_train%*%solve(t(c_train)%*%c_train+diag(10^(-2),ncol(c_train)))%*%t(c_train)

mat=(t(c_train)%*%P%*%c_train);inv_mat=svd(mat)$u%*%diag(1/svd(mat)$d)%*%t(svd(mat)$v)

alpha=inv_mat%*%t(c_train)%*%P%*%y_train

beta=solve(t(x_train)%*%Q%*%x_train+diag(lambda,ncol(x_train)))%*%t(x_train)%*%y_train

y_hat=c_train%*%alpha+x_train%*%beta

A_par=c_train%*%solve(t(c_train)%*%c_train+diag(10^(-2),ncol(c_train)))%*%t(c_train)%*%P

A_npar=(x_train%*%solve(t(x_train)%*%x_train+diag(lambda,ncol(x_train))))%*%t(x_train)%*%Q

tr_A=c(diag(A_par+A_npar))

error=t(c(y_train-y_hat))

GCV=error%*%t(error)/((1-mean(tr_A))^2)

lambda_dat$GCV[j]=GCV

}

lambda=lambda_dat$lambda[lambda_dat$GCV==min(lambda_dat$GCV)]



plot(lambda_dat$lambda,lambda_dat$GCV)

P=diag(1,nrow(x_train))-(x_train%*%solve(t(x_train)%*%x_train+diag(lambda,ncol(x_train))))%*%t(x_train)
  
Q=diag(1,nrow(x_train))-c_train%*%solve(t(c_train)%*%c_train+diag(10^(-2),ncol(c_train)))%*%t(c_train)

mat=(t(c_train)%*%P%*%c_train);inv_mat=svd(mat)$u%*%diag(1/svd(mat)$d)%*%t(svd(mat)$v)

alpha=inv_mat%*%t(c_train)%*%P%*%y_train

beta=solve(t(x_train)%*%Q%*%x_train+diag(lambda,ncol(x_train)))%*%t(x_train)%*%y_train

y_hat=c_train%*%alpha+x_train%*%beta

A_par=c_train%*%solve(t(c_train)%*%c_train+diag(10^(-2),ncol(c_train)))%*%t(c_train)%*%P

A_npar=(x_train%*%solve(t(x_train)%*%x_train+diag(lambda,ncol(x_train))))%*%t(x_train)%*%Q

tr_A=c(diag(A_par),diag(A_npar))

error=t(c(y_train-y_hat))

GCV=error%*%t(error)



plot(c(1:length(y_hat)),y_hat,xlim=c(1,length(y_hat)),ylim=c(min(c(y_hat,y_train)),max(c(y_hat,y_train))),xlab="index",ylab="values",col=2)

par(new=T)

plot(c(1:length(y_hat)),y_train,xlim=c(1,length(y_hat)),ylim=c(min(c(y_hat,y_train)),max(c(y_hat,y_train))),xlab="index",ylab="values",col=3)



```

```{r}

#nonparametric regression methods for logitudal data analysis p.237

#semiparametric population mean model(smoothing spline method)

library(dplyr)


data=data.frame(USArrests)

colnames(data)=c("y","x1","x2","x3")

data=data.frame(t=c(1:nrow(data)),data)

N=nrow(data)


#data=data.frame(num=1:n,y=sample(c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),n,replace=T),x1=sample(c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),n,replace=T),x2=sample(c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),n,replace=T))

n=floor(nrow(data)*0.7)

U=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Z=data$y

kernel=function(x,y){
  
return(exp(-abs(x-y)))  
  
}

S=matrix(1,nrow=nrow(data),ncol=nrow(data))

for(j in 1:ncol(S)){
for(i in 1:nrow(S)){  
  
S[j,i]=kernel(data$t[i],data$t[j]) 
   
}
}  

S=S/sum(S)

delta=solve(t(U)%*%(diag(1,N)-S)%*%U)%*%t(U)%*%(diag(1,N)-S)%*%Z

f=S%*%(Z-U%*%delta)

result=data.frame(num=1:nrow(data),y=Z,predict=U%*%delta+f)

plot(result$num,result$y,xlim=c(1,nrow(result)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=2,main=paste0("cor:",cor(result$y,result$predict)))
par(new=T)
plot(result$num,result$predict,xlim=c(1,nrow(result)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=3,main=paste0("cor:",cor(result$y,result$predict)))


```



```{r}

#線形混合モデル2 一般化

#non-parametric regression method for logitudal data analysis p.24

n=100;p1=15;p2=19;

class=3

N=class*n





XX=array(0,dim=c(1,p1+1))

ZZ=array(0,dim=c(1,p2+1))

YY=array(0,dim=c(1,2))


for(i in 1:class){
  
X=array(0,dim=c(n,p1))

Z=array(0,dim=c(n,p2))

for(j in 1:n){

  X[j,]=rpois(p1,sample(1:20,1))  

  Z[j,]=rpois(p2,sample(20:80,1))

  }
  
X=cbind(rep(i,n),X);Z=cbind(rep(i,n),Z);Y=cbind(rep(i,n),sample(1:100,n,replace=T))

XX=rbind(XX,X);ZZ=rbind(ZZ,Z);YY=rbind(YY,Y)

}      

XX=tail(XX,nrow(XX)-1);ZZ=tail(ZZ,nrow(ZZ)-1)

YY=tail(YY,nrow(YY)-1)

XX=data.frame(XX);ZZ=data.frame(ZZ);YY=data.frame(YY)

sigma=1;beta=rep(1,p1)

V_box=array(0,dim=c(n,N))

D=diag(1,p2)

sigma=1;

ite=10000

for(i in 1:ite){
  
V_box=array(0,dim=c(n,1))  

beta_sub1=array(0,dim=c(p1,p1))

beta_sub2=array(0,dim=c(p1,1))

b_vec=array(0,dim=c(p2,class))

sigma_vec=0;D_mat=array(0,dim=c(p2,p2))
  
for(j in 1:(class)){ 
  
Z=ZZ %>% filter(X1==j) %>% select(-X1)

X=XX %>% filter(X1==j) %>% select(-X1)

Y=YY %>% filter(X1==j) %>% select(-X1)

Z=as.matrix(Z);X=as.matrix(X);Y=as.matrix(Y)

V=Z%*%D%*%t(Z)+diag(sigma,n)

V_box=cbind(V_box,V)

beta_sub1=beta_sub1+t(X)%*%solve(V)%*%X

beta_sub2=beta_sub2+t(X)%*%solve(V)%*%Y

#beta=solve(t(X1)%*%solve(V1)%*%X1+t(X2)%*%solve(V2)%*%X2)%*%(t(X1)%*%solve(V1)%*%Y1+t(X2)%*%solve(V2)%*%Y2)  

b=D%*%t(Z)%*%solve(V)%*%(Y-X%*%beta)

b_vec[,j]=b
 
epsilon=Y-X%*%beta-Z%*%b;

sigma_vec=sigma_vec+as.numeric(t(epsilon)%*%epsilon+sigma*(n-sigma*sum(diag(solve(V)))))/N

#sigma=as.numeric(t(epsilon1)%*%epsilon1+t(epsilon2)%*%epsilon2+sigma*(n-sigma*sum(diag(solve(V1))))+sigma*(n-sigma*sum(diag(solve(V2)))))/N

D_mat=D_mat+(b%*%t(b)+(D-D%*%t(Z)%*%solve(V)%*%Z%*%D))/class

#D=((b1)%*%t(b1)+(b2)%*%t(b2)+(D-D%*%t(Z1)%*%solve(V1)%*%Z1%*%D)+(D-D%*%t(Z2)%*%solve(V2)%*%Z2%*%D) )/2  

#R1=diag(sigma,n);R2=diag(sigma,n)

#GLL=t(Y1-X1%*%beta-Z1%*%b1)%*%solve(R1)%*%(Y1-X1%*%beta-Z1%*%b1)+t(Y2-X2%*%beta-Z2%*%b2)%*%solve(R2)%*%(Y2-X2%*%beta-Z2%*%b2)+t(b2)%*%solve(D)%*%b2+t(b1)%*%solve(D)%*%b1+2*log(det(D))+log(det(R1))+log(det(R2))



}  
  
beta=solve(beta_sub1)%*%beta_sub2

sigma=sigma_vec;D=D_mat

GLL=0

for(j in 1:(class)){ 
  
Z=ZZ %>% filter(X1==j) %>% select(-X1)

X=XX %>% filter(X1==j) %>% select(-X1)

Y=YY %>% filter(X1==j) %>% select(-X1)

Z=as.matrix(Z);X=as.matrix(X);Y=as.matrix(Y)

b=b_vec[,j]

GLL=GLL+as.numeric(t(Y-X%*%beta-Z%*%b)%*%diag(1/sigma,n)%*%(Y-X%*%beta-Z%*%b)+t(b)%*%solve(D)%*%b+log(det(D))+log(det(diag(sigma,n))))

}

print(GLL)  

}


```

```{r}

#non-parametric regression method for logitudal data analysis

#Naive Local Polynomial kernel Method p.73

library(dplyr)

t=3

n=10;N=2

class=5

mat=array(0,dim=c(1,N+2))

Y=c()

t_vec=c()

for(i in 1:class){
  

ts=sample(c((t+1):100),n,replace=T)

t_vec=c(t_vec,ts-t)

t_mat=array(0,dim=c(n,N))

for(j in 1:nrow(t_mat)){
  
t_mat[j,]=(ts[j]-t)^c(1:N)  
  

}

t_mat=cbind(rep(i,n),rep(1,n),t_mat)

mat=rbind(mat,t_mat)

Y=c(Y,apply(t(t_mat)*c(sample(1:1000,n,replace=T)),2,sum))

}

mat=tail(mat,nrow(mat)-1)

K=diag(c(exp(-abs(t_vec))/sum(exp(-abs(t_vec)))))

X=mat[,3:ncol(mat)]


beta=solve(t(X)%*%K%*%X)%*%t(X)%*%K%*%Y

sum(c(Y-X%*%beta)*diag(K))

res=data.frame(Y=Y,predict=X%*%beta)

est.Y=as.numeric(X%*%beta)

plot(c(1:length(Y)),Y,col=2,xlim=c(1,length(Y)),ylim=c(min(c(Y,est.Y)),max(c(Y,est.Y))),xlab="index",ylab="values",type="o")

par(new=T)

plot(c(1:length(Y)),est.Y,col=3,xlim=c(1,length(Y)),ylim=c(min(c(Y,est.Y)),max(c(Y,est.Y))),xlab="index",ylab="values",type="o")

```
```{r}

#non-parametric regression method for logitudal data analysis
#p.193


library(dplyr)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=as.numeric(data$y);N=length(Y);X=cbind(rep(1,N),X)

G=diag(1,ncol(X))



GCV=function(s){
  
y_hat=X%*%solve(t(X)%*%X+s*G)%*%t(X)%*%Y  

df=t(c(y_hat-Y))

W=diag(1/N,N)
  
GCV=df%*%W%*%t(df)/((1-sum(diag(solve(t(X)%*%X+s*G)%*%t(X)%*%X))/N)^2 ) 
  
return(GCV)

}

GCV_data=data.frame(lambda=seq(450,550,0.001),GCV=0)

for(j in 1:nrow(GCV_data)){
  
GCV_data$GCV[j]=as.numeric(GCV(as.numeric(GCV_data$lambda[j])) ) 
  
}

plot(GCV_data$lambda,GCV_data$GCV)


lambda=GCV_data$lambda[GCV_data$GCV==min(GCV_data$GCV)]

beta=solve(t(X)%*%X+lambda*G)%*%t(X)%*%Y

plot(c(1:N),Y,xlim=c(1,N),ylim=c(min(c(Y,X%*%beta)),max(c(Y,X%*%beta))),type="o",col=2,xlab="index",ylab="values",main=paste0("cor:",cor(Y,c(X%*%beta))))

par(new=T)

plot(c(1:N),c(X%*%beta),xlim=c(1,N),ylim=c(min(c(Y,X%*%beta)),max(c(Y,X%*%beta))),type="o",col=3,xlab="index",ylab="values",main=paste0("cor:",cor(Y,c(X%*%beta))))


t_eta_data=data.frame(t=seq(1,50,1),eta=0,V_eta=0,under=0,upper=0)

sigma_hat=sum((Y-c(X%*%beta))^2)/(N-sum(diag(solve(t(X)%*%X+lambda*G)%*%t(X)%*%X)))


for(j in 1:nrow(t_eta_data)){

t=t_eta_data$t[j]

pthi=t^(c(1:ncol(X))-1)

eta=as.numeric(t(pthi)%*%solve(t(X)%*%X+lambda*G)%*%t(X)%*%Y)

t_eta_data$eta[j]=eta

var_eta=as.numeric(as.numeric(sigma_hat)*t(pthi)%*%solve(t(X)%*%X+lambda*G)%*%t(X)%*%X%*%solve(t(X)%*%X+lambda*G)%*%pthi)

t_eta_data$V_eta[j]=var_eta

t_eta_data$under[j]=eta-qnorm(1-0.05/2)*sqrt(var_eta)

t_eta_data$upper[j]=eta+qnorm(1-0.05/2)*sqrt(var_eta)

}

plot(c(t_eta_data$t),c(t_eta_data$eta),xlim=c(min(t_eta_data$t),max(t_eta_data$t)),ylim=c(min(t_eta_data$under),max(t_eta_data$upper)),xlab="index",ylab="values",type="o",col=2)

par(new=T)

plot(c(t_eta_data$t),c(t_eta_data$under),xlim=c(min(t_eta_data$t),max(t_eta_data$t)),ylim=c(min(t_eta_data$under),max(t_eta_data$upper)),xlab="index",ylab="values",type="o",col=3)

par(new=T)

plot(c(t_eta_data$t),c(t_eta_data$upper),xlim=c(min(t_eta_data$t),max(t_eta_data$t)),ylim=c(min(t_eta_data$under),max(t_eta_data$upper)),xlab="index",ylab="values",type="o",col=3,main="eta_nps")


```
```{r}


#線形混合モデル3 一般化

#non-parametric regression method for logitudal data analysis

#p.170 REML-Based EM-Algorithm

n=1000;p1=30;

class=9

N=class*n





XX=array(0,dim=c(1,p1+1))

YY=array(0,dim=c(1,2))


for(i in 1:class){

X=array(0,dim=c(n,p1))

for(j in 1:n){

  X[j,]=rpois(p1,sample(1:20,1))  

  }

X=cbind(rep(i,n),X);Y=cbind(rep(i,n),sample(1:100,n,replace=T))

XX=rbind(XX,X);YY=rbind(YY,Y)

}      

XX=tail(XX,nrow(XX)-1);YY=tail(YY,nrow(YY)-1)

XX=data.frame(XX);YY=data.frame(YY)

sigma=1;eta=rep(1,p1)

V_box=array(0,dim=c(n,N))

D=diag(1,p1);


G=diag(1,p1);lambda=1;lambda_v=1



ite=14

for(i in 1:ite){
  
eta_pre=eta
  
D=solve(solve(D)+lambda_v*G)

V_box=array(0,dim=c(n,1))

v_vec=array(0,dim=c(p1,class))

eta_sub1=array(0,dim=c(p1,p1))

eta_sub2=array(0,dim=c(p1,1))

for(j in 1:class){
  
X=XX %>% filter(X1==j) %>% select(-X1)

Y=YY %>% filter(X1==j) %>% select(-X1)

X=as.matrix(X);Y=as.matrix(Y)
  
V=X%*%D%*%t(X)+diag(sigma,n)

V_box=cbind(V_box,V)

eta_sub1=eta_sub1+t(X)%*%solve(V)%*%X

eta_sub2=eta_sub2+t(X)%*%solve(V)%*%Y

} 

eta=solve(eta_sub1+lambda*G)%*%eta_sub2

sigma_vec=0;D_mat=array(0,dim=c(p1,p1));V_box=V_box[,2:ncol(V_box)]

GLL=0
  
for(j in 1:class){
  
X=XX %>% filter(X1==j) %>% select(-X1)

Y=YY %>% filter(X1==j) %>% select(-X1)

X=as.matrix(X);Y=as.matrix(Y)
  
V=X%*%D%*%t(X)+diag(sigma,n)

v_vec[,j]=c(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta))

eps=Y-X%*%eta-X%*%c(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta))

W=solve(V)-solve(V)%*%X%*%solve(eta_sub1+lambda*G)%*%t(X)%*%solve(V)

sigma_vec=sigma_vec+as.numeric(t(eps)%*%eps+sigma*(n-sigma*sum(diag(W))))

D_mat=D_mat+((D%*%t(X)%*%solve(V)%*%(Y-X%*%eta))%*%t(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta)))+D-D%*%t(X)%*%W%*%X%*%D

GLL=GLL+as.numeric(t(eps)%*%solve(diag(sigma,n))%*%eps+log(det(D)+1)+t(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta))%*%solve(D)%*%(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta))+sum(log(rep(sigma,n)))+lambda_v*t(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta))%*%solve(G)%*%(D%*%t(X)%*%solve(V)%*%(Y-X%*%eta)))
  
}

sigma=sigma_vec/(n*class)

D=D_mat/n;GLL=GLL+as.numeric(lambda*t(eta)%*%G%*%eta)

print(GLL)

}



```

```{r}


#pois B spline

library(dplyr)

data=data.frame(cell=c(1:17),x1=c(rep(1,9),rep(2,8)),x2=c(rep(1,7),2,2,1,2,2,2,2,2,2,2),x3=c(1,2,3,4,5,6,7,1,2,7,1,2,3,4,5,6,7),w=c(62.9,112.9,133.1,376.6,9.4,70.8,4.4,352.1,840.1,14.5,844.8,1296,1214.9,3740.7,109.4,404.7,66.3),y=c(270,62,68,19,0,14,228,148,82,0,111,76,30,15,37,12,15))

lambda=(10);

X=as.matrix(data[,colnames(data) %in% c("w")])/100

Z=sort(X);

B=cbind(Z,log(abs(Z)),abs(Z))

beta=rep(1,ncol(B));omega=diag(1,ncol(B))

Y=c(data$y)/1000;W=diag(1,nrow(data));Y=diag(W)*Y/diag(W)

r=rep(1,ncol(B))




ite=100000

for(j in 1:ite){

r=exp(B%*%beta)

W=W*diag(r)

y=Y/r-1+B%*%beta

beta=solve(t(B)%*%W%*%B+lambda*omega)%*%t(B)%*%W%*%y

#print(beta)

}

cor(y,log(r))

res=data.frame(y=y,log_r=log(r))

plot(c(1:nrow(res)),y,ylim=c(min(c(res$y,res$log_r)),max(c(res$y,res$log_r))),xlab="index",ylab="values",col=2)

par(new=T)

plot(c(1:nrow(res)),log(r),ylim=c(min(c(res$y,res$log_r)),max(c(res$y,res$log_r))),xlab="index",ylab="values",col=3)


```
```{r}

#B-spline pois reg

library(dplyr)

data=data.frame(cell=c(1:20),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),w=c(62.9,112.9,133.1,376.6,9.4,70.8,4.4,352.1,840.1,14.5,844.8,1296,1214.9,3740.7,109.4,404.7,66.3,60,62,70),y=c(270,62,68,19,10,14,228,148,82,10,111,76,30,15,37,12,15,38,50,60))



X=(as.matrix(data[,colnames(data) %in% c("x1","x2")]))

for(j in 1:ncol(X)){
  
X[,j]=(X[,j]-mean(X[,j]))/sd(X[,j])  
  
}

Y=y=c(data$y);w=c(data$w);w=w

B1=function(x){
  
return(sum(exp(-x)))  
  
}

B2=function(x){
  
return(sum(x^2))  
  
}

B3=function(x){
  
return(sum(x))  
  
}

B1_mat=c();B2_mat=c();B3_mat=c()

for(j in 1:nrow(X)){
  
B1_mat=c(B1_mat,B1(X[j,]))  

B2_mat=c(B2_mat,B2(X[j,]))

B3_mat=c(B3_mat,B3(X[j,]))

}

B_mat=cbind(B1_mat,B2_mat,B3_mat)

B_mat=cbind(rep(1,nrow(B_mat)),B_mat)

ite=100

beta=rep(1,ncol(B_mat))

CV_data=data.frame(lambda=c(1:20000),CV=0)

for(i in 1:nrow(CV_data)){

lambda=CV_data$lambda[i];omega=diag(1,ncol(B_mat))

for(j in 1:ite){
  
r_vec=exp(B_mat%*%beta)  

y2=y/r_vec-1+B_mat%*%beta

beta=solve(t(B_mat)%*%diag(c(w*r_vec))%*%B_mat+lambda*omega)%*%(t(B_mat)%*%diag(c(w*r_vec))%*%y2)
  

delta=2*(sum(w*y*log(y))-sum(w*y*(B_mat%*%beta))-sum(w*y)+sum(exp(B_mat%*%beta)))+lambda*sum(beta^2)

A=B_mat%*%solve(t(B_mat)%*%diag(c(w*r_vec))%*%B_mat+lambda*omega)%*%t(B_mat)%*%diag(c(w*r_vec))

#print(delta)

}  

CV=sum(w*(y-exp(B_mat%*%beta))^2/(1-diag(A)))

CV_data$CV[i]=CV

}





beta=rep(1,ncol(B_mat))

lambda=min(CV_data$lambda[CV_data$CV==min(CV_data$CV)])

for(j in 1:ite){

r_vec=exp(B_mat%*%beta)  

y2=y/r_vec-1+B_mat%*%beta

beta=solve(t(B_mat)%*%diag(c(w*r_vec))%*%B_mat+lambda*omega)%*%(t(B_mat)%*%diag(c(w*r_vec))%*%y2)


delta=2*(sum(w*y*log(y))-sum(w*y*(B_mat%*%beta))-sum(w*y)+sum(exp(B_mat%*%beta)))+lambda*sum(beta^2)

print(delta)

}
  
res=data.frame(num=c(1:nrow(X)),predict=exp(B_mat%*%beta),y=y)

plot(res$num,res$predict,type="o",col=2,xlim=c(1,nrow(X)),ylim=c(min(c(res$predict,res$y)),max(c(res$predict,res$y))),xlab="index",ylab="values")

par(new=T)

plot(res$num,res$y,type="o",col=3,xlim=c(1,nrow(X)),ylim=c(min(c(res$predict,res$y)),max(c(res$predict,res$y))),xlab="index",ylab="values")



```