## theorem 6.5(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.order<- function(data, covariate, p0, p1){
  data <- as.matrix(data)
  Z <- as.matrix(covariate)
  n <- ncol(data)
  m <- ncol(Z)

  if(length(p0)==1)  p0 <- rep(0:p0, c(rep(1,p0),(n-p0)))
  if(length(p1)==1)  p1 <- rep(0:p1, c(rep(1,p1),(n-p1)))
    
  p <- p0
  q <- p1-p0

  RSS <- 0
  RSS.M <- 0
  sum.M1 <- 0
  sum.M2 <- 0 
  for(i in 2: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(p[i]==0)
      B <- Z
    else
      B <- cbind(Z, as.matrix(data[, (i-min(p[i], i-1)):(i-1)]) )

    if((p[i]+q[i])==0)
      C <- Z
    else
      C <- cbind(Z, as.matrix(data[, (i-min(p[i]+q[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)))

    x <- min(p[i]+q[i], i-1)-p[i] 
    y <- N-m-min(p[i]+q[i], i-1)
    sum.M1 <- sum.M1 + x
    sum.M2 <- sum.M2 + phi(x, y)

  }
  LRT <- RSS
  MLRT <- sum.M1*RSS.M/sum.M2

  df <- sum(q[2: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)    
}