check.sphere <- function(y,groups=NULL,trial=NULL,id=NULL) # several tests on sphericity: # - Mauchly # (see e.g. Winer (1991): Statistical principles in experimental Design, p 255) # - multisample Mauchly # (Mendoza, Jorge L. (1980): A significance test for multisample sphericity, Psychometrika, Vol 45, No 4 ) # - Circularity (Compound Symmetry) test # (see e.g. Winer (1991): Statistical principles in experimental Design, p 517) # - Likelihood Ratio test # (see e.g. Muirhead, R.J. & Waternaud, C.M. (1980): Asymptotic distributions in canonical # correlation analysis and other multivariate procedures for nonnormal populations. Biometrika, 67) # - Sphericity test for nonnormal multivariate distributions based on the kurtosis # (see e.g. Muirhead, R.J. & Waternaud, C.M. (1980)) # - John's test # (see e.g. Wang, Q. & Yao, J. (2013): On the sphericity test with large-dimensional observations. # Electronic Journal of Statistics, Vol. 7) # - Nagao test # (see e.g. Wang, Q. & Yao, J. (2013) # # data can be supplied either in wide format od in lang format # wide format: # y : dependent variables (data.frame or matrix) # groups: vector of group indices (optional) # e.g. check.sphere(winer518[,3:5],winer518[,2]) # long format: # y : vector of dependent variable # groups: vector of group indices (optional) # trial : vector of repeated measure factor indices # id : vector of case ids # e.g. with(winer518t,check.sphere(score,trial=Zeit,id=Vpn,groups=Geschlecht)) {error <- 0 group <- F if (!is.null(groups)) group <- T if (!is.null(trial) & !is.null(id)) # long format supplied { if (group) df<-data.frame(Id=id,g=groups,y=y,w=trial) else df<-data.frame(Id=id,y=y,w=trial) dfw<-reshape(df,direction="wide",idvar="Id",timevar="w") if (group) {nt<-(dim(dfw)[2]-1)/2 g<-dfw[,2] y<-dfw[,seq(3,2*nt+2,by=2)] } else {nt<-dim(dfw)[2]-1 y<-dfw[,seq(2,nt+1)] } } else g <- groups nv <- dim(y)[2] # n of replications nc <- dim(y)[1] # n of cases if (group) {lg <- length(g) if (lg==nc) {ni <- table(g) # cell counts ng <- length(ni) } # n of groups else error<-2 } else ng <- 1 ano <- data.frame(v1=1:9,v2=1:9,v3=1:9) names(ano) <- c("statistic","df","p value") # covariance matrix xm <- apply(y,2,mean) S <- matrix(0,nv,nv) for (i in 1:nc) for (j1 in 1:nv) for (j2 in 1:j1) S[j1,j2] <- S[j1,j2]+ (y[i,j1]-xm[j1])*(y[i,j2]-xm[j2]) for (j1 in 1:(nv-1)) for (j2 in (j1+1):nv) S[j1,j2] <- S[j2,j1] S <- S/(nc-1) if (group) { xmg <- matrix(0,ng,nv) for (j in 1:nv) xmg[,j] <- tapply(y[,j],g,mean) Sp <- matrix(0,nv,nv) # pooled covariance matrix Sg <- array(0,c(ng,nv,nv)) # groupwise covariance matrices for (i in 1:nc) {ig <- as.integer(g[i]) for (j1 in 1:nv) for (j2 in 1:j1) {Sp[j1,j2] <- Sp[j1,j2] + (y[i,j1]-xmg[ig,j1])*(y[i,j2]-xmg[ig,j2]) Sg[ig,j1,j2] <- Sg[ig,j1,j2]+ (y[i,j1]-xmg[ig,j1])*(y[i,j2]-xmg[ig,j2]) } } for (j1 in 1:(nv-1)) for (j2 in (j1+1):nv) { Sp[j1,j2] <- Sp[j2,j1] # fill up upper triangle for (ig in 1:ng) Sg[ig,j1,j2] <- Sg[ig,j2,j1] } Sp <- Sp/(nc-ng) Ss <- matrix(0,nv,nv) for (ig in 1:ng) {Ss <- Ss + Sg[ig,,] Sg[ig,,] <- Sg[ig,,]/(ni[ig]-1) } Ss<-Ss/(ng) S <- Sp # use pooled matrix in case of grouping } sii <- 0 # variances for (j in 1:nv) sii[j]<-S[j,j] # compound symmetry matrix xvar <- 0 ; xcov <- 0 for (j1 in 1:nv) xvar <- xvar+S[j1,j1] for (j1 in 2:nv) for (j2 in 1:(j1-1)) xcov <- xcov+S[j1,j2] xcov <- xcov/(nv*(nv-1)/2) Scs <- matrix(xcov,nv,nv) for (j1 in 1:nv) Scs[j1,j1] <- xvar/nv # univariate kurtosis kur.v <- 0 # for each variable for (j in 1:nv) {kur.v[j]<-0 for (i in 1:nc) kur.v[j]<-kur.v[j]+(y[i,j]-xm[j])^4 kur.v[j]<-(kur.v[j]/nc)/(3*S[j,j]^2) } # kurt=1 for normal distr kurt.u <- mean(kur.v) # mean of all kurtoses kurt2 <- kurt.u - 1 # deviation from normal C <- t(contr.helmert(nv)) # normalized orthogonal transformation for (j in 1:(nv-1)) {cc<-sqrt(sum(C[j,]^2)) ; C[j,]<-C[j,]/cc } S1 <- C%*%S%*%t(C) S2 <- S1%*%S1 p <- nv-1 # Box epsilon ev <- eigen(S1)$values U <- (sum(diag(S2)))/(sum(diag(S1)))^2 eps <- 1/(p*U) # John's test V<- (nc*p*p/2)*(U - 1/p) pjon<-1-pchisq(V/(1+kurt2),p*(p+1)/2-1) # chisq version + m4 correction ano[5,] <- c(V/(1+kurt2),p*(p+1)/2-1,pjon) row.names(ano)[5] <- "John's test chisq" mx <- (nv-1)/2 # exact beta version if (nv>3) mx <- mx - (nv-2)*(nv-3)/((nc+1)*nv) m1 <- mx*(nv+2) m2 <- mx*nv*(nc-2) T <- ((nv-1)*U-1)/(nv-1) pbet <- 1-pbeta(T,m1/2,m2/2) ano[6,] <- c(T,m1,pbet) row.names(ano)[6] <- "John's test beta" # Nagao T <- V f <- nv*(nv+1)/2-1 ap <- ( nv^3 +3*nv^2 -8*nv -12 -200/nv)/12 # chi**2 approximation bp <- (-2*nv^3 -5*nv^2 +7*nv +12 +420/nv)/8 cp <- ( nv^3 +2*nv^2 -nv -2 -216/nv)/4 dp <- (-2*nv^3 -3*nv^2 +nv +436/nv)/24 P <- pchisq(T,f) + (ap*pchisq(T,f+6) + bp*pchisq(T,f+4) + cp*pchisq(T,f+2) + dp*pchisq(T,f) )/nc pnag <- 1-P ano[7,] <- c(T,f,max(0,pnag)) row.names(ano)[7] <- "John-Nagao" # Mauchly test (as in Winer) detS1<- prod(ev) trS1 <- sum(diag(S1)) logW <- -(nc-1)*(log(detS1)-p*log(trS1/p)) d <- 1 - (2*p^2 + p + 2)/(6*p*(nc-1)) L <- d*logW df <- p*(p+1)/2-1 pmau <- 1-pchisq(L,df) ano[1,] <- c(L,df,pmau) row.names(ano)[1] <- "Mauchly test" # LR test plr <- 1-pchisq(logW,(nv-1)*(nv+2)/2) ano[4,] <- c(logW,(nv-1)*(nv+2)/2,plr) row.names(ano)[4] <- "Likelihood Ratio" # compound symmetry test (as in Winer) M2 <- -(nc-ng)*log(prod(eigen(S)$values)/prod(eigen(Scs)$values)) C2 <- nv*(nv+1)*(nv+1)*(2*nv-3)/(6*(nc-ng)*(nv-1)*(nv*nv+nv-4)) mchi<- (1-C2)*M2 pmgr<- 1-pchisq(mchi,(nv*nv+nv-4)/2) ano[9,] <- c(mchi,(nv*nv+nv-4)/2,pmgr) row.names(ano)[9] <- "compound symmetry" # multisample Mauchly test (Mendoza) if (group) { n1 <- (nc-ng) A1 <- n1*nv*log(n1)/2 A2 <- n1*nv*log(sum(diag(C%*%Ss%*%t(C)))/nv)/2 for (ig in 1:ng) { A1 <- A1 + (ni[ig]-1)*log(prod(eigen(C%*%Sg[ig,,]%*%t(C))$values))/2 A2 <- A2 + (ni[ig]-1)*nv*log(ni[ig]-1)/2 } mapp1 <- (n1*nv^2*(nv+1)*(2*nv+1)-2*n1*nv^2)*sum(1/(ni-1))-4 # chi**2 approximation mapp2 <- 6*n1*nv*(ng*nv*(nv+1)-2) chime <- -2*(A1-A2)*(1-mapp1/mapp2) pmend<- 1-pchisq(chime,ng*nv*(nv+1)/2-1) ano[2,] <- c(chime,ng*nv*(nv+1)/2-1,pmend) row.names(ano)[2] <- "multisample Mauchly" # multisample Mauchly test Harris A1 <- 0 ; A2 <- 0 for (ig in 1:ng) { A <- C%*%Sg[ig,,]%*%t(C) A1 <- A1 + ni[ig]*sum(diag(A%*%A)) A2 <- A2 + ni[ig]*sum(diag(A)) } W <- (nc*p/2)*(nc*(nv-1)*A1/(A2^2)-1) pharr<- 1-pchisq(W,ng*nv*(nv-1)/2-1) ano[3,] <- c(W,ng*nv*(nv-1)/2-1,pharr) row.names(ano)[3] <- "multisample Harris" } # Muirhead & Waternaud test, using kurtosis chi2 <- logW/(1+kurt2) pmw2 <- 1-pchisq(chi2,(nv-1)*(nv+2)/2) ano[8,] <- c(chi2,(nv-1)*(nv+2)/2,max(0,pmw2)) row.names(ano)[8] <- "Muirhead & Waternaud" if (!group) ano<-ano[-(2:3),] # ano <- round(ano,digits=4) list(results=ano, Box.epsilon=eps, k.deviation=kurt2, error=error ) }