一変数傾向分析

```{r}

#カテゴリーごとの尤度比検定法

library(dplyr)

N=100

data1=data.frame(x1=sample(c(1:3),N,replace=TRUE,prob=c(0.3,0.1,0.6)),x2=sample(c(1:3),N,replace=TRUE,prob=c(0.2,0.2,0.6)),x3=sample(c(1:3),N,replace=TRUE,prob=c(0.3,0.3,0.4)))


data2=data.frame(x1=sample(c(1:3),N,replace=TRUE,prob=c(0.1,0.1,0.8)),x2=sample(c(1:3),N,replace=TRUE,prob=c(0.1,0.3,0.6)),x3=sample(c(1:3),N,replace=TRUE,prob=c(0.5,0.3,0.2)))


summary1_x1=data1 %>% group_by(x1) %>% summarise(n=n())

summary1_x2=data1 %>% group_by(x2) %>% summarise(n=n())

summary1_x3=data1 %>% group_by(x3) %>% summarise(n=n())

summary2_x1=data2 %>% group_by(x1) %>% summarise(n=n())

summary2_x2=data2 %>% group_by(x2) %>% summarise(n=n())

summary2_x3=data2 %>% group_by(x3) %>% summarise(n=n())

n11=summary1_x1$n;N11=sum(summary1_x1$n)

p11=n11/N11

n12=summary1_x2$n;N12=sum(summary1_x2$n)

p12=n12/N12

n13=summary1_x3$n;N13=sum(summary1_x3$n)

p13=n13/N13

n21=summary2_x1$n;N21=sum(summary2_x1$n)

p21=n21/N21

n22=summary2_x2$n;N22=sum(summary2_x2$n)

p22=n22/N22

n23=summary2_x3$n;N23=sum(summary2_x3$n)

p23=n23/N23

#x1について

L11=((gamma(N11+1)/prod(gamma(n11+1)))*prod(p11^n11))

L21=((gamma(N21+1)/prod(gamma(n21+1)))*prod(p21^n21))

if(L11>L21){

LR=-2*log(L21/L11)

}else{
  
LR=-2*log(L11/L21)  
  
}

alpha=0.05

LR>qchisq(alpha,length(unique(data1$x1,data2$x1))-1)

#x2について

L12=((gamma(N12+1)/prod(gamma(n12+1)))*prod(p12^n12))

L22=((gamma(N22+1)/prod(gamma(n22+1)))*prod(p22^n22))

if(L12>L22){

LR=-2*log(L22/L12)

}else{
  
LR=-2*log(L12/L22)  
  
}

alpha=0.05

LR>qchisq(alpha,length(unique(data1$x2,data2$x2))-1)

#x3について

L13=((gamma(N13+1)/prod(gamma(n13+1)))*prod(p13^n13))

L23=((gamma(N23+1)/prod(gamma(n23+1)))*prod(p23^n23))

if(L13>L23){

LR=-2*log(L23/L13)

}else{
  
LR=-2*log(L13/L23)  
  
}

alpha=0.05

LR>qchisq(alpha,length(unique(data1$x3,data2$x3))-1)

#一様性の検定

m1=c(sum(n11*p11),sum(n12*p12),sum(n13*p13))

m2=c(sum(n21*p21),sum(n22*p22),sum(n23*p23))

#H:data1のx1についてのカテゴリの発生率がdata2と等しいかどうか

alpha=0.05

2*sum(abs(n11*log(n11/n21)))>qchisq(alpha,N-1)

#H:data1のx2についてのカテゴリの発生率がdata2と等しいかどうか

alpha=0.05

2*sum(abs(n12*log(n12/n22)))>qchisq(alpha,N-1)

#H:data1とdata2のx1における平均頻度が等しいかどうか

alpha=0.05

2*sum(((m2+m1)/2)*abs(log(m2/m1)))>qchisq(alpha,N-1)

#u:主効果

u11_ave=mean(log(n11))

u11_vec=log(n11)-u11_ave

u12_ave=mean(log(n12))

u12_vec=log(n12)-u12_ave

u13_ave=mean(log(n13))

u13_vec=log(n13)-u13_ave

u21_ave=mean(log(n21))

u21_vec=log(n21)-u21_ave

u22_ave=mean(log(n22))

u22_vec=log(n22)-u22_ave

u23_ave=mean(log(n23))

u23_vec=log(n23)-u23_ave


X=array(0,dim=c(ncol(data1)+1,ncol(data1)+1))

diag(X)=rep(1,ncol(data1)+1);X[,1]=rep(1,ncol(data1)+1)

X%*%c(u11_ave,u11_vec)

#u11_vecの(0.05,0.95)％区間の区間推定

alpha=0.05

print(paste0(log(qpois(alpha,n11)*exp(-u11_ave))," , ",log(qpois(1-alpha,n11)*exp(-u11_ave))))

#u12_vecの(0.05,0.95)％区間の区間推定

print(paste0(log(qpois(alpha,n12)*exp(-u12_ave))," , ",log(qpois(1-alpha,n12)*exp(-u12_ave))))

#u13_vecの(0.05,0.95)％区間の区間推定

print(paste0(log(qpois(alpha,n13)*exp(-u13_ave))," , ",log(qpois(1-alpha,n13)*exp(-u13_ave))))

#u21_vecの(0.05,0.95)％区間の区間推定

alpha=0.05

print(paste0(log(qpois(alpha,n21)*exp(-u21_ave))," , ",log(qpois(1-alpha,n21)*exp(-u21_ave))))

#u12_vecの(0.05,0.95)％区間の区間推定

print(paste0(log(qpois(alpha,n22)*exp(-u22_ave))," , ",log(qpois(1-alpha,n22)*exp(-u22_ave))))

#u13_vecの(0.05,0.95)％区間の区間推定

print(paste0(log(qpois(alpha,n23)*exp(-u23_ave))," , ",log(qpois(1-alpha,n23)*exp(-u23_ave))))

```