#' @title Extract moments from a reduced lvm
#' @description Extract moments from a reduced lvm
#'
#' @param x latent variable model
#' @param p Parameter vector used to calculate statistics
#' @param ... additional arguments to be passed to the low level functions
#'
#' @export
moments.lvmfit.reduced <- function (x, p = pars(x), ...) {
xfull <- reduce2lvm(x)
moments(xfull, p = p[coef(xfull)], ...)
}
#' @title Moments for a reduced lvm model
#' @name momentLVMr
#' @description Compute the likelihood, score and hessian of a reduced lvm model
#'
#' @param x,object \code{lvm.reduced}-object
#' @param p current parameter estimate
#' @param data dataset
#' @param n number of observations
#' @param indiv should the individual contribution be returned? Else average over the observations
#' @param type method used to compute the hessian
#' @param implementation default is R. If set to cpp all the computation of the gradient is made in a C++ routine.
#' @param ... additional arguments
#'
#' @details this function assumes that the external parameters in p are at the end of the vector
#'
#' @rdname momentLVMr
#' @export
gaussianLP_method.lvm <- "nlminb2"
#' @rdname momentLVMr
#' @export
gaussianLP_objective.lvm <- function(x, p, data, ...){
l <- gaussianLP_logLik.lvm(x, p=p, data=data, ...)
return(-l)
}
#' @rdname momentLVMr
#' @export
gaussianLP_logLik.lvm <- function(object, p, data, ...) {
if(is.matrix(p)){p <- as.double(p)}
if(is.null(names(p))){names(p) <- coef(object)}
dataLP <- calcLP(data, p = p,
lp.x = lp(object, type = "x", format = "list"),
lp.link = lp(object, type = "link", format = "list"),
lp.name = lp(object, type = "name"))
## from normal_objective.lvm
y <- lava::index(object)$endogenous
ord <- lava::ordinal(object)
status <- rep(0,length(y))
bin <- tryCatch(match(do.call("binary",list(x=object)),y),error=function(x) NULL)
status[match(ord,y)] <- 2
mu <- stats::predict(object,data=dataLP,p=p) ## can be optimized
S <- attributes(mu)$cond.var
class(mu) <- "matrix"
thres <- matrix(0,nrow=length(y),max(1,attributes(ord)$K-1)); rownames(thres) <- y
for (i in seq_len(length(attributes(ord)$fix))) {
nn <- names(attributes(ord)$idx)[i]
ii <- attributes(ord)$idx[[nn]]
val <- (attributes(mu)$e[ii])
thres[nn,seq_len(length(val))] <-
cumsum(c(val[1],exp(val[-1])))
}
yl <- yu <- as.matrix(data[,y,drop=FALSE])
if (!inherits(yl[1,1],c("numeric","integer","logical")) ||
!inherits(yu[1,1],c("numeric","integer","logical")))
stop("Unexpected data (normal_objective)")
l <- sum(mets::loglikMVN(yl = yl, yu = yu, status = status, mu = mu, S = S, thres = thres))
return(l)
}
#' @rdname momentLVMr
#' @export
gaussianLP_gradient.lvm <- function(x, p, data, ...){
val <- -gaussianLP_score.lvm(x, p = p, data = data, ...)
if (!is.null(nrow(val))){
val <- colSums(val)
}
return(val)
}
#' @rdname momentLVMr
#' @export
gaussianLP_score.lvm <- function(x, p, data, indiv = FALSE, implementation = "R", ...) {
if(is.matrix(p)){p <- as.double(p)}
if(is.null(names(p))){names(p) <- coef(x)}
#### extract information
lp.endo <- lp(x, type = "endogeneous")
lp.link <- lp(x, type = "link", format = "list")
lp.x <- lp(x, type = "x", format = "list")
names.data <- colnames(data)
names.p <- names(p)
M <- moments(x,p) # lava:::moments.lvm
D <- callS3methodParent(x, FUN = "deriv", class = "lvm.reduced",p=p)
if(implementation == "cpp"){
score <- scoreLVM(data = as.matrix(data),
p = p,
mu = M$xi%x%rep(1,NROW(data)),
S = M$C,
dmu = D$dxi,
dS = D$dS,
scoreFun = mets::scoreMVN,
indexCoef = lapply(lp.link, function(link){match(link, names.p) - 1}),
indexEndo = lapply(lp.x, function(link){match(link, names.data) - 1}),
indexIntercept = match(lp.endo, names.p) - 1,
indexLP = match(lp(x, type = "name"), names.data) - 1,
indexManifest = match(lava::manifest(x), names.data) - 1,
indiv = indiv
)
colnames(score) <- coef(x)
}else{
#### compute
dataLP <- calcLP(data, p = p,
lp.x = lp.x, lp.link = lp.link, lp.name = lp(x, type = "name"))
score <- mets::scoreMVN(y = dataLP[,lava::manifest(x),drop = FALSE],
mu = M$xi%x%rep(1,NROW(dataLP)),
S= M$C,
dmu = D$dxi,
dS = D$dS)
colnames(score) <- coef(x)#c(name.intercept,name.regression,name.covariance)
## apply chain rule
n.lp <- length(x$lp)
for(iterLP in 1:n.lp){
name.intercept <- lp.endo[iterLP]
## extract data
X <- data[,lp.x[[iterLP]],drop = FALSE]
dlp <- X # - what about dB/db in presence of constrains
score[,lp.link[[iterLP]]] <- apply(dlp, 2, function(j){j*score[,name.intercept]})
}
if(indiv == FALSE){
score <- rbind(apply(score,2,sum))
}
}
return(score)
}
#' @rdname momentLVMr
#' @export
gaussianLP_hessian.lvm <- function(x, p, n, type,...) {
dots <- list(...);
if(type == "E"){
S <- -gaussianLP_score.lvm(x,p=p,data=dots$data,indiv = TRUE)
I <- t(S)%*%S
attributes(I)$grad <- colSums(S)
return(I)
}else if(type=="num"){
myg <- function(p1){ -gaussianLP_score.lvm(x,p=p1,n=n,data=dots$data,indiv=FALSE) }
I <- numDeriv::jacobian(myg,p)
I <- (I+t(I))/2
return( I )
}else if(type == "information"){ ## true part
browser()
S <- -gaussianLP_score.lvm(x,p=p,data=dots$data,indiv = FALSE)
lp.name <- lp(x, type = "name")
lp.link <- lp(x, type = "link", format = "list")
lp.x <- lp(x, type = "x", format = "list")
dataLP <- calcLP(dots$data, p = p,
lp.x = lp.x, lp.link = lp.link, lp.name = lp.name)
## direct
I <- information(x=x, p=p, n=n, data = dataLP)
return(I)
}
}
#' @rdname momentLVMr
#' @export
gaussian1LP_method.lvm <- gaussianLP_method.lvm
#' @rdname momentLVMr
#' @export
gaussian1LP_logLik.lvm <- gaussianLP_logLik.lvm
#' @rdname momentLVMr
#' @export
gaussian1LP_objective.lvm <- gaussianLP_objective.lvm
#' @rdname momentLVMr
#' @export
gaussian1LP_score.lvm <- gaussianLP_score.lvm
#' @rdname momentLVMr
#' @export
gaussian1LP_gradient.lvm <- gaussianLP_gradient.lvm
#' @rdname momentLVMr
#' @export
gaussian1LP_hessian.lvm <- function(x, type, ...){
gaussianLP_hessian.lvm(x, type = "num", ...)
}
#' @rdname momentLVMr
#' @export
gaussian2LP_method.lvm <- gaussianLP_method.lvm
#' @rdname momentLVMr
#' @export
gaussian2LP_logLik.lvm <- gaussianLP_logLik.lvm
#' @rdname momentLVMr
#' @export
gaussian2LP_objective.lvm <- gaussianLP_objective.lvm
#' @rdname momentLVMr
#' @export
gaussian2LP_score.lvm <- gaussianLP_score.lvm
#' @rdname momentLVMr
#' @export
gaussian2LP_gradient.lvm <- gaussianLP_gradient.lvm
#' @rdname momentLVMr
#' @export
gaussian2LP_hessian.lvm <- function(x, type, ...){
gaussianLP_hessian.lvm(x, type = "E", ...)
}
#' @title Compute the linear predictor
#' @description Compute the value of the linear predictors of a LVM and store it into the dataset
#'
#' @inheritParams momentLVMr
#' @param lp.x the name of the exogeneous variables involved in the LP
#' @param lp.link list containing for each LP the name of the links
#' @param lp.name list containing for each LP its name
calcLP <- function(data, p,
lp.x, lp.link, lp.name){
n.lp <- length(lp.name)
for(iterLP in 1:n.lp){
## extract coefficients according to constrains
b <- p[lp.link[[iterLP]]] # x$lp$y1$con
## compute linear predictor for the reduce model
X <- data[,lp.x[[iterLP]],drop = FALSE]
data[,lp.name[iterLP]] <- as.matrix(X) %*% b
}
return(data)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.