## theorem 7.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.covari<-function(data, covariate1, covariate2, p){
  data <- as.matrix(data)
  Z1 <- as.matrix(covariate1)
  Z2 <- as.matrix(covariate2)
  n <- ncol(data)
  m0 <- ncol(Z1)
  m <- ncol(Z2)

  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){
    Z1 <- as.matrix(Z1[!is.na(data[,i]),])    
    Z2 <- as.matrix(Z2[!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 <- Z1
      C <- Z2
    }
    else{
      B <- cbind(Z1, as.matrix(data[, (i-min(p[i], i-1)):(i-1)]) )
      C <- cbind(Z2, 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)))

    x <- m-m0 
    y <- N-m-min(p[i], i-1)
    sum.M <- sum.M + phi(x, y)

  }
  LRT <- RSS
  MLRT <- n*(m-m0)*RSS.M/sum.M

  df <- n*(m-m0)
  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)    
}




race100k<-read.table("race100k.dt")
## delete missing data
race100k<- race100k[-c(3, 10, 50, 55),]  ## 76 competitors
race<-race100k[, 2:11]
## centered age
age.sq<-cbind(rep(1,nrow(race100k)),(race100k[,12]- mean(race100k[,12])),(race100k[,12]- mean(race100k[,12]))^2)
age.lin<-cbind(rep(1,nrow(race100k)),(race100k[,12]- mean(race100k[,12])))
age.sat<-rep(1,nrow(race100k))

p<-c(0,1,1,1,2,1,1,2,3,5)

##table 7.3(matched)
> test.R.covari(race, age.sat, age.sq, p)
$LRT
[1] 34.74391

$pvalue.LRT
[1] 0.02151141

$MLRT
[1] 32.58337

$pvalue.MLRT
[1] 0.03746529

> test.R.covari(race, age.sat, age.lin, p)
$LRT
[1] 18.28406

$pvalue.LRT
[1] 0.05035684

$MLRT
[1] 17.26657

$pvalue.MLRT
[1] 0.06866969

> test.R.covari(race, age.lin, age.sq, p)
$LRT
[1] 16.45984

$pvalue.LRT
[1] 0.08720385

$MLRT
[1] 15.3272

$pvalue.MLRT
[1] 0.1205803