# R/tau3ordered.R In CA3variants: Three-Way Correspondence Analysis Variants

```tau3ordered<-
function (f3, digits = 3)
{
####################################
#f3 three-way absolute frequency table
# ordered response variable: orderedY=T,F
#number of ordered predictor variables: orderedX=0 or 1 or 2 if orderedX=1 then consider X1 ordered and X2 not ordered
#########################
nn <- dim(f3)
ni <- nn[1]
ui <- rep(1, ni)
nj <- nn[2]
nk <- nn[3]
n <- sum(f3)
p3 <- f3/n
if (length(dim(f3)) != 3)
stop("f3 is not a 3 way table \n")
pi <- apply(p3, 1, sum)
pj <- apply(p3, 2, sum)
pk <- apply(p3, 3, sum)
pijk <- pi %o% pj %o% pk #to get the product
p1jk <- ui %o% pj %o% pk #to get the product
devt <- 1 - sum(pi^2) #tau3 denominator
tau3 <- sum(((p3 - pijk)^2/p1jk)) #tau3 numerator
itau3 <- tau3/devt #tau3 index
pij <- apply(p3, c(1, 2), sum)
pik <- apply(p3, c(1, 3), sum)
pjk <- apply(p3, c(2, 3), sum)
p1j <- ui %o% pj
p1k <- ui %o% pk
p2ij <- pi %o% pj
p2ik <- pi %o% pk
p2jk<-pj %o% pk
tauij <- sum(((pij - p2ij)^2/p1j))
itauij <- tauij/devt
tauik <- sum(((pik - p2ik)^2/p1k))
itauik <- tauik/devt
khjk <- 1/ni * (sum((pjk - (pj %o% pk))^2/(pj %o% pk)))
ikhjk <- khjk/devt
khin3 <- tau3 - tauij - tauik - khjk
ikhin3 <- khin3/devt
cost<-(n - 1) * (ni - 1) *  (1/devt)
#cat("Numerator Values of partial and total indices\n")
nom <- c("Term-IJ", "Term-IK", "Term-JK", "Term-IJK", "Term-total")
# nomI <- c("TauIJ", "TauIK", "TauJK", "TauIJK", "TauM")
x <- c(tauij, tauik, khjk, khin3, tau3)
y <- (100 * x)/tau3
dres <- (ni - 1) * (nj - 1) * (nk - 1)
dij <- (ni - 1) * (nj - 1)
dik <- (ni - 1) * (nk - 1)
djk <- (nj - 1) * (nk - 1)
dtot <- dij + dik + djk + dres
CM <- (n - 1) * (ni - 1) * itau3
Cij <- (n - 1) * (ni - 1) * itauij
Cik <- (n - 1) * (ni - 1) * itauik
Cjk <- (n - 1) * (ni - 1) * ikhjk
Cijk <- (n - 1) * (ni - 1) * ikhin3
zz <- c(itauij, itauik, ikhjk, ikhin3, itau3)
zz2 <- c(Cij, Cik, Cjk, Cijk, CM)
zz3 <- c(dij, dik, djk, dres, dtot)
pvalue<-1-pchisq(zz2,zz3)
x2<-zz2/zz3
z <- rbind(x, zz, y, zz2, zz3,pvalue,x2)
nomr <- c("Tau-Numerator", "Tau-index",
"% of Inertia", "CM-Statistic", "df","p-value","CM-statistic/df")
dimnames(z) <- list(nomr, nom)
z<-round(z, digits = digits)
# cat("C-statistic of Marcotorchino\n")
# print(CM)
# cat("Tau3 denominator\n")
# print(devt)
############################################ computing polynomials
mi<-c(1:ni)
Apoly <- emerson.poly(mi, pi)	# Emerson orthogonal polynomials
Apoly <- Apoly[, - c(1)  ] #delete zero-vec but with trivial solution to reconstruct tau
poli <- diag(sqrt(pi)) %*% Apoly	#poly with weight Di
mj <- c(1:nj)
Bpoly <- emerson.poly(mj, pj)
Bpoly <- Bpoly[,-c(1) ]
polj <- diag(sqrt(pj)) %*%(Bpoly) #poly with weight Dj
mk <- c(1:nk)
Cpoly <- emerson.poly(mk, pk)
Cpoly <- Cpoly[,-c(1) ]
polk <- diag(sqrt(pk)) %*%(Cpoly) #poly with weight Dk
#polk <- t(Cpoly)
#xs<-rstand3(f3) #standardized residuals
#xs<-p3
##########################################################partial term IJ
fij <- (pij - p2ij)/sqrt(p1j)
z3par<-   t(poli)%*%fij%*%polj
tauij<-(sum(z3par^2)/(devt))*(ni-1)*(n-1)
pvalijtot<-1 - pchisq(tauij, (ni-1)*(nj-1))
#  cat("index numerator of tau_ij reconstructed by 2 polynomials \n")
#  print(tauij)
z3ij<-(z3par^2/(devt))*(ni-1)*(n-1)
#cat("partition of index numerator of tau_ij reconstructed by 2 polynomials \n")
tauijcol<-apply(z3ij,2,sum)
tauijrow<-apply(z3ij,1,sum)
#browser()
pval<-c()
for (i in 1:(ni)){
pval[i]<-1 - pchisq(tauijrow[i], nj-1)
}
pvaltauijrow<-pval
pval<-c()
for (i in 1:(nj)){
pval[i]<-1 - pchisq(tauijcol[i], ni-1)
}
pvaltauijcol<-pval
#row and column poly component
dfi<-rep((nj-1),(ni))
perci<-tauijrow/tauij*100
zi<-cbind(tauijrow,perci,dfi,pvaltauijrow)
zitot<-c(apply(zi[,1:3],2,sum),pvalijtot)
zi<-rbind(zi,zitot)
dfj<-rep((ni-1),(nj))
percj<-tauijcol/tauij*100
zj<-cbind(tauijcol,percj,dfj,pvaltauijcol)
zjtot<-c(apply(zj[,1:3],2,sum),pvalijtot)
zj<-rbind(zj,zjtot)
zij<-rbind(zi,zj)
nomi<-paste("poly-row",0:(ni-1),sep="")
nomj<-paste("poly-col",0:(nj-1),sep="")
dimnames(zij)<-list(c(nomi,"Tau-IJ",nomj,"Tau-IJ"),c("Term-IJ-poly","%inertia","df","p-value"))
zij<-round(zij,digits=digits)
#################################################partial term IK
fik <- (pik - p2ik)/sqrt(p1k)
# z3par<-   t(poli[,-1])%*%fik%*%polk[,-1]
z3par<-   t(poli)%*%fik%*%polk
tauik<-(sum(z3par^2)/devt)*(ni-1)*(n-1)
pvaliktot<-1 - pchisq(tauik, (ni-1)*(nk-1))
z3ik<-  (z3par^2/devt)*(ni-1)*(n-1)
tauikcol<-apply(z3ik,2,sum)
tauikrow<-apply(z3ik,1,sum)
#---
pval<-c()
for (i in 1:(ni)){
pval[i]<-1 - pchisq(tauikrow[i], nk-1)
}
pvaltauikrow<-pval
pval<-c()
for (i in 1:(nk)){
pval[i]<-1 - pchisq(tauikcol[i], ni-1)
}
pvaltauikcol<-pval
#row and tube poly component
dfi<-rep((nk-1),(ni))
perci<-tauikrow/tauik*100
zi<-cbind(tauikrow,perci,dfi,pvaltauikrow)
zitot<-c(apply(zi[,1:3],2,sum),pvaliktot)
zi<-rbind(zi,zitot)
dfk<-rep((ni-1),(nk))
perck<-tauikcol/tauik*100
zk<-cbind(tauikcol,perck,dfk,pvaltauikcol)
zktot<-c(apply(zk[,1:3],2,sum),pvaliktot)
zk<-rbind(zk,zktot)
zik<-rbind(zi,zk)
nomi<-paste("poly-row",0:(ni-1),sep="")
nomk<-paste("poly-tube",0:(nk-1),sep="")
dimnames(zik)<-list(c(nomi,"Tau-IK",nomk,"Tau-IK"),c("Term-IK-poly","%inertia","df","p-value"))
zik<-round(zik,digits=digits)
################################################partial term JK
fjk <- (pjk - p2jk)/sqrt(p2jk)
z3par<-   t(polj)%*%fjk%*%polk
z3jk<- 1/ni * (z3par^2/devt)*(ni-1)*(n-1)
taujk=(sum(z3par^2)/devt)*(ni-1)*(n-1)
pvaljktot<-1 - pchisq(taujk, (nj-1)*(nk-1))
taujkcol<-apply(z3jk,2,sum)
taujkrow<-apply(z3jk,1,sum)
#---
pval<-c()
for (i in 1:(nj)){
pval[i]<-1 - pchisq(taujkrow[i], nk-1)
}
pvaltaujkrow<-pval
pval<-c()
for (i in 1:(nk)){
pval[i]<-1 - pchisq(taujkcol[i], nj-1)
}
pvaltaujkcol<-pval
#row and column poly component
dfj<-rep((nk-1),(nj))
percj<-taujkrow/taujk*100
zj<-cbind(taujkrow,percj,dfj,pvaltaujkrow)
zjtot<-c(apply(zj[,1:3],2,sum),pvaljktot)
zj<-rbind(zj,zjtot)
dfk<-rep((nj-1),(nk))
perck<-taujkcol/taujk*100
zk<-cbind(taujkcol,perck,dfk,pvaltaujkcol)
zktot<-c(apply(zk[,1:3],2,sum),pvaljktot)
zk<-rbind(zk,zktot)
zjk<-rbind(zj,zk)
nomj<-paste("poly-col",0:(nj-1),sep="")
nomk<-paste("poly-tube",0:(nk-1),sep="")
dimnames(zjk)<-list(c(nomj,"Tau-JK",nomk,"Tau-JK"),c("Term-JK-poly","%inertia","df","p-value"))
zjk<-round(zjk,digits=digits)
#################################three ordered variables and interaction term
fijk=(p3 - pijk)/sqrt(p1jk)
#browser()
p2ijk<-p2ij%o%pk
p2ikj<-p1k%o%pj
p2ikj<-aperm(p2ikj,c(1,3,2))
p2jki<-p2jk%o%ui
p2jki<-aperm(p2jki,c(3,1,2))
#z3n<- sqrt((ni-1)*(n-1))* t(poli[,-1])%*%flatten(fijk)%*%Kron(polj[,-1],polk[,-1])/sqrt(devt)
z3n<-sqrt((ni-1)*(n-1))*  ( t(poli)%*%flatten(fijk)%*%Kron(polj,polk))/sqrt(devt)
tauint<-(z3n^2)
dim(tauint)<-c(ni,nj,nk)
#dim(tauint)<-c(ni-1,nj-1,nk-1)
tautot<-(sum(z3n^2)/devt)*(ni-1)*(n-1)
pvalinttot<-1 - pchisq(tautot, (ni-1)*(nj-1)*(nk-1))
#browser()
tauintrow<-apply(tauint,1,sum)
dfi<-rep((nj-1)*(nk-1),(ni))
#          dfi<-rep((nj-1)*(nk-1),(ni-1))
pvaltauintrow<-1 - pchisq(tauintrow, ni-1)
#----------------------------------
tauintcol<-apply(tauint,2,sum)
pvaltauintcol<-1 - pchisq(tauintcol, nj-1)
#-----------------------------------
tauinttub<-apply(tauint,3,sum)
pvaltauinttub<-1 - pchisq(tauinttub, nk-1)
#-----------------------------------------
#row column and tube poly component on three-way intereraction term
perci<-tauintrow/tautot*100
zi<-cbind(tauintrow,perci,dfi,pvaltauintrow)
zitot<-c(apply(zi[,1:3],2,sum),pvalinttot)
zi<-rbind(zi,zitot)
#----------
dfj<-rep((ni-1)*(nk-1),(nj))
#    dfj<-rep((ni-1)*(nk-1),(nj-1))
percj<-tauintcol/tautot*100
zj<-cbind(tauintcol,percj,dfj,pvaltauintcol)
zjtot<-c(apply(zj[,1:3],2,sum),pvalinttot)
zj<-rbind(zj,zjtot)
#--------
#    dfk<-rep((ni-1)*(nj-1),(nk-1))
dfk<-rep((ni-1)*(nj-1),(nk))
perck<-tauinttub/tautot*100
zk<-cbind(tauinttub,perck,dfk,pvaltauinttub)
zktot<-c(apply(zk[,1:3],2,sum),pvalinttot)
zk<-rbind(zk,zktot)
zijk<-rbind(zi,zj,zk)
#    nomi<-paste("poly",1:(ni-1),sep="")
#    nomj<-paste("poly",1:(nj-1),sep="")
#    nomk<-paste("poly",1:(nk-1),sep="")
nomi<-paste("poly-row",0:(ni-1),sep="")
nomj<-paste("poly-col",0:(nj-1),sep="")
nomk<-paste("poly-tube",0:(nk-1),sep="")
#browser()
dimnames(zijk)<-list(c(nomi, "Tau-IJK",nomj,"Tau-IJK",nomk,"Tau-IJK"),c("Term-IJK-poly","%inertia","df","p-value"))
zijk<-round(zijk,digits=digits)
#============================================================
return(list(z=z,zij=zij,zik=zik,zjk=zjk,zijk=zijk,pij=pij,pik=pik,pjk=pjk,cost=cost))
}
```

## Try the CA3variants package in your browser

Any scripts or data that you put into this service are public.

CA3variants documentation built on Oct. 10, 2022, 5:07 p.m.