計測データ（薬剤）

#統計検定量から多次元尺度(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)))



}

}
}  

#torgerson method

library(dplyr)

mat=kai_mat

d=mat

ave_ncol_d=apply(d^2,1,mean)

ave_nrow_d=apply(d^2,2,mean)

B=array(0,dim=c(nrow(mat),ncol(mat)))

for(i in 1:nrow(B)){
for(j in 1:ncol(B)){  

B[i,j]=(-1/2)*((d[i,j]^2)-ave_ncol_d[i]-ave_nrow_d[j]+mean(d^2))  

}
}

eigen_values_B=eigen(B)$values

eigen_vectors_B=eigen(B)$vectors

eigen_values_B=eigen_values_B[1:2]

eigen_vectors_B=eigen_vectors_B[,1:2]

X=eigen_vectors_B%*%diag(sqrt(eigen_values_B),length(eigen_values_B))

plot(X[,1],X[,2],col=2,type="p",xlab="principal component 1",ylab="principal component 2", main="torgerson method")