R/LRTval.R

In elara7/VCM: A R package for the paper written by HAIQIANG CHEN, YING FANG, AND YINGXING LI

#' Testing whether need to use varying coefficience
#'
#' @title Testing whether need to use varying coefficience
#' @param X predictors, a vector
#' @param Y dependent variable, a vector
#' @param Z assist variable, a vector
#' @param p order of B-splines, p = 1 by default
#' @param nkonts the number of knots
#' @param df the numbers of constraint of fixed coefficient, df = 1 by default, can't be larger than p
#' @param sigma.output covariance matrix of disturbance term u
#' @export

LRTval=function(X, Z, Y, p=1,nknots=10,df = 1,sigma.output=diag(rep(1, length(Y)))){
  require(lme4)
  require(RLRsim)
  require(nlme)
  if(df > p){
    stop("df should not larger than p")
  }
  n=length(Y)
  x1 <- diag(X)
  z1 <- rep(1,n)
  for (i in 1:p){z1=cbind(z1, z^i)}
  myknots = quantile(unique(z),seq(0,1,length=(nknots+2))[-c(1,(nknots+2))])
  z2 <- outer(z, myknots, FUN="-")
  z3 <- z2*(z2>0)
  if(p>1) {z3=z3^p}
  A1 <- x1 %*% z1
  A2 <- x1 %*% z3
  Xnames = paste("X",1:ncol(A1),sep="")
  Znames = paste("Z",1:ncol(A2),sep="")
  fixed.model = as.formula(paste("DATA.temp ~ -1+",
                                 paste(paste("X",1:ncol(A1),sep=""),collapse="+")))
  fixed.model2 = as.formula(paste("DATA.temp ~ -1+",
                                  paste(paste("X",1:(ncol(A1)-df),sep=""),collapse="+")))
  random.model = as.formula(paste("~-1+",paste(paste("Z",1:ncol(A2),sep=""),collapse="+")))

  DATA.output = as.vector( sigma.output  %*% Y )
  hat1 = sigma.output %*% A1 ;
  hat2 = sigma.output %*% A2

  DATA.temp= DATA.output
  colnames(hat1)=Xnames
  colnames(hat2)=Znames
  subject<-rep(1,n)
  ALLDATA  = data.frame(cbind(subject,DATA.temp, hat1, hat2))
  mA = lme(fixed=fixed.model, data=ALLDATA,
           random=list(subject = pdIdent(random.model)), method="ML")
  m0 = lm(fixed.model2, data=ALLDATA)
  obs.LRT = as.numeric(2*(logLik(mA)-logLik(m0)))
  pvalue = pchisq(obs.LRT,1,lower.tail = FALSE)
  res <- list(LRT =obs.LRT,p.value = pvalue)
  res
}