R/predictContinuous.R

In CopulaGAMM: Copula-Based Mixed Regression Models

#' @title Conditional expectation for a copula-based estimation of mixed regression models for continuous response
#' @description Compute the conditional expectation of a copula-based  2-level hierarchical model for continuous response.
#' @param object  Object of class ``EstContinuous`` generated by EstContinuous.
#' @param newdata List of variables for be predicted (``clu`` for clusters, ``xc`` for the copula covariates, and ``xm`` for the margins covariates). The covariates can be NULL.
#' @param nq      number of nodes and weighted for Gaussian quadrature of the product of conditional copulas; default is 25.
#'
#' @return \item{mest}{Conditional expectations}
#'
#' @references Krupskii, Nasri & Remillard (2023). On factor copula-based mixed regression models
#' @author Pavel Krupskii and Bruno N. Remillard, January 20, 2023
#' @import statmod, matrixStats
#' @examples
#' data(out.normal)
#' newdata=list(clu=c(1:50),xm=rep(0.4,50))
#' pred= predictContinuous(out.normal,newdata)
#' @export

predictContinuous=function(object,newdata=NULL,nq=25)
{

  family  = object$family
  thC0    = object$thC0
  dfC     = object$dfC
  V       = object$V
  par     = object$coefficients
  cluster = object$cluster
  disc    = object$disc
  rot     = object$rot
  model   = object$model
  dfM     = object$dfM

  gl=statmod::gauss.quad.prob(nq)
  nl=gl$nodes
  wl=gl$weights

  fun1 = function(q,a){q*a}
  fun2 = function(u,a){ qweibull(u,a)}



  if(is.null(newdata))
  {

    thF  = object$thF
    clu  = object$clu
    switch(model,
           "normal"  = {
             ql = qnorm(nl)
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
             mu = thF[,1]
           },

           "t"       = {
             ql=qt(nl,dfM)
             mu = thF[,1]
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
           },
           "laplace" = {
             ql=qlap(nl)
             mu = thF[,1]
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
           },
           "exponential" = {
             ql=qexp(nl)
             thF = cbind(0,thF)
             mu = thF[,1]
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
           },
           "weibull"     ={
             thF = cbind(0,thF)
             mu = thF[,1]
             sigma = thF[,3]
             shape=thF[,2]
             M0=outer(nl,shape,FUN=fun2)
             dd=dim(M0)
             sig=rep(sigma,dd[1])
             msig = matrix(sig,nrow=dd[1],byrow=TRUE)
             M = M0*msig}
    )

    nx=length(clu)
    mest=rep(0,nx)
    for(i in 1:nx)
    {
      k = which(cluster ==clu[i] )
      mest[i] = mu[i] + sum(wl*M[,i]*dcop(nl,V[k],family,rot,thC0[i],dfC))
    }
  }else{

    clu  = newdata$clu
    nx   = length(clu)
    xc   = newdata$xc
    xm   = newdata$xm


    if(is.null(xc))
    {Matxc = matrix(1,nrow=nx,ncol=1)}else
    {Matxc = cbind(1,xc)}

    if(is.null(xm))
    {Matxm = matrix(1,nrow=nx,ncol=1)}else
    {Matxm = cbind(1,xm)}

    k1 = ncol(Matxc)
    k2 = ncol(Matxm)
    nx = nrow(Matxm)

    thC  = par[1:k1]

    thF = colSums(par[k1+(1:k2)]*t(Matxm))
    thF = cbind(thF, par[k1+k2+1])
    switch(model,
           "normal"  = {
             ql = qnorm(nl)
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
             mu = thF[,1]
           },

           "t"       = {
             ql=qt(nl,dfM)
             mu = thF[,1]
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
           },
           "laplace" = {
             ql=qlap(nl)
             mu = thF[,1]
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
           },
           "exponential" = {
             ql=qexp(nl)
             thF = cbind(0,thF)
             mu = thF[,1]
             sigma = thF[,2]
             M=outer(ql,sigma,FUN=fun1)
           },
           "weibull"     ={
             thF = cbind(0,thF)
             mu = thF[,1]
             sigma = thF[,3]
             shape=thF[,2]
             M0=outer(nl,shape,FUN=fun2)
             dd=dim(M0)
             sig=rep(sigma,dd[1])
             msig = matrix(sig,nrow=dd[1],byrow=TRUE)
             M = M0*msig}
    )



    mest=rep(0,nx)
    for(i in 1:nx)
    {
      k = clu[i]
      Matxck = Matxc[i,]
      thCk = sum(thC*Matxck)

      thC0 = linkCop(thCk,family)$cpar



      mest[i] = mu[i] + sum(wl*M[,i]*dcop(nl,V[k],family,rot,thC0,dfC))

    }
  }

  return(mest)

}