社会環境情報の計数データの実践的解析法　浅野長一郎・竹内光悦先生　共立出版

```{r}

#p.12

#Ridit 尺度

data=data.frame(category=c("非常に悪い","悪い","やや悪い","どちらでもない","やや良い","良い","非常に良い"),G=paste0(rep("G",7),c(1:7)),freq=c(17,54,60,19,9,6,14))

N=sum(data$freq)

freq1=data$freq/2

freq2=c(0,cumsum(data$freq)[-length(data$freq)])

freq=(freq1+freq2)/N


0.5-sqrt(1/(3*N))

0.5+sqrt(1/(3*N))

```

```{r}

#probit 変換

#p.13

x=c(1:5)

n=c(20,20,30,20,20)

r=c(0,2,16,15,19)

p=r/n

a=1;b=1

ite=10^(5)

eta=0.0001

h=0.01

for(l in 1:ite){
  
da=(pnorm((a+h)+b*x)-pnorm(a+b*x))/h

db=(pnorm(a+(b+h)*x)-pnorm(a+b*x))/h

P=pnorm(a+b*x)

a=a+eta*sum(n*(p-P)*da/(P*(1-P)))
  
b=b+eta*sum(n*(p-P)*db/(P*(1-P)))  

loglik=sum(r*log(P)+(n-r)*log(1-P))

print(loglik)

}

sig=1/b

mu=(5-a)*sig


```
```{r}

#社会環境情報の計数データの実践的解析法　p,22

library(dplyr)

X=matrix(c(3,2,2,3,2,3,2,1,2,3,2,2,1,3,3,3,1,2,1,3,2,3,2,3,2,3,2,1,2,1,3,4,2,3,2,3,2,4,3,2,4,3,2,3,3,4,3,2,1,2,1,1,3,2,3,3,3,2,2,1,2,2,2,2,3,3,3,3,2,3,4,3,2,3,4,4,3,4,3,3),ncol=8)

k=ncol(X);r=nrow(X)

rank=sort(unique(c(X)))

mat1=array(0,dim=c(nrow(X),length(rank)))

mat2=(array(0,dim=c(ncol(X),length(rank))))

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

for(i in rank){  

  mat1[j,i]=sum(X[j,]==i)  

}

}  

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

for(i in rank){  

  mat2[j,i]=sum(X[,j]==i)  

}

}  

n..=apply(mat1,2,sum)

B=diag(n..)

A=t(t(n..))%*%t(n..)

W=((k*r*B)-A)

C=t(mat1)%*%mat1

D=t(mat2)%*%mat2

V=r*C-A+k*D-A

#W%*%solve(V)






#再急降下法

f=function(x){

val=t(c(x))%*%(W-V)%*%t(t(c(x)))

return(c(val))

}

X=c(0,rep(0.5,2),1)

ite=10000

h=0.001

eta=10^(-6)

for(l in 1:ite){
  
vec=X

for(j in 2:(length(X)-1)){
  
vec_sub=vec;vec_sub[j]=vec_sub[j]+h

X[j]=X[j]-eta*(f(vec_sub)-f(vec))/h
  
}

print(f(X))

}


```

```{r}


#第3章

#順位による方法　p.25

#符号検定

#ビールAとビールBを好んだ人（A,B)とどちらでもないCがいたとき、

A=31;B=16;C=3;n=A+B+C

k=min(A,B)

#統計量

static=2*(ifelse(abs(k+0.5-n/2)>abs(k-0.5-n/2),k-0.5-n/2,k+0.5-n/2))/sqrt(n)

#嗜好に差がある

abs(static)>qnorm(1-0.05/2)



#3.2　線上の連の数による検定

#2つの母集団が同じものとみなせるのかどうか（同等性）

A=c(41,52,45,47,43)

B=c(62,50,49,59,54,46)

seq=sort(c(A,B))

nA=length(A);nB=length(B)

r=min(length(A),length(B))

u=c(2:r)

p=c()

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

if(u[j]/2==floor(u[j]/2)){

k=u[j]/2

p=c(p,2*choose(nA-1,k-1)*choose(nB-1,k-1)/choose(nA+nB,nA))

}else{

k=(u[j]+1)/2

p=c(p,(choose(nA-1,k-1)*choose(nB-1,k-2)+choose(nB-1,k-1)*choose(nA-1,k-2))/choose(nA+nB,nA))

}

}

sum(p)

#median test p.31

mat=matrix(c(26,28,24,22),ncol=2)

N=sum(mat)

t=N*(abs(det(mat))-N/2)^2/prod(c(apply(mat,1,sum),apply(mat,2,sum)))

qchisq(1-0.05,df=1)

#p.32 拡張median test

a1=c(2,1,0,2,0,3,2,2,1,0,0);n1=length(a1)

a2=c(2,1,3,0,3,0,3,1,4,2,2,1,2,3,4);n2=length(a2)

a3=c(1,2,2,2,1,2,2,1,3,0,2,1,1);n3=length(a3)

a4=c(0,2,1,2,0,1,0,0,1,2,0);n4=length(a4)

samples=c(a1,a2,a3,a4);ns=c(n1,n2,n3,n4)

m=c(length(a1[a1<median(samples)]),length(a2[a2<median(samples)]),length(a3[a3<median(samples)]),length(a4[a4<median(samples)]))

N=length(samples)

a=ifelse(floor(N/2)==N/2,N/2,(N-1)/2)

kai2=N*(N-1)*sum((1/ns)*(m-a*ns/N)^2)/(a*(a-1))

qchisq(1-0.05,3)

```

```{r}

#Thurstone-Mosteller の方法

user=8

n=matrix(c(NA,6,5,7,6,2,NA,3,5,6,3,5,NA,4,5,1,3,4,NA,4,2,2,3,4,NA),ncol=5)

p=n/(user)

d=qnorm(p)

d2=d;diag(d2)=0

mu=c()

for(j in 1:nrow(d2)){
  
mu=c(mu,mean(d2[j,]))  
  
}

#適合度の評価

d_hat=array(0,dim=c(length(mu),length(mu)))

for(j in 1:length(mu)){
  
d_hat[j,]=mu[j]-mu  
  
}

p_hat=pnorm(d_hat)

diag(p_hat)=0

sita=asin(sqrt(p))

sita_hat=asin(sqrt(p_hat))

sita_diff=(sita-sita_hat)^2;diag(sita_diff)=0

kai2=sum(sita_diff)*user/821

#自由度
pthi=(nrow(n)-1)*(nrow(n)-2)/2

qchisq(1-0.05,df=pthi)

```

```{r}

#scheffeの一対比較法

mat=matrix(c(0,0,0,1,0,1,0,3,0,2,0,2,1,1,2,1,1,0,2,0,1,0,1,1,1,0,1,0,2,0),ncol=5,nrow=6)

points=c(-2,-1,0,1,2)


total_num=4

mu=mat%*%points/total_num

nums=data.frame(t1=c(1,2,1,3,2,3),t2=c(2,1,3,1,3,2),mu)

k=3

pis=array(0,dim=c(k,k))

delta=array(0,dim=c(k,k))

alpha=c()

for(i in 1:k){
for(j in 1:k){  
if(i<j){  
  
pis[i,j]=(nums$mu[nums$t1==i & nums$t2==j]-nums$mu[nums$t1==j & nums$t2==i])/2  
  
delta[i,j]=(nums$mu[nums$t1==i & nums$t2==j]+nums$mu[nums$t1==j & nums$t2==i])/2

}}}

pis=pis-t(pis)

alpha=apply(pis,1,sum)/k

r=array(0,dim=c(k,k))

for(i in 1:k){
  
r[i,]=pis[i,]-(alpha[i]-alpha)  
  
}

#全平方和

S0=sum(apply(mat,2,sum)*points^2)

#muの項

S_mu=sum(total_num*mu^2)

#誤差項

Se=S0-S_mu

#pisの項

S_pis=2*total_num*sum(pis[upper.tri(pis)]^2)

#主効果

S_alpha=2*total_num*k*sum(alpha^2)

#組み合わせ効果

Sr=S_pis-S_alpha

#順序効果

S_delta=S_mu-S_pis



```