#統計検定量から多次元尺度(torgerson method)へ
library(dplyr)
data=data.frame(処理=c(rep("M",12),rep("P",11)),sex=c(1,1,1,1,1,1,2,2,2,2,2,2,1,1,1,1,1,1,2,2,2,2,2),t1=c(317,186,377,229,276,272,219,260,284,365,298,274,232,367,253,230,190,290,337,283,325,266,338),t2=c(280,189,395,258,310,250,210,245,256,304,321,245,205,354,256,218,188,263,337,279,257,258,343),t3=c(275,190,368,282,306,250,236,264,241,294,341,262,244,358,247,245,212,291,383,277,288,253,307),t4=c(270,135,334,272,309,255,239,268,242,287,342,263,197,333,228,215,201,312,318,264,326,284,274),t5=c(274,197,338,264,300,228,242,317,243,311,357,235,218,338,237,230,169,299,361,269,293,245,262),t6=c(266,205,334,265,264,250,221,314,241,302,335,246,233,355,235,207,179,279,341,271,275,263,309))
k=ncol(data)-2;n=nrow(data)
Y=as.matrix(data[,!(colnames(data) %in% c("処理","sex"))])
tD=array(0,dim=c(k-1,k))
for(j in 1:(k-1)){
tD[j,j:(j+1)]=c(-1,1)
}
D=t(tD)
diff=solve(t(D)%*%D)%*%t(D)
d_mat=array(0,dim=c(n,n))
mat=array(0,dim=c(n,n))
kai_mat=array(0,dim=c(n,n))
for(j in 1:n){
for(i in 1:n){
if(i!=j){
d_mat[i,j]=sum(((diff)%*%c(Y[i,]-Y[j,]) )^2)
mat[i,j]=sum(((diff)%*%c(Y[i,]-Y[j,]) ))
kai=0
for(l in 1:(k-1)){
kai=kai+(k*l)*(sum(Y[i,1:l]-Y[j,1:l])/l-(mean(Y[i,])-mean(Y[j,])))^2/(k-l)
}
kai_mat[i,j]=kai/2
}
}
}
cor_mat=array(0,dim=c(n,n))
for(j in 1:n){
for(i in 1:n){
if(i!=j){
coef=(diff)%*%c(Y[i,]-Y[j,])
cor_mat[i,j]=cor(Y[i,]-Y[j,],D%*%coef)
print(paste0("columns:",i,",",j))
print((Y[i,]-Y[j,]-mean(Y[i,]-Y[j,]))/sd(Y[i,]-Y[j,]))
print((c(D%*%coef)-mean(c(D%*%coef)))/sd(c(D%*%coef)))
}
}
}