R/sCorrect-iid2.R

In lavaSearch2: Tools for Model Specification in the Latent Variable Framework

### iid2.R --- 
#----------------------------------------------------------------------
## author: Brice Ozenne
## created: okt 12 2017 (13:16) 
## Version: 
## last-updated: jan 18 2022 (09:38) 
##           By: Brice Ozenne
##     Update #: 625
#----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
#----------------------------------------------------------------------
## 
### Code:

## * Documentation - iid2
#' @title  Influence Function With Small Sample Correction.
#' @description  Extract the influence function from a latent variable model.
#' It is similar to \code{lava::iid} but with small sample correction.
#' @name iid2
#'
#' @param object,x a \code{lvmfit} or \code{lvmfit2} object (i.e. output of \code{lava::estimate} or \code{lavaSearch2::estimate2}).
#' @param cluster [integer vector] the grouping variable relative to which the observations are iid.
#' @param robust [logical] if \code{FALSE}, the influence function is rescaled such its the squared sum equals the model-based standard error (instead of the robust standard error).
#' Do not match the model-based correlation though.
#' @param as.lava [logical] if \code{TRUE}, uses the same names as when using \code{stats::coef}.
#' @param ssc [character] method used to correct the small sample bias of the variance coefficients (\code{"none"}, \code{"residual"}, \code{"cox"}). Only relevant when using a \code{lvmfit} object. 
#' @param ... additional argument passed to \code{estimate2} when using a \code{lvmfit} object. 
#'
#' @details When argument object is a \code{lvmfit} object, the method first calls \code{estimate2} and then extract the variance-covariance matrix.
#'
#' @seealso \code{\link{estimate2}} to obtain \code{lvmfit2} objects.
#'
#' @return A matrix containing the 1st order influence function relative to each sample (in rows)
#' and each model coefficient (in columns).
#' 
#' @examples
#' #### simulate data ####
#' n <- 5e1
#' p <- 3
#' X.name <- paste0("X",1:p)
#' link.lvm <- paste0("Y~",X.name)
#' formula.lvm <- as.formula(paste0("Y~",paste0(X.name,collapse="+")))
#'
#' m <- lvm(formula.lvm)
#' distribution(m,~Id) <- Sequence.lvm(0)
#' set.seed(10)
#' d <- sim(m,n)
#'
#' #### latent variable model ####
#' e.lvm <- estimate(lvm(formula.lvm),data=d)
#' iid.tempo <- iid2(e.lvm)
#'
#'
#' @concept extractor
#' @keywords smallSampleCorrection
#' @export
`iid2` <-
  function(object, ...) UseMethod("iid2")

## * iid2.lvmfit
#' @rdname iid2
#' @export
iid2.lvmfit <- function(object, robust = TRUE, cluster = NULL, as.lava = TRUE, ssc = lava.options()$ssc,...){

    return(lava::iid(estimate2(object, ssc = ssc, ...), robust = robust, cluster = cluster, as.lava = as.lava))

}

## * iid2.lvmfit2
#' @rdname iid2
#' @export
iid2.lvmfit2 <- function(object, robust = TRUE, cluster = NULL, as.lava = TRUE, ...){

    dots <- list(...)
    if(length(dots)>0){
        warning("Argument(s) \'",paste(names(dots),collapse="\' \'"),"\' not used by ",match.call()[1],". \n")
    }

    ## ** compute iid
    object.score <- score2(object, indiv = TRUE, cluster = cluster, as.lava = FALSE)
    object.vcov <- vcov2(object, as.lava = FALSE)
    object.iid <- object.score %*% object.vcov

    if(robust == FALSE){
        object.chol.vcov <- matrixPower(object.vcov, power = 1/2, symmetric = TRUE)
        object.chol.rvcov <- matrixPower(crossprod(object.iid[rowSums(!is.na(object.iid))>0,,drop=FALSE]), power = -1/2, symmetric = TRUE)

        object.iid <- object.iid %*% object.chol.rvcov %*% object.chol.vcov
    }

    ## restaure names
    colnames(object.iid) <- colnames(object.score)
    if(!is.null(cluster)){
        attr(object.iid,"cluster") <- attr(object.score,"cluster")
    }

    ## ** export
    if(as.lava){
        object.iid <- object.iid[,names(object$sCorrect$skeleton$originalLink2param),drop=FALSE]
        colnames(object.iid) <- as.character(object$sCorrect$skeleton$originalLink2param)
    }
    return(object.iid)

}

## * iid.lvmfit2
#' @rdname iid2
#' @export
iid.lvmfit2 <- function(x, robust = TRUE, cluster = NULL, as.lava = TRUE, ...){
    return(iid2(x, robust = robust, cluster = cluster, as.lava = as.lava, ...)) ## necessary as first argument of iid must be x 
}

## * Documentation - iid2plot
#' @title  Display the i.i.d. Decomposition
#' @description  Extract the i.i.d. decomposition and display it along with the corresponding coefficient.
#'
#' @param object a \code{lvmfit} or \code{lvmfit2} object (i.e. output of \code{lava::estimate} or \code{lavaSearch2::estimate2}).
#' @param param [character] name of one of the model parameters.
#'
#' @export
iid2plot <- function(object, param){

    all.param <- coef2(object)
    if(length(param) != 1){
        stop("Incorrect argument \'param\'. \n",
             "Should have length 1. \n")
    }
    if(param %in% names(all.param) == FALSE){
        stop("Incorrect argument \'param\'. \n",
             "Should be one of the model parameters. \n")
    }

    object.iid <- iid2(object)[,param]
    h <- graphics::hist(object.iid, main = paste0(param,"=",all.param[param]), xlab = "iid")

    return(invisible(h))
}

##----------------------------------------------------------------------
### iid2.R ends here