## theorem 7.1(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, mu0, p){
  data <- as.matrix(data)
  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){
    data <- data[!is.na(data[,i]),]      #delete missing data
    N <- nrow(data)  #Ni

    A1 <- data[, i]-mu0[i]
    A2 <- data[, i]

    if(i==1 || p[i]==0){
      RSS1 <- log(sum(A1^2))
      RSS2 <- log((summary(lm(A2~1))[6]$sigma)^2*(N-1-min(p[i],i-1)))  
    }

    else{  
      B1 <- as.matrix(data[, (i-min(p[i], i-1)):(i-1)])-matrix(rep(c(mu0[(i-min(p[i], i-1)):(i-1)]), N), nrow=N, byrow=TRUE)
      B2 <- as.matrix(data[, (i-min(p[i], i-1)):(i-1)])    
      RSS1 <- log((summary(lm(A1~B1-1))[6]$sigma)^2*(N-min(p[i],i-1)))  
      RSS2 <- log((summary(lm(A2~B2))[6]$sigma)^2*(N-1-min(p[i],i-1))) 
    }
    RSS <- RSS + N*(RSS1-RSS2)
    RSS.M <- RSS.M + (RSS1-RSS2)

    y <- N-1-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)    
}



##“Race Data”
race100k<-read.table("race100k.dt")  ##complete set of 80 competitors
race<-race100k[, 2:11]
mu0<-c(48, 51, 50, 54, 54, 60, 63, 68, 68, 68)

##(p.202) matched
> test.R.mean(race, mu0, c(0,1,1,1,2,1,1,2,3,5))
$LRT
[1] 20.46500

$pvalue.LRT
[1] 0.02514893

$MLRT
[1] 19.6394    

$pvalue.MLRT
[1] 0.03285387