R/clprobit.R

In lava.tobit: Latent Variable Models with Censored and Binary Outcomes

##' Composite Likelihood for probit latent variable models
##'
##' Estimate parameters in a probit latent variable model via a composite
##' likelihood decomposition.
##' @param x \code{lvm}-object
##' @param data data.frame
##' @param k Size of composite groups
##' @param type Determines number of groups. With \code{type="nearest"} (default)
##' only neighboring items will be grouped, e.g. for \code{k=2}
##' (y1,y2),(y2,y3),... With \code{type="all"} all combinations of size \code{k}
##' are included
##' @param pairlist A list of indices specifying the composite groups. Optional
##' argument which overrides \code{k} and \code{type} but gives complete
##' flexibility in the specification of the composite likelihood
##' @param silent Turn output messsages on/off
##' @param \dots Additional arguments parsed on to lower-level functions
##' @return An object of class \code{clprobit} inheriting methods from \code{lvm}
##' @author Klaus K. Holst
##' @seealso \link{estimate}
##' @keywords models regression
##' @export
clprobit <- function(x,data,k=2,type=c("nearest","all"),pairlist,silent=TRUE,
                     ...) {
  y <- lava::endogenous(x)
  binsurv <- rep(FALSE,length(y))
  for (i in 1:length(y)) {
    z <- data[,y[i]]
    ## binsurv[i] <- is.Surv(z) | (is.factor(z) && length(levels(z))==2)
    binsurv[i] <- survival::is.Surv(z) | (is.factor(z))
  }

  binsurv <- unique(c(y[binsurv],binary(x),lava::ordinal(x)))
  ##  binsurvpos <- which(colnames(data)%in%binsurv)
  if (!missing(pairlist)) {
    binsurvpos <- which(colnames(data)%in%lava::endogenous(x))
  } else {
    binsurvpos <- which(colnames(data)%in%binsurv)
  }
  
  if (missing(pairlist)) {
    if (length(binsurv)<(k+1)) stop("No need for composite likelihood analysis.")

    if (type[1]=="all") {
      mypar <- utils::combn(length(binsurv),k) ## all pairs (or multiplets), k=2: k*(k-1)/2
    } else {
      mypar <- sapply(0:(length(binsurv)-k), function(x) x+1:k)
    }
  } else {
    mypar <- pairlist
  }  
  
  if (is.matrix(mypar)) {
    mypar0 <- mypar; mypar <- c()
    for (i in seq(ncol(mypar0)))
      mypar <- c(mypar, list(mypar0[,i]))
  }
  
  nblocks <- length(mypar)
  mydata0 <- data[c(),,drop=FALSE]  
  mydata <-  as.data.frame(matrix(NA, nblocks*nrow(data), ncol=ncol(data)))
  names(mydata) <- names(mydata0)
  for (i in 1:ncol(mydata)) {
    if (is.factor(data[,i])) {
      mydata[,i] <- factor(mydata[,i],levels=levels(mydata0[,i]))
    }
    if (survival::is.Surv(data[,i])) {
      S <- data[,i]
      for (j in 2:nblocks) S <- rbind(S,data[,i])
      S[,1] <- NA
      mydata[,i] <- S
    }
  }
  
  for (ii in 1:nblocks) {    
    data0 <- data;
    for (i in binsurvpos[-mypar[[ii]]]) {
      if (survival::is.Surv(data[,i])) {
        S <- data0[,i]; S[,1] <- NA
        data0[,i] <- S
      } else {
        data0[,i] <- NA
        if (is.factor(data[,i])) data0[,i] <- factor(data0[,i],levels=levels(data[,i]))
      }
    }
    mydata[(1:nrow(data))+(ii-1)*nrow(data),] <- data0
    ##mydata <- rbind(mydata,data0)
  }

  ## suppressMessages(browser())
  suppressWarnings(e0 <- lava::estimate(x,data=mydata,missing=TRUE,silent=silent,
              ...))

  S <- lava::score(e0,indiv=TRUE)
  nd <- nrow(data)
  block1 <- which((1:nd)%in%(rownames(S)))
  blocks <- sapply(1:nblocks, function(x) 1:length(block1)+length(block1)*(x-1))
  Siid <- matrix(0,nrow=length(block1),ncol=ncol(S))
  for (j in 1:ncol(blocks)) {
    Siid <- Siid+S[blocks[,j],]
  }
  iI <- stats::vcov(e0);
  J <- t(Siid)%*%(Siid)
  e0$iidscore <- Siid
  e0$blocks <- blocks
  e0$vcov <- iI%*%J%*%iI ## thetahat-theta0 :=(asymp) I^-1*S => var(thetahat) = iI*var(S)*iI 
  cc <- e0$coef; cc[,2] <- sqrt(diag(e0$vcov))
  cc[,3] <- cc[,1]/cc[,2]; cc[,4] <- 2*(1-stats::pnorm(abs(cc[,3])))
  e0$coef <- cc
  class(e0) <- c("clprobit",class(e0))
  return(e0)
}

score.clprobit <- function(x,indiv=FALSE,...) {
  if (!indiv)
    return(colSums(x$iidscore))
  x$iidscore
}