RAの新診断基準をいじってみる

• シミュレーションで臨床情報とＳＮＰ多型データを作って、それを、「木」にしてみる
• package "sphere"が必要(こちら)

RegularSpherize2<-
function (O = matrix(round(runif(2 * 3) * 100, 0), 2, 3)) 
{
    EandM <- tableExpAndMarginals(O)
    N <- length(O[, 1])
    M <- length(O[1, ])
    E <- EandM$etable
    D <- EandM$dtable
    K <- sum(D^2/E)
    X <- CategoryVector2D(N, M)
    P <-  (length(D[, 1]) - 1) * (length(D[1, ]) - 1)/(length(D[, 1]) * 
        length(D[1, ])) * c(t(X) %*% matrix(c(t(D)), length(D), 
        1))
    V <-  t(X) %*% diag(1/(c(t(E) * K))) %*% X

    Xe <- X/c(t(sqrt(E)))
    EigenOut <- eigen(V)
    U <- EigenOut$vectors
    Ui <- solve(U)
    Pi <- Ui %*% P
    Mi <- diag(sqrt(EigenOut$values))
    if (length(EigenOut$values) == 1) {
        Mi <- matrix(1, 1, 1)
    }
    Pii <- Mi %*% Pi
    list(O = O, E = E, D = D, K = K, X = X, V = V, Xe = Xe, P = P, 
        EigenOut = EigenOut, Ui = Ui, Mi = Mi, Pi = Pi, Pii = Pii)
}

#############
# PはRA診断基準の情報
# P[,1]は腫脹大関節数
# P[,2]は腫脹小関節数
# P[,3]はリウマトイド因子(陰性:0,陽性:1,強陽性:2)
# P[,4]はACPA(陰性:0,陽性:1,強陽性:2)
# P[,5]は関節炎の持続期間(単位週)
# P[,6]はESR(正常:0,異常:1)
# P[,7]はCRP(正常:0,異常:1)
#############


RAscore<-function(P){
 Score<-Joint<-Autoantibody<-Inflammation<-Duration<-rep(0,length(P[,1]))
 Joint<-calcJoint(P[,1:2])
 Autoantibody<-apply(P[,3:4],1,FUN="max")
 Autoantibody[which(Autoantibody==2)]<-3
 Autoantibody[which(Autoantibody==1)]<-2
 Inflammation<-apply(P[,6:7],1,FUN="max")
 
 Duration[which(P[,5]>=6)]<-1
 
 Score<-Joint+Autoantibody+Inflammation+Duration
 list(Score=Score,Joint=Joint,Autoantibody=Autoantibody,Inflammation=Inflammation,Duration=Duration)
}

calcJoint<-function(Pj){
 Psum<-Pj[,1]+Pj[,2]
 ret<-rep(0,length(Pj[,1]))
 ret[which(Psum>=1 & Pj[,2]==0)]<-1
 ret[which(Pj[,2]>=1 & Pj[,2]<=3 & Psum<=10)]<-2
 ret[which(Pj[,2]>=4 & Pj[,2]<=10 & Psum<=10)]<-3
 ret[which(Psum>10 & Pj[,2]>=1)]<-5
 ret
}


Np<-7
Ng<-1

Ns<-500

#P<-matrix(sample(c(0,1,2,3),Ns*Np,replace=TRUE),Ns,Np)
P<-matrix(0,Ns,Np)
P[,1]<-sample(c(0,1,2,3),Ns,replace=TRUE)
P[,2]<-rpois(Ns,5)
P[,3]<-sample(c(0,1,2),Ns,replace=TRUE)
P[,4]<-sample(c(0,1,2),Ns,replace=TRUE)
P[,5]<-runif(Ns)*10
P[,6]<-sample(c(0,1),Ns,replace=TRUE)
P[,7]<-sample(c(0,1),Ns,replace=TRUE)

#G<-matrix(sample(c(0,1,2),Ns*Ng,replace=TRUE),Ns,Ng)
G<-matrix(0,Ns,Ng)
RAsc<-RAscore(P)

G[which(RAsc$Score>=3)]<-1
G[which(RAsc$Score>=6)]<-2
O<-matrix(0,3,Ns)
for(i in 0:2){
 O[i+1,which(G==i)]<-1
}

Ws<-matrix(0,Np,Ns*3)
for(i in 1:length(Ws[,1])){
 for(j in 2:3){
  Ws[i,(Ns*(j-1)+1):(Ns*j)]<-(j-1)*P[,i]
 }
}
Sp<-RegularSpherize2(O)
df<-length(Sp$Pii)
Ts<-matrix(0,length(Ws[,1]),df)
for(i in 1:length(Ws[,1])){
 Ts[i,]<-TestDf1(O,Ws[i,],Sp=Sp)$Tiist
}

SpPst<-Sp$Pii/nrm(Sp$Pii)
PiTs<-cbind(SpPst,t(Ts))

corrr<-t(PiTs)%*%PiTs
corrr[which(corrr>1,arr.ind=TRUE)]<-1

acoscor<-acos(corrr)
library(ape)
njout<-nj(acoscor)

plot(njout,type="u")

#GenMatOut1<-GeneralMatrixTest4(O=O,Ws=Ws,permutation=FALSE)
#GenMatOut2<-GeneralMatrixTest4(O=O,Ws=Ws,permutation=TRUE)


##################
NL<-0:15
NS<-0:15

NLS<-expand.grid(NL,NS)

NLSA<-cbind(NLS,NLS[,1]+NLS[,2])


col<-rep("black",length(NLSA[,1]))

col[which(NLSA[,3]<=1)]<-"gray"
col[which(NLSA[,2]==0 & NLSA[,3]>=1)]<-"blue"
col[which(NLSA[,2]>=1 & NLSA[,2]<=3 & NLSA[,3]<=10)]<-"purple"
col[which(NLSA[,2]>=4 & NLSA[,2]<=10 & NLSA[,3]<=10)]<-"red"
col[which(NLSA[,3]>10& NLSA[,2]>=1)]<-"black"
plot(NLSA[,1],NLSA[,2],col=col,pch=15,xlab="NL",ylab="NS")