## Theorem 6.6(a)(for monotone missing data & constant p) ## (a)+Modified LRT(6.13) 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.homo<-function(data, covariate, p){ data <- as.matrix(data) Z <- as.matrix(covariate) n <- ncol(data)-1 #last column of the data is the class level (1~G) m <- ncol(Z) G <- as.numeric(data[nrow(data), n+1]) RSS <- 0 sum <- 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==0) B <- Z else B <- cbind(Z, as.matrix(data[, (i-min(p, i-1)):(i-1)]) ) delta1 <- N*log(var(residuals(lm(A~B)))*(N-1)/N) delta2 <- 0 sum1 <- 0 for(g in 1:G){ data1 <- data[data[, n+1]==g,] Z1 <- Z[data[, n+1]==g,] N1 <- nrow(data1) #Ni(g) Z1.g <- NULL for(k in 1:m){ if(sum(Z1[,k])!=0) Z1.g <- cbind(Z1.g, Z1[,k]) } A <- data1[, i] if(i==1 || p==0) B <- Z1.g else B <- cbind(Z1.g, as.matrix(data1[, (i-min(p, i-1)):(i-1)]) ) delta2 <- delta2+N1*log(var(residuals(lm(A~B)))*(N1-1)/N1) x <- N-N1-G+1 y <- N1-min(p, i-1)-1 sum1 <- sum1 + N1*(phi(x, y)-log(N/N1)) } sum <- sum+sum1 RSS <- RSS+(delta1-delta2) } LRT <- RSS MLRT <- (((G-1)*(2*n-p)*(p+1)/2)/sum)*RSS ## (6.13) df <- (G-1)*(2*n-p)*(p+1)/2 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) }