## theorem 7.3(the code works for balanced\monotone missing data & constant\variable order p) phi <- function(x, y){ if(x == 0) r <- 0 else if(x == 1) r <- (2*y-1)/(2*y*(y-1)) else if(x > 0 && x%%2 == 0){ s <- 0 for(l in 1: (x/2)) s <- s + 1/(y+2*l-2) r <- 2*s } else if(x > 1 && x%%2 == 1) r <- phi(1, y) + phi(x-1, y+1) r } test.R.mean<-function(data, covariate, p){ data <- as.matrix(data) Z <- as.matrix(covariate) n <- ncol(data) if(length(p)==1) p <- rep(0:p, c(rep(1,p),(n-p))) RSS <- 0 RSS.M <- 0 sum.M <- 0 for(i in 1:n){ Z <- as.matrix(Z[!is.na(data[,i]),]) data <- data[!is.na(data[,i]),] #delete missing data N <- nrow(data) #Ni A <- data[, i] if(i==1 || p[i]==0){ B <- rep(1,N) C <- Z } else{ B <- as.matrix(data[, (i-min(p[i], i-1)):(i-1)]) C <- cbind(Z, as.matrix(data[, (i-min(p[i], i-1)):(i-1)]) ) } RSS <- RSS + N*(log(var(residuals(lm(A~B)))*(N-1)) - log(var(residuals(lm(A~C)))*(N-1))) RSS.M <- RSS.M + (log(var(residuals(lm(A~B)))*(N-1)) - log(var(residuals(lm(A~C)))*(N-1))) y <- N-2-min(p[i], i-1) sum.M <- sum.M + phi(1, y) } LRT <- RSS MLRT <- n*RSS.M/sum.M df <- n pvalue.LRT <- 1-pchisq(LRT, df) pvalue.MLRT <- 1-pchisq(MLRT, df) list(LRT=LRT, pvalue.LRT=pvalue.LRT, MLRT=MLRT, pvalue.MLRT=pvalue.MLRT) } ##ex of monotone missing data speech<-read.table("cochlear09.dt") covariate<-rbind(t(matrix(rep(c(1,0),20), nrow=2)), t( matrix(rep(c(0,1),21), nrow=2))) ##Table 7.2(matched for Ho: Profiles equal ) > test.R.mean(speech[,2:5], covariate, 1) $LRT [1] 5.921606 $pvalue.LRT [1] 0.2050798 $MLRT [1] 6.208932 $pvalue.MLRT [1] 0.1840790 > test.R.mean(speech[,2:5], covariate, c(0,1,1,1)) $LRT [1] 5.921606 $pvalue.LRT [1] 0.2050798 $MLRT [1] 6.208932 $pvalue.MLRT [1] 0.1840790 > test.R.mean(speech[,2:5], covariate, c(0,1,0,3)) $LRT [1] 7.304985 $pvalue.LRT [1] 0.1206225 $MLRT [1] 7.3772 $pvalue.MLRT [1] 0.1172478