ロジスティック回帰






```{r}

#bernoulli-logistic regression p.105

#one dimension

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

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

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x2")]))

pis=data$y

Y=as.numeric(data$sign)

alpha=0.01;beta=0.01

ite=10

alpha_beta_data=data.frame(ite=1:ite,alpha=0,beta=0)

for(j in 1:ite){

V=array(0,dim=c(nrow(data),nrow(data)))
  
diag(V)=((1/(1+exp(-(X%*%c(alpha,beta)))))*(1-1/(1+exp(-(X%*%c(alpha,beta))))))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(pis-(1/(1+exp(-(X%*%c(alpha,beta))))))

alpha=alpha+mat[1];beta=beta+mat[2]

print(alpha);print(beta)

alpha_beta_data$alpha[j]=alpha

alpha_beta_data$beta[j]=beta

}

data=data %>% dplyr::mutate(predict=1/(1+exp(-(X%*%c(alpha,beta)))))

library(mvtnorm)

normal_dis_data=data.frame(x=0,y=0,dis=0)

for(j in seq(10,12,0.01)){
  
normal_dis_data=rbind(normal_dis_data,data.frame(x=j,y=seq(-2,2,0.01),dis=0))  

}

for(j in 1:nrow(normal_dis_data)){
  
normal_dis_data$dis[j]=dmvnorm(c(normal_dis_data$x[j],normal_dis_data$y[j]),c(alpha,beta),(t(X)%*%V%*%(X)))  
  
}

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

library(rgl)

plot3d(normal_dis_data$x,normal_dis_data$y,normal_dis_data$dis,col=1,type="h",xlab="x",ylab="y",zlab="density_func")

#H0:beta=0

log_lambda=-2*(sum(Y*log(mean(Y)))+sum((1-Y)*log(1-mean(Y)))-sum(Y*(X%*%c(alpha,beta)))+sum(log(1+exp(X%*%c(alpha,beta)))))

alpha=0.05

if(log_lambda>qchisq(1-alpha,df=1)){
  
print(paste0("Hypothesis is rejected"))  
  
}

data=data.frame(data,predict_sign=ifelse(data$predict>0.49,1,0))

#likelihood ratio

L_real=prod((pis^Y)*((1-pis)^(1-Y)))

L_predict=prod(((data$predict)^data$predict_sign)*(1-data$predict)^(1-data$predict_sign))

likelihood_ratio=L_predict/L_real

```

```{r}

#bernoulli-logistic regression p.105

#multi dimension

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

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

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

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

pis=data$y

Y=as.numeric(data$sign)

alpha=rep(0.01,(ncol(X)))

ite=10

alpha_beta_data=array(0,dim=c(nrow(X),ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(data),nrow(data)))
  
diag(V)=((1/(1+exp(-(X%*%alpha))))*(1-1/(1+exp(-(X%*%alpha)))))

print(sum(V))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(pis-(1/(1+exp(-(X%*%alpha)))))

alpha=alpha+mat

alpha_beta_data[j,]=alpha

}

data=data %>% dplyr::mutate(predict=1/(1+exp(-(X%*%alpha))))

#H0:alpha=0_vector(Wald test)

a=0.05

t(alpha)%*%solve(t(X)%*%V%*%X)%*%alpha < qchisq(1-a,length(alpha)-1)

#H0:alpha=0_vector(LR test)

LR=function(s){

log_estimate=log(prod(((1/(1+exp(-(X%*%alpha))))^data$sign)*(1-1/(1+exp(-(X%*%alpha))))^abs(data$sign)))

value=log(prod(((1/(1+exp(-(X%*%s))))^data$sign)*(1-1/(1+exp(-(X%*%s))))^abs(data$sign)))

return(-2*(value-log_estimate))

}

LR(c(4,0.2,-0.2))<qchisq(1-a,length(alpha)-1)

```



```{r}

library(dplyr)

#bernoulli-logistic regression 

#multi dimension

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

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

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

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

pis=as.numeric(data$sign)

n=sample(c(10:20),nrow(data),replace=T)

d=c()

for(j in 1:length(n)){

d=c(d,sample(c(1:n[j]),1))

}

alpha=rep(0.01,(ncol(X)))

ite=10000

eta=0.01

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

p=(1/(1+exp(-(X%*%alpha))))

V=diag(c(n*p*(1-p)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(d-n*p)

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

log_lik=sum(d*log(p))+sum((n-d)*log(1-p))

print(log_lik)

}

data=data %>% dplyr::mutate(predict=1/(1+exp(-(X%*%alpha))))

#Pearson kai2 static

#H:ロジスティク回帰モデルがデータにフィットしている

X2=sum((d-n*p)^2/(n*p*(1-p)))

qchisq(1-0.01,df=nrow(X)-ncol(X)-1)

```


```{r}

#bernoulli-logistic(通常)の場合

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

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

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

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

pis=data$y

Y=as.numeric(data$sign)

alpha=rep(0.01,(ncol(X)))

ite=10

alpha_beta_data=array(0,dim=c(nrow(X),ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(data),nrow(data)))

diag(V)=((1/(1+exp(-(X%*%alpha))))*(1-1/(1+exp(-(X%*%alpha)))))

print(sum(V))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(pis-(1/(1+exp(-(X%*%alpha)))))

alpha=alpha+mat

alpha_beta_data[j,]=alpha

}

data=data %>% dplyr::mutate(predict=1/(1+exp(-(X%*%alpha))))

#H0:alpha=0_vector(Wald test)

a=0.05

t(alpha)%*%solve(t(X)%*%V%*%X)%*%alpha < qchisq(1-a,length(alpha)-1)

#H0:alpha=0_vector(LR test)

LR=function(s){

log_estimate=log(prod(((1/(1+exp(-(X%*%alpha))))^data$sign)*(1-1/(1+exp(-(X%*%alpha))))^abs(data$sign)))

value=log(prod(((1/(1+exp(-(X%*%s))))^data$sign)*(1-1/(1+exp(-(X%*%s))))^abs(data$sign)))

return(-2*(value-log_estimate))

}

LR(c(4,0.2,-0.2))<qchisq(1-a,length(alpha)-1)
```

```{r}

#plot(実際の標本:黒、最小二乗法:赤、ベルヌーイロジスティック:緑)

data=data.frame(num=c(1:nrow(linear_logit)),y=linear_logit$y,linear_predict=linear_logit$predict,bernoulli_predict=bernoulli_logit$predict)

plot(data$num,data$y,type="o",col=1,xlim=c(1,nrow(data)),ylim=c(min(data$y,data$linear_predict,data$bernoulli_predict),max(data$y,data$linear_predict,data$bernoulli_predict)),xlab="index",ylab="prob")

par(new=T)

plot(data$num,data$linear_predict,type="o",col=2,xlim=c(1,nrow(data)),ylim=c(min(data$y,data$linear_predict,data$bernoulli_predict),max(data$y,data$linear_predict,data$bernoulli_predict)),xlab="index",ylab="prob")

par(new=T)

plot(data$num,data$bernoulli_predict,type="o",col=3,xlim=c(1,nrow(data)),ylim=c(min(data$y,data$linear_predict,data$bernoulli_predict),max(data$y,data$linear_predict,data$bernoulli_predict)),xlab="index",ylab="prob")

```


```{r}

#SMR model p.119

I=5;J=9

d_mat=array(0,dim=c(I,J))

n_mat=array(0,dim=c(I,J))

a=rep(1,I);b=rep(1,J)

for(i in 1:I){

val=rpois(J,sample(10:100,1))  
  
n_mat[i,]=n_mat[i,]+val

}

for(j in 1:J){
  
val=rpois(I,sample(10:100,1))  
  
n_mat[,j]=n_mat[,j]+val  
  
}

for(i in 1:I){
for(j in 1:J){  
  
d_mat[i,j]=n_mat[i,j]-sample(1:(n_mat[i,j]-1),1)
    
}
}

ite=100

mu=1

for(l in 1:ite){
  
a_pre=a;b_pre=b  
  
di_plus=apply(d_mat,1,sum)

a=log(di_plus/n_mat%*%exp(b))

dj_plus=apply(d_mat,2,sum)

b=log(dj_plus/t(n_mat)%*%exp(a))

print(sum(abs(a_pre-a))+sum(abs(b_pre-b)))

pij=d_mat/n_mat;log_pij=log(pij)

predict=array(0,dim=c(I,J))

for(i in 1:I){
for(j in 1:J){  
 
predict[i,j]=a[i]+b[j]   
  
}
}

mu=mean(log_pij-predict)

predict=predict+mu

print(sum(abs(pij-exp(predict))))

}




```

```{r}

#binomial logistic p.183

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

n_vec=rpois(nrow(data),sample(10:30,1))

d_vec=c()

for(j in 1:length(n_vec)){
  
d_vec=c(d_vec,sample(1:(n_vec[j]-1),1))  
  
}

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

pis=data$y

alpha=rep(0.01,(ncol(X)))

ite=10

alpha_beta_data=array(0,dim=c(nrow(X),ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(data),nrow(data)))

diag(V)=(n_vec*(1/(1+exp(-(X%*%alpha))))*(1-1/(1+exp(-(X%*%alpha)))))

print(sum(V))



mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(d_vec-n_vec*(1/(1+exp(-(X%*%alpha)))))

alpha=alpha+mat

alpha_beta_data[j,]=alpha

}

data=data %>% dplyr::mutate(d=d_vec,predict=n_vec*(1/(1+exp(-(X%*%alpha)))))


#solution2 p.193

V=array(0,dim=c(nrow(data),nrow(data)))

diag(V)=(n_vec*(1/(1+exp(-(X%*%alpha))))*(1-1/(1+exp(-(X%*%alpha)))))

q_vec=n_vec*(1/(1+exp(-(X%*%alpha))))

z=X%*%alpha+solve(V)%*%(d_vec-q_vec)

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

X2=sum((d_vec-q_vec)^2/diag(V))

```

```{r}

#層別解析

library(dummies)

data=data.frame(num=1:20,y=sample(0:1,20,replace=T,prob=c(0.3,0.7)),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),e=c(rep(1,10),rep(2,5),rep(3,5)))

E=dummy(data$e)

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

X=t((t(X)-apply(X,2,mean))/apply(X,2,sd))

Y=data$y

beta=rep(1,ncol(X));r=rep(1,ncol(E))

ite=300000;eta=0.00001

loglik_vec=c()

for(l in 1:ite){
  
for(i in 1:length(beta)){
  
beta_sub=beta;

beta_sub[i]=beta_sub[i]+0.01  
  
loglik1=as.numeric(beta_sub%*%t(X)%*%Y+r%*%t(E)%*%Y)-sum(log(1+exp(beta_sub%*%t(X)+r%*%t(E))))

loglik2=as.numeric(beta%*%t(X)%*%Y+r%*%t(E)%*%Y)-sum(log(1+exp(beta%*%t(X)+r%*%t(E))))

beta[i]=beta[i]+eta*(loglik1-loglik2)/0.01  
  
}  
  
for(i in 1:length(r)){
  
r_sub=r;

r_sub[i]=r_sub[i]+0.01  
  
loglik1=as.numeric(beta%*%t(X)%*%Y+r_sub%*%t(E)%*%Y)-sum(1+exp(beta%*%t(X)+r_sub%*%t(E)))

loglik2=as.numeric(beta%*%t(X)%*%Y+r%*%t(E)%*%Y)-sum(1+exp(beta%*%t(X)+r%*%t(E)))

r[i]=r[i]+eta*(loglik1-loglik2)/0.01  
  
}  
  
log_lik=as.numeric(beta%*%t(X)%*%Y+r%*%t(E)%*%Y)-sum(1+exp(beta%*%t(X)+r%*%t(E)))  
  
#print(log_lik)

loglik_vec=c(loglik_vec,log_lik)

}

```

```{r}

#logisticにおけるクラス分類

library(dplyr);library(dummies)

data=data.frame(iris)

data$Species=as.integer(data$Species)

Y=dummy(data$Species)

X=data[,!(colnames(data) %in% c("Species"))]

alpha_box=array(0,dim=c(length(unique(data$Species)),ncol(X)+1))

log_lik=c()

X=as.matrix(cbind(rep(1,nrow(X)),X))

p_value=array(0,dim=c(nrow(data),ncol(Y)))

ite=20

for(i in 1:length(unique(data$Species))){

alpha=rep(0.01,ncol(X));pis=Y[,i]

for(j in 1:ite){

p=1/(1+exp(-X%*%alpha))

V=array(0,dim=c(nrow(data),nrow(data)))

diag(V)=p*(1-p)

print(sum(V))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(pis-p)

alpha=alpha+mat


}

p_vec=p;inv_p_vec=1-p_vec

p_vec[p_vec==0]=(10^(-20));inv_p_vec[inv_p_vec==0]=(10^(-20))

log_lik=c(log_lik,as.numeric(pis*log(p_vec)+(1-pis)*log(inv_p_vec)))

alpha_box[i,]=alpha

p_value[,i]=p

}

colnames(p_value)=c("p1","p2","p3")

data=data.frame(data,p_value)

```
```{r}

#多変量解析入門　小西

#多群非線形ロジスティク判別

library(dummies)

data(iris)

X=as.matrix(iris[,-ncol(iris)]);Y=dummy(iris[,ncol(iris)])

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

w=rep(1,(ncol(Y)-1)*ncol(X))



ite=10000

eta=0.001

for(l in 1:ite){
  
w_mat=matrix(w,ncol=ncol(Y)-1)  
  
pis2=cbind(exp(X%*%w_mat)/(1+apply(exp(X%*%w_mat),1,sum)),1/(1+apply(exp(X%*%w_mat),1,sum)))  
  
pis=cbind(exp(X%*%w_mat)/(1+apply(exp(X%*%w_mat),1,sum)))

mat=Y[,-ncol(Y)]-pis

dw=w_mat

for(j in 1:(ncol(Y)-1)){
  
X2=X*c(mat[,j]);X2=apply(X2,2,sum)  
  
dw[,j]=c(X2)  

}

dw2=array(0,dim=c(ncol(X)*(ncol(Y)-1),ncol(X)*(ncol(Y)-1)))

for(i in 1:(ncol(Y)-1)){
for(j in 1:(ncol(Y)-1)){ 
  
if(i==j){
  
dw2[((i-1)*ncol(X)+1):(i*ncol(X)),((i-1)*ncol(X)+1):(i*ncol(X))]=-t(X*c(pis[,i]))%*%(X*c(1-pis[,j]))
  
}else{
  
dw2[((i-1)*ncol(X)+1):(i*ncol(X)),((j-1)*ncol(X)+1):(j*ncol(X))]=t(X*c(pis[,i]))%*%(X*c(pis[,j]))  
  
}

}}

w=w-eta*solve(dw2)%*%c(dw)

#print(sum(abs(Y-pis2)))
  
}

```
```{r}

#情報量基準　p.93

#logistic

#multi dimension

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

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

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]/10000))

pis=data$y

Y=as.numeric(data$sign)


GIC_data=data.frame(lam=seq(0.01,100,0.01),GIC=0)

alpha_mat=array(0,dim=c(nrow(GIC_data),ncol(X)))

#for(s in 1:nrow(GIC_data)){

alpha=rep(0.01,(ncol(X)))

ite=1000;s=100

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

lam=GIC_data$lam[s];eta=0.01

K=diag(1,length(alpha))

GIC_val_pre=10^8;GIC_val=10^8-1

while(abs(GIC_val-GIC_val_pre)>10^(-6))({

GIC_val_pre=GIC_val  

p=1/(1+exp(-X%*%alpha))

da=t(X)%*%diag(c(pis-p))%*%t(t(rep(1,nrow(X))))-nrow(X)*lam*K%*%alpha

da2=-t(X)%*%diag(c(p))%*%(diag(1,nrow(X))-diag(c(p)))%*%X-nrow(X)*lam*K

alpha=alpha-eta*solve(da2)%*%da

alpha_beta_data[j,]=alpha

log_lik=sum(pis*log(p)+(1-pis)*log(1-p))-nrow(X)*lam*as.numeric(t(alpha)%*%K%*%alpha)/2

R=t(X)%*%diag(c(p))%*%(diag(1,nrow(X))-diag(c(p)))%*%X/nrow(X)+lam*K

Q=t(X)%*%diag(c(pis-p)^2)%*%X-lam*K%*%alpha%*%t(rep(1,nrow(X)))%*%diag(c(p))%*%X

Q=Q/nrow(X)

GIC=-2*sum(pis*log(p)+(1-pis)*log(1-p))+2*sum(diag(solve(R)%*%Q))

GIC_val=GIC

})

GIC_data$GIC[s]=GIC

alpha_mat[s,]=c(alpha)

print(paste0(s,":",lam,",",GIC))

#}

```