## Theorem 6.6(a)(for balanced data and constant p): call the functions of theorem 5.6 (a)(b)(c)(d)

## (a)+Modified(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)
  N <- nrow(data)
  m <- ncol(Z)
  G <- as.numeric(data[N, n+1])

  ## theorem5.6(a)
  if(p==0)
    delta <- REMLE.para(data[,1:n], Z, p, REMLE=FALSE)$cov.p

  ## theorem5.6(d)
  else
    delta <- REMLE.phi.delta(data[,1:n], Z, p, REMLE=FALSE)

  RSS <- 0
  sum <- 0   
  for(i in 1:n){   

    delta1 <- N*log(delta[i,i])   ## RSS/N=delta(MLE)
      
    delta2 <- 0
    sum1 <- 0
    for(g in 1:G){
      data1 <- data[data[, n+1]==g,]   
      Z1 <- Z[data[, n+1]==g,]
      N1 <- nrow(data1)  #N(g)

      Z1.g <- NULL
      for(k in 1:m){
        if(sum(Z1[,k])!=0)
          Z1.g <- cbind(Z1.g, Z1[,k])
      }

      if(p==0) 
        delta2 <- delta2+N1*log(REMLE.para(data1[,1:n], Z1.g, p, REMLE=FALSE)$cov.p[i,i])  ## thm5.6(a) RSS/N=delta(MLE)
      else
        delta2 <- delta2+N1*log(REMLE.phi.delta(data1[,1:n], Z1.g, p, REMLE=FALSE)[i,i])  ## thm5.6(d) RSS/N=delta(MLE)
    
      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)    
}



####################################################################################
## the following functions are of theorem 5.6 (a)(b)(c)(d)

###############Theorem 5.6(a)###############
## set REMLE=TRUE for REML estimation
## set REMLE=FALSE for ML estimation

REMLE.para<-function(data, covariate, p, REMLE){
  data <- as.matrix(data)
  Z <- as.matrix(covariate)
  N <- nrow(data)
  n <- ncol(data)
  m <- ncol(Z)
  V <- matrix(rep(0, n*n), nrow = n)

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

  if(m==1){
    beta.hat <- as.vector(apply(data,2, mean))
    A <- as.matrix(cov(data))*(N-m)/N
  }

  else{
    Y <- NULL
    for(i in 1:N)
      Y <- cbind(Y, t(data[i,]))
    Y <- t(Y) 
    A <- matrix(rep(0, n*n), nrow = n)
 
    ## (5.30) p.134
    beta.hat <- kronecker(solve(t(Z)%*%Z)%*%t(Z), diag(n))%*%Y
  
    for(i in 1:N){
      M <- matrix(data[i,], nr=n)-kronecker(matrix(Z[i,], nr=1), diag(n))%*%beta.hat
      A <- A+M%*%t(M)
    }
    A <- A/N 
  }

  ## (5.31, 5.32)
  if(REMLE==TRUE)  cov <- (N/(N-m))*A       
  else  cov <- A 
         
  for(i in 1:n){
    for(j in (i-p[i]):i)
      V[j,i] <- V[i,j] <- cov[i,j]
  }

  if(m==1)   
    list(muhat=beta.hat, cov.p=V)
  else   
    list(betahat=beta.hat, cov.p=V)
}



###############Theorem 5.6(b)###############   

## calculate the correlations on off diagonals
cov.cor<-function(V){
  for(i in 2: nrow(V)){
    for(j in 1:(i-1))
      V[j,i] <- V[i,j] <- V[i,j]/sqrt(V[i,i]*V[j,j])
  }
  V
}

## recursive formula for mle in(2.32)
recur.mle<-function(i, p, cov){
  kasi <- matrix(t(cov[c((i-p):(i-1)), c(1:(i-p-1))]), nr=(i-p-1)) %*%
    solve(t(cov[c((i-p):(i-1)), c((i-p):(i-1))])) %*% matrix(t(cov[i, c((i-p):(i-1))]), nr=p)
  kasi
}

REMLE.cov<-function(data, covariate, p, REMLE){
  n <- ncol(data)

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

  cov <- REMLE.para(data, covariate, p, REMLE=REMLE)$cov.p

  for(i in 2:n){
    if((i-p[i]) > 1 && (i-p[i]) < i){
      kasi <- recur.mle(i,p[i], cov) ## recursive formula in (2.44)
      for(j in 1:(i-p[i]-1))
        cov[j,i] <- cov[i,j] <- as.vector(kasi[j])
    }
  }     
  cov
}


###############Theorem 5.6(c)###############

## inverse
inv <- function(cov){
  solve(cov)
}

## pi>0 for i>1
## intervenor-adjusted partial correlation

## recursion formula for partial correlations in (2.6)
recur.pcor<-function(index1, index2, index3, V){
  index1 <- index1
  index2 <- index2
  index3 <- index3
  index <- c(index1, index2, index3)

  if(length(index3)==1){
    rxy.z <- (V[index[1], index[2]]-V[index[1], index[3]]*V[index[2], index[3]])/
            (sqrt(1-V[index[1], index[3]]^2)*sqrt(1-V[index[2], index[3]]^2))
    return(rxy.z)
  }
  else{
    index1 <- index[1]
    index2 <- index[2]
    index30 <- index[3]    
    index3c <- index[-c(1,2,3)]

    rxy.zc <- recur.pcor(index1, index2, index3c, V)
    rxz0.zc <- recur.pcor(index1,index30,index3c, V)
    ryz0.zc <- recur.pcor(index2,index30,index3c, V)
    rxy.z <- (rxy.zc - rxz0.zc*ryz0.zc)/(sqrt(1-rxz0.zc^2)*sqrt(1-ryz0.zc^2))
    return(rxy.z)
  }                       
}  

partial.inter<-function(cor){
  V <- cor
  n <- ncol(cor)
  m <- n-2
  for(k in 1:m){
    if(k==m)
       V[1, n] <- V[n, 1] <- recur.pcor(1, n, c(2:(n-1)), cor)

     else
       diag(V[1:(n-k-1), (k+2):n]) <- diag(V[(k+2):n, 1:(n-k-1)]) <- sapply(1:(m+1-k), function(i) if(k==1) 
                    recur.pcor(i, i+k+1, i+1, cor) else recur.pcor(i, i+k+1, c((i+1):(i+k)), cor))
  }
  V
}


###############Theorem 5.6(d)###############
## pi>0 for i>1

REMLE.phi.delta<-function(data, covariate, p, REMLE){
  N <- nrow(data)
  n <- ncol(data)
  V <- matrix(rep(0, n*n), nrow = n)
  T <- matrix(rep(0, n*n), nrow = n)
  diag(T) <- rep(1,n)
  D <- matrix(rep(0, n*n), nrow = n)
  A <- matrix(rep(0, n*n), nrow = n)

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

  cov <- REMLE.para(data, covariate, p, REMLE=REMLE)$cov.p

  D[1,1] <- cov[1,1]
  A[1,1] <- solve(T[1,1])*D[1,1]*solve(T[1,1])
  V[1,1] <- D[1,1]    

  for(i in 2:n){ 
    phi <- matrix(t(cov[(i-p[i]):(i-1),i]), nr=1) %*% solve(A[(i-p[i]):(i-1), (i-p[i]):(i-1)])
    delta <- cov[i,i]-matrix(t(cov[(i-p[i]):(i-1), i]), nr=1) %*% 
       solve(A[(i-p[i]):(i-1), (i-p[i]):(i-1)]) %*% matrix(cov[(i-p[i]):(i-1), i], nc=1)
    V[i,i] <- delta
    V[i, (i-p[i]):(i-1)] <- phi
    T[i, (i-p[i]):(i-1)] <- -phi
    D[i,i] <- delta      
    if(i!=n) 
      A <- solve(T[(1:i), (1:i)]) %*% D[(1:i), (1:i)] %*% solve(t(T[(1:i), (1:i)]))
  }
  V
}