Nothing
R/gh_LL.R
In REMixed: Regularized Estimation in Mixed Effect Model
Defines functions
decomp
lambda.max
lambda.max.ind
amgauss.hermite
agauss.hermite
gh.LL.ind
gh.LL
Documented in
gh.LL
#' Adaptive Gauss-Hermite approximation of log-likelihood derivatives
#'
#' @description
#' Computes Adaptive Gauss-Hermite approximation of the log-likelihood and its derivatives in NLMEM with latent observation processes, see \code{\link{REMixed-package}} for details on the model.
#'
#' @details
#' Based on notation introduced \code{\link{REMixed-package}}. The log-likelihood of the model \eqn{LL(\theta,\alpha_1)} for a set of population parameters \eqn{\theta} and regulatization parameters \eqn{\alpha_1} is estimated using Adaptative Gausse-Hermite quadrature, using conditional distribution estimation to locate the mass of the integrand. If the project has been initialized as a Monolix project, the user can use \code{\link{readMLX}} function to retrieve all the project information needed here.
#'
#' @param dynFUN function computing the dynamics of interest for a set of parameters. This function need to contain every sub-function that it may needs (as it is called in a foreach loop). The output of this function need to return a data.frame with \code{time} : as first columns and named dynamics in other columns. It must take in input : \itemize{\item\code{y} : a named vector with the initial condition. The names are the dynamics names.
#' \item\code{parms} : a named vector of parameter.
#' \item\code{time} : vector a timepoint.}
#'
#' See \code{\link{dynFUN_demo}}, \code{\link{model.clairon}}, \code{\link{model.pasin}} or \code{\link{model.pk}} for examples.
#' @param y initial condition of the mechanism model, conform to what is asked in \code{dynFUN}.
#' @param mu list of individuals random effects estimation (vector of r.e. need to be named by the parameter names), use to locate the density mass; (optional, see description).
#' @param Omega list of individuals estimated standard deviation diagonal matrix (matrix need to have rows and columns named by the parameter names), use to locate the density mass; (optional, see description).
#' @param theta list of model parameters containing (see details) \itemize{\item\code{phi_pop} : named vector with the population parameters with no r.e. \eqn{(\phi_{l\ pop})_{l\leq L}} (NULL if none) ;
#' \item\code{psi_pop} : named vector with the population parameters with r.e. \eqn{(\psi_{l\ pop})_{l\leq m}} ;
#' \item\code{gamma} : named list (for each parameters) of named vector (for each covariates) of covariate effects from parameters with no r.e. ;
#' \item\code{beta} : named list (for each parameters) of named vector (for each covariates) of covariate effects from parameters with r.e..
#' \item\code{alpha0} : named vector of \eqn{(\alpha_{0k})_{k\leq K}} parameters (names are identifier of the observation model, such as in a Monolix project);
#' \item\code{omega} : named vector of estimated r.e. standard deviation;}
#' (optional, see description).
#' @param alpha1 named vector of regulatization parameters \eqn{(\alpha_{1k})_{k\leq K}}, with identifier of observation model as names, (optional, see description).
#' @param covariates matrix of individual covariates (size N x n). Individuals must be sorted in the same order than in \code{mu} and \code{Omega}, (optional, see description).
#' @param ParModel.transfo named list of transformation functions \eqn{(h_l)_{l\leq m}} and \eqn{(s_k)_{k\leq K}} for the individual parameter model (names must be consistent with \code{phi_pop} and \code{psi_pop}, missing entries are set by default to the identity function ; optional, see description).
#' @param ParModel.transfo.inv Named list of inverse transformation functions for the individual parameter model (names must be consistent with \code{phi_pop} and \code{psi_pop} ; optional, see description).
#' @param Sobs list of individuals trajectories for the direct observation models \eqn{(Y_{pi})_{p \leq P,i\leq N}}. Each element \eqn{i\leq N} of the list, is a list of \eqn{p\leq P} data.frame with time \eqn{(t_{pij})_{j\leq n_{ip}}} and observations \eqn{(Y_{pij})_{j\leq n_{ip}}}. Each data.frame is named with the observation model identifiers.
#' @param Robs list of individuals trajectories for the latent observation models \eqn{(Z_{ki})_{k \leq K,i\leq N}}. Each element \eqn{i\leq N} of the list, is a list of \eqn{k\leq K} data.frame with time \eqn{(t_{kij})_{j\leq n_{ik}}} and observations \eqn{(Z_{kij})_{j\leq n_{ik}}}. Each data.frame is named with the observation model identifiers.
#' @param Serr named vector of the estimated error mocel constants \eqn{(\varsigma_p)_{p\leq P}} with observation model identifiers as names.
#' @param Rerr named vector of the estimated error mocel constants \eqn{(\sigma_k)_{k\leq K}} with observation model identifiers as names.
#' @param ObsModel.transfo list containing two lists of transformations and two vectors linking each transformations to their observation model name in the Monolix project. The list should include identity transformations and be named \code{S} and \code{R}. The two vectors should be named \code{linkS} and \code{linkR}.
#'
#' Both \code{S} (for the direct observation models) and \code{linkS}, as well as \code{R} (for latent process models) and \code{linkR}, must have the same length.
#'
#' \itemize{
#' \item\code{S}: a list of transformations for the direct observation models. Each transformation corresponds to a variable \eqn{Y_p=h_p(S_p)}, where the name indicates which dynamic is observed (from \code{dynFUN}); \item\code{linkS} : a vector specifying the observation model names (that is used in the monolix project, \code{alpha1}, etc.) for each transformation, in the same order as in \code{S};
#'
#' \item\code{R}: similarly, a list of transformations for the latent process models. Although currently there is only one latent dynamic, each \eqn{s_k, k\leq K} transformation corresponds to the same dynamic but may vary for each \eqn{Y_k} observed. The names should match the output from \code{dynFUN}; \item \code{linkR} : a vector specifying the observation model names for each transformation, in the same order as in \code{R}.
#' }
#' @param data output from \code{\link{readMLX}} containing parameters "\code{mu}", "\code{Omega}", "\code{theta}", "\code{alpha1}", "\code{covariates}", "\code{ParModel.transfo}", "\code{ParModel.transfo.inv}", "\code{Sobs}", "\code{Robs}", "\code{Serr}", "\code{Rerr}", "\code{ObsModel.transfo}" extract from a monolix project.
#' @param n number of points per dimension to use for the Gauss-Hermite quadrature rule.
#' @param prune integer between 0 and 1, percentage of pruning for the Gauss-Hermite quadrature rule (default NULL).
#' @param parallel logical, if computation should be done in parallel.
#' @param ncores number of cores to use for parallelization, default will detect the number of cores available.
#' @param onlyLL logical, if only the log-likelihood should be computed (and not \eqn{\partial_{\alpha_1} LL} or \eqn{\partial_{\alpha_1}^2 LL}).
#' @param verbose logical, if progress bar should be printed through the computation.
#'
#' @return A list with the approximation by Gauss-Hermite quadrature of the likelihood \code{L}, the log-likelihood \code{LL}, the gradient of the log-likelihood \code{dLL}, and the Hessian of the log-likelihood \code{ddLL} at the point \eqn{\theta, \alpha} provided.
#' @export
#'
#' @examples
#'\dontrun{
#' project <- getMLXdir()
#'
#'
#' ObsModel.transfo = list(S=list(AB=log10),
#' linkS="yAB",
#' R=rep(list(S=function(x){x}),5),
#' linkR = paste0("yG",1:5))
#'
#' alpha=list(alpha0=NULL,
#' alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))
#'
#' data <- readMLX(project,ObsModel.transfo,alpha)
#'
#' LL <- gh.LL(dynFUN = dynFUN_demo,
#' y = c(S=5,AB=1000),
#' ObsModel.transfo=ObsModel.transfo,
#' data = data)
#'
#' print(LL)
#' }
gh.LL <- function(
dynFUN,
y,
mu=NULL,
Omega=NULL,
theta=NULL,
alpha1=NULL,
covariates=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs=NULL,
Robs=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
n=NULL,
prune=NULL,
parallel=TRUE,
ncores=NULL,
onlyLL=FALSE,
verbose=TRUE){
if(is.null(data)){
test <- sapply(c("mu","Omega","theta","alpha1","covariates","ParModel.transfo","ParModel.transfo.inv","Sobs","Robs","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu) || is.null(Omega)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
for(d in 1:length(data)){
test <- eval(parse(text=(paste0("is.null(",names(data)[d],")"))))
if(test){
eval(parse(text=paste0(names(data[d]),"<- data[[d]]")))
}
}
}
if(is.null(n)){
if(length(theta$psi_pop)==1){
n <- 100
}else if(length(theta$psi_pop)==2){
n <- 10
}else{
n <- 7
}
}
if(parallel){
if(!is.null(ncores)){
cluster <- snow::makeCluster(ncores)
}else{
cluster <- snow::makeCluster(parallel::detectCores())
}
doSNOW::registerDoSNOW(cluster)
}
N=length(mu)
i = 1
ntasks <- N
if(verbose){
pb <- txtProgressBar(max = ntasks, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)
}else{
opts <- NULL
}
res = foreach::foreach(i = 1:N,.packages = "REMixed",.options.snow=opts)%dopar%{
if(0 %in% diag(Omega[[i]])){
diag(Omega[[i]])[diag(Omega[[i]])==0] <- 10**(-5)
}
LLi = gh.LL.ind(mu_i = mu[[i]],
Omega_i = Omega[[i]],
theta = theta,
alpha1 = alpha1,
dynFUN = dynFUN,
y = y,
covariates_i = covariates[i,,drop=F],
ParModel.transfo = ParModel.transfo,
ParModel.transfo.inv = ParModel.transfo.inv,
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==i,]}),
Robs_i = lapply(Robs,FUN=function(R){R[R$id==i,]}),
Serr = Serr,
Rerr = Rerr,
ObsModel.transfo = ObsModel.transfo,
ind = i,
n = n,
prune = prune,
onlyLL = onlyLL)
return(LLi)
}
L = prod(sapply(res,FUN=function(ri){ri$Li}))
LL = sum(sapply(res,FUN=function(ri){ri$LLi}))
if(verbose)
close(pb)
if(parallel){
snow::stopCluster(cluster)
}
if(!onlyLL){
if(length(alpha1)!=1){
dLL = rowSums(sapply(res,FUN=function(ri){ri$dLLi}))
ddLL = matrix(rowSums(sapply(res,FUN=function(ri){ri$ddLLi})),ncol=length(dLL))
}else{
dLL = sum(sapply(res,FUN=function(ri){ri$dLLi}))
ddLL = matrix(sum(sapply(res,FUN=function(ri){ri$ddLLi})),ncol=length(dLL))
}
return(list(L=L,LL=LL,dLL=dLL,ddLL=ddLL))
}else{
return(LL)
}
}
gh.LL.ind <- function(
dynFUN,
y,
mu_i=NULL,
Omega_i=NULL,
theta=NULL,
alpha1=NULL,
covariates_i=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs_i=NULL,
Robs_i=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
ind = NULL,
n = NULL,
prune=NULL,
onlyLL=FALSE){
mu <- Omega <- Sobs <- Robs <- covariates <- NULL
if(is.null(data)){
test <- sapply(c("mu_i","Omega_i","theta","alpha1","covariates_i","ParModel.transfo","ParModel.transfo.inv","Sobs_i","Robs_i","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu_i) || is.null(Omega_i)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
argsNONind = c("theta","alpha1","ParModel.transfo","ParModel.transfo.inv","Serr","Rerr","ObsModel.transfo")
for(d in 1:length(argsNONind)){
test <- eval(parse(text=paste0("is.null(",argsNONind[d],")")))
if(test){
eval(parse(text=paste0(argsNONind[d],"<- data[[argsNONind[d]]]")))
}
}
argsInd = c("mu","Omega","Robs","Sobs","covariates")
for(d in 1:length(argsInd)){
test <- eval(parse(text=paste0("is.null(",paste0(argsInd[d],"_i"),")")))
if(test){
eval(parse(text=paste0(argsInd[d],"<- data[[argsInd[d]]]")))
}
}
if((is.null(mu_i) || is.null(Omega_i) || is.null(Robs_i) || is.null(Sobs_i) || is.null(covariates_i)) & is.null(ind)){
stop("Please provide individual information (mu_i,Omega_i,Sobs_i,Robs_i,covariates_i), or the individual id ind.")
}
if(is.null(mu_i)){
mu_i = mu[[ind]]
}
if(is.null(Omega_i)){
Omega_i = Omega[[ind]]
}
if(is.null(Sobs_i)){
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==ind,]})
}
if(is.null(Robs_i)){
Robs_i = lapply(Robs,FUN=function(R){R[R$id==ind,]})
}
if(is.null(covariates_i)){
covariates_i = covariates[ind,,drop=F]
}
}
if(is.null(n)){
if(length(theta$psi_pop)==1){
n <- 100
}else if(length(theta$psi_pop)==2){
n <- 10
}else{
n <- 7
}
}
dm = length(mu_i)
if(dm!=1){
mh.parm <- amgauss.hermite(n,mu=mu_i,Omega=Omega_i,prune=prune)
detsq.omega = prod(theta$omega)
root.omega = diag(1/theta$omega**2)
}else{
mh.parm <- agauss.hermite(n,mu = mu_i,sd=sqrt(as.numeric(Omega_i)))
detsq.omega = as.numeric(theta$omega)
root.omega = 1/as.numeric(theta$omega)**2
}
nd = length(mh.parm$Weights)
R.sz = length(Rerr)
S.sz = length(Serr)
all.tobs = sort(union(unlist(lapply(Robs_i,FUN=function(x){x$time})),unlist(lapply(Sobs_i,FUN=function(x){x$time}))))
# Need to compute, for each eta_i x individual i, the dynamics of the model
# eta_i = split(mh.parm$Points,1:nd)[[1]]
dyn <- setNames(lapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
PSI_i = indParm(theta[c("phi_pop","psi_pop","gamma","beta")],covariates_i,setNames(eta_i,colnames(Omega_i)),ParModel.transfo,ParModel.transfo.inv)
dyn_eta_i <- dynFUN(all.tobs,y,unlist(unname(PSI_i)))
return(dyn_eta_i)
}),paste0("eta_",1:nd))
# Compute marginal latent density for each individual and value of RE
R.margDensity <- setNames(lapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
# .Machine$double.xmin
decomp_marg = lapply(1:R.sz,FUN=function(k){
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
aux <- 1/(sig*sqrt(2*pi))*exp(-1/2*((Zki-a0-a1*Rki)/sig)**2)
if(any(aux==0)){
pw <- Rmpfr::mpfr(-1/2*((Zki-a0-a1*Rki)/sig)**2,precBits = 32)
aux = 1/(sig*sqrt(2*pi))*exp(pw)
decomp_aux = setNames(sapply(decomp(prod(aux)),as.numeric),c("exponent","mantissa"))
}else{
decomp_aux <- lapply(aux,decomp)
decomp_intermediaire = decomp(prod(sapply(decomp_aux,FUN=function(decomp){decomp["mantissa"]})))
mantissa_aux = decomp_intermediaire["mantissa"]
exponent_aux = sum(sapply(decomp_aux,FUN=function(decomp){decomp["exponent"]})) + decomp_intermediaire["exponent"]
decomp_aux <- c(exponent_aux,mantissa_aux)
}
return(decomp_aux)
})
if(any(sapply(decomp_marg,function(x){x["mantissa"]})==0)){
stop("[Error] in log-likelihood computation, latent part of the marginal likelihood is equal to zero.")
}
# decomp_marg = lapply(little_marg,decomp)
decomp_intermediaire = decomp(prod(sapply(decomp_marg,FUN=function(decomp){decomp["mantissa"]})))
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = sum(sapply(decomp_marg,FUN=function(decomp){decomp["exponent"]})) + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
}),paste0("eta_",1:nd))
# Compute marginal density for each individual and value of RE
if(S.sz!=0){
S.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
eta_i = split(mh.parm$Points,1:nd)[[ei]]
dyn.ei = dyn[[paste0("eta_",ei)]]
res = prod(sapply(1:S.sz,FUN=function(p){
Yp = ObsModel.transfo$linkS[p]
sig = Serr[[Yp]]
tpi = Sobs_i[[Yp]]$time
Ypi = Sobs_i[[Yp]][,Yp]
trs = ObsModel.transfo$S[which(ObsModel.transfo$linkS==Yp)]
Spi = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tpi,names(trs)])
prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Ypi-Spi)/sig)**2))
}))*(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}else{
S.margDensity = setNames(sapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}
# Compute individual log-Likelihood
Li.aux = lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(Li.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(Li.aux,function(li){
return(li[["mantissa"]]*10**(li[["exponent"]] - max_exponent))
})))
Li = c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
# sum(mh.parm$Weights*R.margDensity*S.margDensity)
invLi.aux = decomp(1/Li[["mantissa"]])
invLi = c(exponent=-Li[["exponent"]]+invLi.aux[["exponent"]],mantissa=invLi.aux[["mantissa"]])
# ln(p*10**o)=ln(p)+ln(10**o)=ln(p)+o*ln(10)
LLi = log(Li[["mantissa"]])+Li[["exponent"]]*log(10)
if(!onlyLL){
# Compute gradient of individual log-Likelihood
# dim.alpha = length(alpha1) # = K = length(Robs_i )
dfact = lapply(1:R.sz,FUN=function(k){
f = setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))
}),paste0("eta_",1:nd))
})
dLLi = sapply(dfact,FUN=function(f){
res.aux <- lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*f[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(res.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(res.aux,function(ri){
return(ri[["mantissa"]]*10**(ri[["exponent"]] - max_exponent))
})))
aux =c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
return(invLi[["mantissa"]]*aux[["mantissa"]]*10**(invLi[["exponent"]]+aux[["exponent"]]))
# 1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
})
# Compute hessienne of individual log-Likelihood
ddLLi = matrix(0,ncol=R.sz,nrow=R.sz)
# diagonal term
ddfact = lapply(1:R.sz,FUN=function(k){
setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(
(1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))**2 - 1/(sig**2)*(sum(Rki**2))
)
}),paste0("eta_",1:nd))
})
diag(ddLLi) = sapply(ddfact,FUN=function(f){
res.aux <- lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*f[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(res.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(res.aux,function(ri){
return(ri[["mantissa"]]*10**(ri[["exponent"]] - max_exponent))
})))
aux =c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
return(invLi[["mantissa"]]*aux[["mantissa"]]*10**(invLi[["exponent"]]+aux[["exponent"]]))
# return(1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
})-(dLLi)**2
# non diagonal term
for(k1 in 1:R.sz){
for(k2 in 1:R.sz){
if(k1<k2){
ddfact = setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
prod(sapply(c(k1,k2),FUN=function(k){
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(
1/sig**2*(sum(Rki*(Zki-a0-a1*Rki)))
)
}))
}),paste0("eta_",1:nd))
# ddLLi[k1,k2] <- ddLLi[k2,k1] <- 1/Li*sum(mh.parm$Weights*ddfact*R.margDensity*S.margDensity)-prod(dLLi[c(k1,k2)])
res.aux <- lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*ddfact[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(res.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(res.aux,function(ri){
return(ri[["mantissa"]]*10**(ri[["exponent"]] - max_exponent))
})))
aux =c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
ddLLi[k1,k2] <- ddLLi[k2,k1] <- (invLi[["mantissa"]]*aux[["mantissa"]]*10**(invLi[["exponent"]]+aux[["exponent"]])) - prod(dLLi[c(k1,k2)])
}
}
}
}else{
dLLi <- ddLLi <- NULL
}
return(list(Li=Li[["mantissa"]]*10**(Li[["exponent"]]),LLi=LLi,dLLi=dLLi,ddLLi=ddLLi))
}
# One dimensional Adaptative Hermite Gauss ------------------------------------------
agauss.hermite <- function(n,mu=0,sd=1){
gh <- fastGHQuad::gaussHermiteData(n)
names(gh) <- c("Points","Weights")
ww <- gh$Weights * sqrt(2) * sd * exp(gh$Points**2)
xx <- mu + sd * sqrt(2) * gh$Points
return(data.frame(Points=xx,Weights=ww))
}
# Multi-dimensional Adaptative Hermite Gauss -----------------------------------------
amgauss.hermite <- function(n,mu=rep(0,ncol(Omega)),Omega=diag(rep(1,length(mu))),prune=NULL){
dm = length(mu)
gh <- as.data.frame(fastGHQuad::gaussHermiteData(n))
idx <- as.matrix(expand.grid(rep(list(1:n),dm)))
pts <- matrix(gh[idx,1],nrow(idx),dm)
wts <- apply(matrix(gh[idx,2],nrow(idx),dm), 1, prod)
# plot(pts, cex=-5/log(wts), pch=19)
if(!is.null(prune)) {
qwt <- quantile(wts, probs=prune)
pts <- pts[wts > qwt,]
wts <- wts[wts > qwt]
}
eig <- eigen(Omega)
S = eig$vectors
L = diag(sqrt(eig$values))
rot <- S %*% L
ppts <-t( sqrt(2) * rot %*% t(pts) + mu)
wwts <- wts * sqrt(2)**dm * abs(det(rot)) * exp(apply(pts,1,FUN=function(x){t(x)%*%x}))
return(list(Points=ppts, Weights=wwts))
}
# individual contribution to log-lik gradient in alpha1=0 ----------------------------
lambda.max.ind <- function(
dynFUN,
y,
mu_i=NULL,
Omega_i=NULL,
theta=NULL,
alpha1=NULL,
covariates_i=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs_i=NULL,
Robs_i=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
ind = NULL,
n = NULL,
prune=NULL){
mu <- Omega <- Sobs <- Robs <- covariates <- NULL
if(is.null(data)){
test <- sapply(c("mu_i","Omega_i","theta","alpha1","covariates_i","ParModel.transfo","ParModel.transfo.inv","Sobs_i","Robs_i","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu_i) || is.null(Omega_i)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
argsNONind = c("theta","alpha1","ParModel.transfo","ParModel.transfo.inv","Serr","Rerr","ObsModel.transfo")
for(d in 1:length(argsNONind)){
test <- eval(parse(text=paste0("is.null(",argsNONind[d],")")))
if(test){
eval(parse(text=paste0(argsNONind[d],"<- data[[argsNONind[d]]]")))
}
}
argsInd = c("mu","Omega","Robs","Sobs","covariates")
for(d in 1:length(argsInd)){
test <- eval(parse(text=paste0("is.null(",paste0(argsInd[d],"_i"),")")))
if(test){
eval(parse(text=paste0(argsInd[d],"<- data[[argsInd[d]]]")))
}
}
if((is.null(mu_i) || is.null(Omega_i) || is.null(Robs_i) || is.null(Sobs_i) || is.null(covariates_i)) & is.null(ind)){
stop("Please provide individual information (mu_i,Omega_i,Sobs_i,Robs_i,covariates_i), or the individual id ind.")
}
if(is.null(mu_i)){
mu_i = mu[[ind]]
}
if(is.null(Omega_i)){
Omega_i = Omega[[ind]]
}
if(is.null(Sobs_i)){
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==ind,]})
}
if(is.null(Robs_i)){
Robs_i = lapply(Robs,FUN=function(R){R[R$id==ind,]})
}
if(is.null(covariates_i)){
covariates_i = covariates[ind,,drop=F]
}
}
if(is.null(n)){
n <- floor(100**(1/length(theta$psi_pop)))
}
dm = length(mu_i)
if(dm!=1){
mh.parm <- amgauss.hermite(n,mu=mu_i,Omega=Omega_i,prune=prune)
detsq.omega = prod(theta$omega)
root.omega = diag(1/theta$omega**2)
}else{
mh.parm <- agauss.hermite(n,mu = mu_i,sd=sqrt(as.numeric(Omega_i)))
detsq.omega = as.numeric(theta$omega)
root.omega = 1/as.numeric(theta$omega)**2
}
nd = length(mh.parm$Weights)
R.sz = length(Rerr)
S.sz = length(Serr)
all.tobs = sort(union(unlist(lapply(Robs_i,FUN=function(x){x$time})),unlist(lapply(Sobs_i,FUN=function(x){x$time}))))
dyn <- setNames(lapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
PSI_i = indParm(theta[c("phi_pop","psi_pop","gamma","beta")],covariates_i,setNames(eta_i,colnames(Omega_i)),ParModel.transfo,ParModel.transfo.inv)
dyn_eta_i <- dynFUN(all.tobs,y,unlist(unname(PSI_i)))
return(dyn_eta_i)
}),paste0("eta_",1:nd))
# Compute marginal latent density for each individual and value of RE
R.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
res = prod(sapply(1:R.sz,FUN=function(k){
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Zki-a0-a1*Rki)/sig)**2))
}))
}),paste0("eta_",1:nd))
# A(Y|theta0)
if(S.sz!=0){
S.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
eta_i = split(mh.parm$Points,1:nd)[[ei]]
dyn.ei = dyn[[paste0("eta_",ei)]]
res = prod(sapply(1:S.sz,FUN=function(p){
Yp = ObsModel.transfo$linkS[p]
sig = Serr[[Yp]]
tpi = Sobs_i[[Yp]]$time
Ypi = Sobs_i[[Yp]][,Yp]
trs = ObsModel.transfo$S[which(ObsModel.transfo$linkS==Yp)]
Spi = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tpi,names(trs)])
prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Ypi-Spi)/sig)**2))
}))*(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}else{
S.margDensity = setNames(sapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}
# s'k(Ri(krij))A(Y|theta0)
S2.margDensity <- lapply(1:R.sz,FUN=function(k){
setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(1/sig**2*sum((Zki-a0)*Rki)*S.margDensity[[ei]])
}),paste0("eta_",1:nd))})
indi.contrib <- sapply(1:R.sz,FUN=function(k){
sum(mh.parm$Weights*S2.margDensity[[k]])
})/sum(mh.parm$Weights*S.margDensity)
return(indi.contrib)
}
# log-lik gradient in alpha1=0 for maximum penalization parameter to reach ----------
lambda.max <- function(
dynFUN,
y,
mu=NULL,
Omega=NULL,
theta=NULL,
alpha1=NULL,
covariates=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs=NULL,
Robs=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
n = NULL,
prune=NULL,
parallel = TRUE,
ncores=NULL,
onlyLL=FALSE,
verbose=TRUE){
if(is.null(data)){
test <- sapply(c("mu","Omega","theta","alpha1","covariates","ParModel.transfo","ParModel.transfo.inv","Sobs","Robs","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu) || is.null(Omega)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
for(d in 1:length(data)){
test <- eval(parse(text=(paste0("is.null(",names(data)[d],")"))))
if(test){
eval(parse(text=paste0(names(data[d]),"<- data[[d]]")))
}
}
}
if(is.null(n)){
n <- floor(100**(1/length(theta$psi_pop)))
}
if(parallel){
if(!is.null(ncores)){
cluster <- parallel::makeCluster(ncores)
}else{
cluster <- parallel::makeCluster(parallel::detectCores())
}
doSNOW::registerDoSNOW(cluster)
}
N=length(mu)
i = 1
if(verbose){
ntasks <- N
pb <- txtProgressBar(max = ntasks, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)
}else{
opts <- NULL
}
res = foreach::foreach(i = 1:N,.packages = "REMixed",.export = c("lambda.max.ind","amgauss.hermite"),.options.snow=opts)%dopar%{
if(0 %in% diag(Omega[[i]])){
diag(Omega[[i]])[diag(Omega[[i]])==0] <- 10**(-5)
}
lambda.max.ind(mu_i = mu[[i]],
Omega_i = Omega[[i]],
theta = theta,
alpha1 = alpha1,
dynFUN = dynFUN,
y = y,
covariates_i = covariates[i,,drop=F],
ParModel.transfo = ParModel.transfo,
ParModel.transfo.inv = ParModel.transfo.inv,
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==i,]}),
Robs_i = lapply(Robs,FUN=function(R){R[R$id==i,]}),
Serr = Serr,
Rerr = Rerr,
ObsModel.transfo = ObsModel.transfo,
n = n,
prune = prune)
}
if(verbose)
close(pb)
if(parallel)
snow::stopCluster(cluster)
return(max(Reduce("+",res)))
}
decomp <- function(x){
if(x==0){
return(c(exponent=0,mantissa=0))
}
exponent_x <- floor(log10(abs(x))) # Exponent
exponent_lim <- abs(floor(log10(.Machine$double.xmin)))
if(abs(exponent_x)>exponent_lim){
mantissa_x <- x*10**(-sign(exponent_x)*(exponent_lim-1))
mantissa_x <- mantissa_x*10**(-(exponent_x+(exponent_lim-1)))
}else{
mantissa_x <- x*10**(-exponent_x)
}
return(c(exponent=exponent_x,mantissa=mantissa_x))
}
#
# gh.LL.ind_OLD <- function(
# dynFUN,
# y,
# mu_i=NULL,
# Omega_i=NULL,
# theta=NULL,
# alpha1=NULL,
# covariates_i=NULL,
# ParModel.transfo=NULL,
# ParModel.transfo.inv=NULL,
# Sobs_i=NULL,
# Robs_i=NULL,
# Serr=NULL,
# Rerr=NULL,
# ObsModel.transfo=NULL,
# data=NULL,
# ind = NULL,
# n = NULL,
# prune=NULL,
# onlyLL=FALSE){
#
# mu <- Omega <- Sobs <- Robs <- covariates <- NULL
#
# if(is.null(data)){
# test <- sapply(c("mu_i","Omega_i","theta","alpha1","covariates_i","ParModel.transfo","ParModel.transfo.inv","Sobs_i","Robs_i","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
# if(any(test))
# stop("Please provide all necessary arguments.")
# }else{
# if((is.null(mu_i) || is.null(Omega_i)) & !all(c("mu","Omega") %in% names(data))){
# stop("Please provide mu and Omega if these are missing from data.")
# }
# argsNONind = c("theta","alpha1","ParModel.transfo","ParModel.transfo.inv","Serr","Rerr","ObsModel.transfo")
# for(d in 1:length(argsNONind)){
# test <- eval(parse(text=paste0("is.null(",argsNONind[d],")")))
# if(test){
# eval(parse(text=paste0(argsNONind[d],"<- data[[argsNONind[d]]]")))
# }
# }
# argsInd = c("mu","Omega","Robs","Sobs","covariates")
# for(d in 1:length(argsInd)){
# test <- eval(parse(text=paste0("is.null(",paste0(argsInd[d],"_i"),")")))
# if(test){
# eval(parse(text=paste0(argsInd[d],"<- data[[argsInd[d]]]")))
#
# }
# }
# if((is.null(mu_i) || is.null(Omega_i) || is.null(Robs_i) || is.null(Sobs_i) || is.null(covariates_i)) & is.null(ind)){
# stop("Please provide individual information (mu_i,Omega_i,Sobs_i,Robs_i,covariates_i), or the individual id ind.")
# }
# if(is.null(mu_i)){
# mu_i = mu[[ind]]
# }
# if(is.null(Omega_i)){
# Omega_i = Omega[[ind]]
# }
# if(is.null(Sobs_i)){
# Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==ind,]})
# }
# if(is.null(Robs_i)){
# Robs_i = lapply(Robs,FUN=function(R){R[R$id==ind,]})
# }
# if(is.null(covariates_i)){
# covariates_i = covariates[ind,,drop=F]
# }
# }
#
# if(is.null(n)){
# n <- floor(100**(1/length(theta$psi_pop)))
# }
#
# dm = length(mu_i)
#
# if(dm!=1){
# mh.parm <- amgauss.hermite(n,mu=mu_i,Omega=Omega_i,prune=prune)
# detsq.omega = prod(theta$omega)
# root.omega = diag(1/theta$omega**2)
# }else{
# mh.parm <- agauss.hermite(n,mu = mu_i,sd=sqrt(as.numeric(Omega_i)))
# detsq.omega = as.numeric(theta$omega)
# root.omega = 1/as.numeric(theta$omega)**2
# }
#
# nd = length(mh.parm$Weights)
# R.sz = length(Rerr)
# S.sz = length(Serr)
# all.tobs = sort(union(unlist(lapply(Robs_i,FUN=function(x){x$time})),unlist(lapply(Sobs_i,FUN=function(x){x$time}))))
#
# # Need to compute, for each eta_i x individual i, the dynamics of the model
# # eta_i = split(mh.parm$Points,1:nd)[[1]]
# dyn <- setNames(lapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
# PSI_i = indParm(theta[c("phi_pop","psi_pop","gamma","beta")],covariates_i,setNames(eta_i,colnames(Omega_i)),ParModel.transfo,ParModel.transfo.inv)
# dyn_eta_i <- dynFUN(all.tobs,y,unlist(unname(PSI_i)))
#
# return(dyn_eta_i)
# }),paste0("eta_",1:nd))
#
# # Compute marginal latent density for each individual and value of RE
# R.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
#
#
# res = prod(sapply(1:R.sz,FUN=function(k){
# yGk = ObsModel.transfo$linkR[k]
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Zki-a0-a1*Rki)/sig)**2))
# }))
# }),paste0("eta_",1:nd))
#
# # Compute marginal density for each individual and value of RE
# if(S.sz!=0){
# S.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
# eta_i = split(mh.parm$Points,1:nd)[[ei]]
# dyn.ei = dyn[[paste0("eta_",ei)]]
#
# res = prod(sapply(1:S.sz,FUN=function(p){
# Yp = ObsModel.transfo$linkS[p]
# sig = Serr[[Yp]]
# tpi = Sobs_i[[Yp]]$time
# Ypi = Sobs_i[[Yp]][,Yp]
#
# trs = ObsModel.transfo$S[which(ObsModel.transfo$linkS==Yp)]
# Spi = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tpi,names(trs)])
#
# prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Ypi-Spi)/sig)**2))
# }))*(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
# }),paste0("eta_",1:nd))
# }else{
# S.margDensity = setNames(sapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
# (1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
# }),paste0("eta_",1:nd))
# }
#
# # Compute individual log-Likelihood
# Li = sum(mh.parm$Weights*R.margDensity*S.margDensity)
# LLi = log(Li)
#
# if(!onlyLL){
# # Compute gradient of individual log-Likelihood
# # dim.alpha = length(alpha1) # = K = length(Robs_i )
# dfact = lapply(1:R.sz,FUN=function(k){
# setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
# yGk = ObsModel.transfo$linkR[k]
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# return(1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))
# }),paste0("eta_",1:nd))
# })
#
# dLLi = sapply(dfact,FUN=function(f){
# return(1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
# })
#
# # Compute hessienne of individual log-Likelihood
# ddLLi = matrix(0,ncol=R.sz,nrow=R.sz)
# # diagonal term
# ddfact = lapply(1:R.sz,FUN=function(k){
# setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
# yGk = ObsModel.transfo$linkR[k]
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# return(
# (1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))**2 - 1/(sig**2)*(sum(Rki**2))
# )
# }),paste0("eta_",1:nd))
# })
#
# diag(ddLLi) = sapply(ddfact,FUN=function(f){
# return(1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
# })-(dLLi)**2
#
# # non diagonal term
# for(k1 in 1:R.sz){
# for(k2 in 1:R.sz){
# if(k1<k2){
# ddfact = setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
# prod(sapply(c(k1,k2),FUN=function(k){
# yGk = ObsModel.transfo$linkR[k]
#
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# return(
# 1/sig**2*(sum(Rki*(Zki-a0-a1*Rki)))
# )
# }))
# }),paste0("eta_",1:nd))
#
# ddLLi[k1,k2] <- ddLLi[k2,k1] <- 1/Li*sum(mh.parm$Weights*ddfact*R.margDensity*S.margDensity)-prod(dLLi[c(k1,k2)])
# }
# }
# }
# }else{
# dLLi <- ddLLi <- NULL
# }
# return(list(Li=Li,LLi=LLi,dLLi=dLLi,ddLLi=ddLLi))
# }
Try the REMixed package in your browser
Any scripts or data that you put into this service are public.
REMixed documentation built on June 28, 2025, 1:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions decomp lambda.max lambda.max.ind amgauss.hermite agauss.hermite gh.LL.ind gh.LL
Documented in gh.LL
#' Adaptive Gauss-Hermite approximation of log-likelihood derivatives
#'
#' @description
#' Computes Adaptive Gauss-Hermite approximation of the log-likelihood and its derivatives in NLMEM with latent observation processes, see \code{\link{REMixed-package}} for details on the model.
#'
#' @details
#' Based on notation introduced \code{\link{REMixed-package}}. The log-likelihood of the model \eqn{LL(\theta,\alpha_1)} for a set of population parameters \eqn{\theta} and regulatization parameters \eqn{\alpha_1} is estimated using Adaptative Gausse-Hermite quadrature, using conditional distribution estimation to locate the mass of the integrand. If the project has been initialized as a Monolix project, the user can use \code{\link{readMLX}} function to retrieve all the project information needed here.
#'
#' @param dynFUN function computing the dynamics of interest for a set of parameters. This function need to contain every sub-function that it may needs (as it is called in a foreach loop). The output of this function need to return a data.frame with \code{time} : as first columns and named dynamics in other columns. It must take in input : \itemize{\item\code{y} : a named vector with the initial condition. The names are the dynamics names.
#' \item\code{parms} : a named vector of parameter.
#' \item\code{time} : vector a timepoint.}
#'
#' See \code{\link{dynFUN_demo}}, \code{\link{model.clairon}}, \code{\link{model.pasin}} or \code{\link{model.pk}} for examples.
#' @param y initial condition of the mechanism model, conform to what is asked in \code{dynFUN}.
#' @param mu list of individuals random effects estimation (vector of r.e. need to be named by the parameter names), use to locate the density mass; (optional, see description).
#' @param Omega list of individuals estimated standard deviation diagonal matrix (matrix need to have rows and columns named by the parameter names), use to locate the density mass; (optional, see description).
#' @param theta list of model parameters containing (see details) \itemize{\item\code{phi_pop} : named vector with the population parameters with no r.e. \eqn{(\phi_{l\ pop})_{l\leq L}} (NULL if none) ;
#' \item\code{psi_pop} : named vector with the population parameters with r.e. \eqn{(\psi_{l\ pop})_{l\leq m}} ;
#' \item\code{gamma} : named list (for each parameters) of named vector (for each covariates) of covariate effects from parameters with no r.e. ;
#' \item\code{beta} : named list (for each parameters) of named vector (for each covariates) of covariate effects from parameters with r.e..
#' \item\code{alpha0} : named vector of \eqn{(\alpha_{0k})_{k\leq K}} parameters (names are identifier of the observation model, such as in a Monolix project);
#' \item\code{omega} : named vector of estimated r.e. standard deviation;}
#' (optional, see description).
#' @param alpha1 named vector of regulatization parameters \eqn{(\alpha_{1k})_{k\leq K}}, with identifier of observation model as names, (optional, see description).
#' @param covariates matrix of individual covariates (size N x n). Individuals must be sorted in the same order than in \code{mu} and \code{Omega}, (optional, see description).
#' @param ParModel.transfo named list of transformation functions \eqn{(h_l)_{l\leq m}} and \eqn{(s_k)_{k\leq K}} for the individual parameter model (names must be consistent with \code{phi_pop} and \code{psi_pop}, missing entries are set by default to the identity function ; optional, see description).
#' @param ParModel.transfo.inv Named list of inverse transformation functions for the individual parameter model (names must be consistent with \code{phi_pop} and \code{psi_pop} ; optional, see description).
#' @param Sobs list of individuals trajectories for the direct observation models \eqn{(Y_{pi})_{p \leq P,i\leq N}}. Each element \eqn{i\leq N} of the list, is a list of \eqn{p\leq P} data.frame with time \eqn{(t_{pij})_{j\leq n_{ip}}} and observations \eqn{(Y_{pij})_{j\leq n_{ip}}}. Each data.frame is named with the observation model identifiers.
#' @param Robs list of individuals trajectories for the latent observation models \eqn{(Z_{ki})_{k \leq K,i\leq N}}. Each element \eqn{i\leq N} of the list, is a list of \eqn{k\leq K} data.frame with time \eqn{(t_{kij})_{j\leq n_{ik}}} and observations \eqn{(Z_{kij})_{j\leq n_{ik}}}. Each data.frame is named with the observation model identifiers.
#' @param Serr named vector of the estimated error mocel constants \eqn{(\varsigma_p)_{p\leq P}} with observation model identifiers as names.
#' @param Rerr named vector of the estimated error mocel constants \eqn{(\sigma_k)_{k\leq K}} with observation model identifiers as names.
#' @param ObsModel.transfo list containing two lists of transformations and two vectors linking each transformations to their observation model name in the Monolix project. The list should include identity transformations and be named \code{S} and \code{R}. The two vectors should be named \code{linkS} and \code{linkR}.
#'
#' Both \code{S} (for the direct observation models) and \code{linkS}, as well as \code{R} (for latent process models) and \code{linkR}, must have the same length.
#'
#' \itemize{
#' \item\code{S}: a list of transformations for the direct observation models. Each transformation corresponds to a variable \eqn{Y_p=h_p(S_p)}, where the name indicates which dynamic is observed (from \code{dynFUN}); \item\code{linkS} : a vector specifying the observation model names (that is used in the monolix project, \code{alpha1}, etc.) for each transformation, in the same order as in \code{S};
#'
#' \item\code{R}: similarly, a list of transformations for the latent process models. Although currently there is only one latent dynamic, each \eqn{s_k, k\leq K} transformation corresponds to the same dynamic but may vary for each \eqn{Y_k} observed. The names should match the output from \code{dynFUN}; \item \code{linkR} : a vector specifying the observation model names for each transformation, in the same order as in \code{R}.
#' }
#' @param data output from \code{\link{readMLX}} containing parameters "\code{mu}", "\code{Omega}", "\code{theta}", "\code{alpha1}", "\code{covariates}", "\code{ParModel.transfo}", "\code{ParModel.transfo.inv}", "\code{Sobs}", "\code{Robs}", "\code{Serr}", "\code{Rerr}", "\code{ObsModel.transfo}" extract from a monolix project.
#' @param n number of points per dimension to use for the Gauss-Hermite quadrature rule.
#' @param prune integer between 0 and 1, percentage of pruning for the Gauss-Hermite quadrature rule (default NULL).
#' @param parallel logical, if computation should be done in parallel.
#' @param ncores number of cores to use for parallelization, default will detect the number of cores available.
#' @param onlyLL logical, if only the log-likelihood should be computed (and not \eqn{\partial_{\alpha_1} LL} or \eqn{\partial_{\alpha_1}^2 LL}).
#' @param verbose logical, if progress bar should be printed through the computation.
#'
#' @return A list with the approximation by Gauss-Hermite quadrature of the likelihood \code{L}, the log-likelihood \code{LL}, the gradient of the log-likelihood \code{dLL}, and the Hessian of the log-likelihood \code{ddLL} at the point \eqn{\theta, \alpha} provided.
#' @export
#'
#' @examples
#'\dontrun{
#' project <- getMLXdir()
#'
#'
#' ObsModel.transfo = list(S=list(AB=log10),
#' linkS="yAB",
#' R=rep(list(S=function(x){x}),5),
#' linkR = paste0("yG",1:5))
#'
#' alpha=list(alpha0=NULL,
#' alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))
#'
#' data <- readMLX(project,ObsModel.transfo,alpha)
#'
#' LL <- gh.LL(dynFUN = dynFUN_demo,
#' y = c(S=5,AB=1000),
#' ObsModel.transfo=ObsModel.transfo,
#' data = data)
#'
#' print(LL)
#' }
gh.LL <- function(
dynFUN,
y,
mu=NULL,
Omega=NULL,
theta=NULL,
alpha1=NULL,
covariates=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs=NULL,
Robs=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
n=NULL,
prune=NULL,
parallel=TRUE,
ncores=NULL,
onlyLL=FALSE,
verbose=TRUE){
if(is.null(data)){
test <- sapply(c("mu","Omega","theta","alpha1","covariates","ParModel.transfo","ParModel.transfo.inv","Sobs","Robs","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu) || is.null(Omega)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
for(d in 1:length(data)){
test <- eval(parse(text=(paste0("is.null(",names(data)[d],")"))))
if(test){
eval(parse(text=paste0(names(data[d]),"<- data[[d]]")))
}
}
}
if(is.null(n)){
if(length(theta$psi_pop)==1){
n <- 100
}else if(length(theta$psi_pop)==2){
n <- 10
}else{
n <- 7
}
}
if(parallel){
if(!is.null(ncores)){
cluster <- snow::makeCluster(ncores)
}else{
cluster <- snow::makeCluster(parallel::detectCores())
}
doSNOW::registerDoSNOW(cluster)
}
N=length(mu)
i = 1
ntasks <- N
if(verbose){
pb <- txtProgressBar(max = ntasks, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)
}else{
opts <- NULL
}
res = foreach::foreach(i = 1:N,.packages = "REMixed",.options.snow=opts)%dopar%{
if(0 %in% diag(Omega[[i]])){
diag(Omega[[i]])[diag(Omega[[i]])==0] <- 10**(-5)
}
LLi = gh.LL.ind(mu_i = mu[[i]],
Omega_i = Omega[[i]],
theta = theta,
alpha1 = alpha1,
dynFUN = dynFUN,
y = y,
covariates_i = covariates[i,,drop=F],
ParModel.transfo = ParModel.transfo,
ParModel.transfo.inv = ParModel.transfo.inv,
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==i,]}),
Robs_i = lapply(Robs,FUN=function(R){R[R$id==i,]}),
Serr = Serr,
Rerr = Rerr,
ObsModel.transfo = ObsModel.transfo,
ind = i,
n = n,
prune = prune,
onlyLL = onlyLL)
return(LLi)
}
L = prod(sapply(res,FUN=function(ri){ri$Li}))
LL = sum(sapply(res,FUN=function(ri){ri$LLi}))
if(verbose)
close(pb)
if(parallel){
snow::stopCluster(cluster)
}
if(!onlyLL){
if(length(alpha1)!=1){
dLL = rowSums(sapply(res,FUN=function(ri){ri$dLLi}))
ddLL = matrix(rowSums(sapply(res,FUN=function(ri){ri$ddLLi})),ncol=length(dLL))
}else{
dLL = sum(sapply(res,FUN=function(ri){ri$dLLi}))
ddLL = matrix(sum(sapply(res,FUN=function(ri){ri$ddLLi})),ncol=length(dLL))
}
return(list(L=L,LL=LL,dLL=dLL,ddLL=ddLL))
}else{
return(LL)
}
}
gh.LL.ind <- function(
dynFUN,
y,
mu_i=NULL,
Omega_i=NULL,
theta=NULL,
alpha1=NULL,
covariates_i=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs_i=NULL,
Robs_i=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
ind = NULL,
n = NULL,
prune=NULL,
onlyLL=FALSE){
mu <- Omega <- Sobs <- Robs <- covariates <- NULL
if(is.null(data)){
test <- sapply(c("mu_i","Omega_i","theta","alpha1","covariates_i","ParModel.transfo","ParModel.transfo.inv","Sobs_i","Robs_i","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu_i) || is.null(Omega_i)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
argsNONind = c("theta","alpha1","ParModel.transfo","ParModel.transfo.inv","Serr","Rerr","ObsModel.transfo")
for(d in 1:length(argsNONind)){
test <- eval(parse(text=paste0("is.null(",argsNONind[d],")")))
if(test){
eval(parse(text=paste0(argsNONind[d],"<- data[[argsNONind[d]]]")))
}
}
argsInd = c("mu","Omega","Robs","Sobs","covariates")
for(d in 1:length(argsInd)){
test <- eval(parse(text=paste0("is.null(",paste0(argsInd[d],"_i"),")")))
if(test){
eval(parse(text=paste0(argsInd[d],"<- data[[argsInd[d]]]")))
}
}
if((is.null(mu_i) || is.null(Omega_i) || is.null(Robs_i) || is.null(Sobs_i) || is.null(covariates_i)) & is.null(ind)){
stop("Please provide individual information (mu_i,Omega_i,Sobs_i,Robs_i,covariates_i), or the individual id ind.")
}
if(is.null(mu_i)){
mu_i = mu[[ind]]
}
if(is.null(Omega_i)){
Omega_i = Omega[[ind]]
}
if(is.null(Sobs_i)){
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==ind,]})
}
if(is.null(Robs_i)){
Robs_i = lapply(Robs,FUN=function(R){R[R$id==ind,]})
}
if(is.null(covariates_i)){
covariates_i = covariates[ind,,drop=F]
}
}
if(is.null(n)){
if(length(theta$psi_pop)==1){
n <- 100
}else if(length(theta$psi_pop)==2){
n <- 10
}else{
n <- 7
}
}
dm = length(mu_i)
if(dm!=1){
mh.parm <- amgauss.hermite(n,mu=mu_i,Omega=Omega_i,prune=prune)
detsq.omega = prod(theta$omega)
root.omega = diag(1/theta$omega**2)
}else{
mh.parm <- agauss.hermite(n,mu = mu_i,sd=sqrt(as.numeric(Omega_i)))
detsq.omega = as.numeric(theta$omega)
root.omega = 1/as.numeric(theta$omega)**2
}
nd = length(mh.parm$Weights)
R.sz = length(Rerr)
S.sz = length(Serr)
all.tobs = sort(union(unlist(lapply(Robs_i,FUN=function(x){x$time})),unlist(lapply(Sobs_i,FUN=function(x){x$time}))))
# Need to compute, for each eta_i x individual i, the dynamics of the model
# eta_i = split(mh.parm$Points,1:nd)[[1]]
dyn <- setNames(lapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
PSI_i = indParm(theta[c("phi_pop","psi_pop","gamma","beta")],covariates_i,setNames(eta_i,colnames(Omega_i)),ParModel.transfo,ParModel.transfo.inv)
dyn_eta_i <- dynFUN(all.tobs,y,unlist(unname(PSI_i)))
return(dyn_eta_i)
}),paste0("eta_",1:nd))
# Compute marginal latent density for each individual and value of RE
R.margDensity <- setNames(lapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
# .Machine$double.xmin
decomp_marg = lapply(1:R.sz,FUN=function(k){
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
aux <- 1/(sig*sqrt(2*pi))*exp(-1/2*((Zki-a0-a1*Rki)/sig)**2)
if(any(aux==0)){
pw <- Rmpfr::mpfr(-1/2*((Zki-a0-a1*Rki)/sig)**2,precBits = 32)
aux = 1/(sig*sqrt(2*pi))*exp(pw)
decomp_aux = setNames(sapply(decomp(prod(aux)),as.numeric),c("exponent","mantissa"))
}else{
decomp_aux <- lapply(aux,decomp)
decomp_intermediaire = decomp(prod(sapply(decomp_aux,FUN=function(decomp){decomp["mantissa"]})))
mantissa_aux = decomp_intermediaire["mantissa"]
exponent_aux = sum(sapply(decomp_aux,FUN=function(decomp){decomp["exponent"]})) + decomp_intermediaire["exponent"]
decomp_aux <- c(exponent_aux,mantissa_aux)
}
return(decomp_aux)
})
if(any(sapply(decomp_marg,function(x){x["mantissa"]})==0)){
stop("[Error] in log-likelihood computation, latent part of the marginal likelihood is equal to zero.")
}
# decomp_marg = lapply(little_marg,decomp)
decomp_intermediaire = decomp(prod(sapply(decomp_marg,FUN=function(decomp){decomp["mantissa"]})))
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = sum(sapply(decomp_marg,FUN=function(decomp){decomp["exponent"]})) + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
}),paste0("eta_",1:nd))
# Compute marginal density for each individual and value of RE
if(S.sz!=0){
S.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
eta_i = split(mh.parm$Points,1:nd)[[ei]]
dyn.ei = dyn[[paste0("eta_",ei)]]
res = prod(sapply(1:S.sz,FUN=function(p){
Yp = ObsModel.transfo$linkS[p]
sig = Serr[[Yp]]
tpi = Sobs_i[[Yp]]$time
Ypi = Sobs_i[[Yp]][,Yp]
trs = ObsModel.transfo$S[which(ObsModel.transfo$linkS==Yp)]
Spi = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tpi,names(trs)])
prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Ypi-Spi)/sig)**2))
}))*(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}else{
S.margDensity = setNames(sapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}
# Compute individual log-Likelihood
Li.aux = lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(Li.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(Li.aux,function(li){
return(li[["mantissa"]]*10**(li[["exponent"]] - max_exponent))
})))
Li = c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
# sum(mh.parm$Weights*R.margDensity*S.margDensity)
invLi.aux = decomp(1/Li[["mantissa"]])
invLi = c(exponent=-Li[["exponent"]]+invLi.aux[["exponent"]],mantissa=invLi.aux[["mantissa"]])
# ln(p*10**o)=ln(p)+ln(10**o)=ln(p)+o*ln(10)
LLi = log(Li[["mantissa"]])+Li[["exponent"]]*log(10)
if(!onlyLL){
# Compute gradient of individual log-Likelihood
# dim.alpha = length(alpha1) # = K = length(Robs_i )
dfact = lapply(1:R.sz,FUN=function(k){
f = setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))
}),paste0("eta_",1:nd))
})
dLLi = sapply(dfact,FUN=function(f){
res.aux <- lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*f[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(res.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(res.aux,function(ri){
return(ri[["mantissa"]]*10**(ri[["exponent"]] - max_exponent))
})))
aux =c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
return(invLi[["mantissa"]]*aux[["mantissa"]]*10**(invLi[["exponent"]]+aux[["exponent"]]))
# 1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
})
# Compute hessienne of individual log-Likelihood
ddLLi = matrix(0,ncol=R.sz,nrow=R.sz)
# diagonal term
ddfact = lapply(1:R.sz,FUN=function(k){
setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(
(1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))**2 - 1/(sig**2)*(sum(Rki**2))
)
}),paste0("eta_",1:nd))
})
diag(ddLLi) = sapply(ddfact,FUN=function(f){
res.aux <- lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*f[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(res.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(res.aux,function(ri){
return(ri[["mantissa"]]*10**(ri[["exponent"]] - max_exponent))
})))
aux =c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
return(invLi[["mantissa"]]*aux[["mantissa"]]*10**(invLi[["exponent"]]+aux[["exponent"]]))
# return(1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
})-(dLLi)**2
# non diagonal term
for(k1 in 1:R.sz){
for(k2 in 1:R.sz){
if(k1<k2){
ddfact = setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
prod(sapply(c(k1,k2),FUN=function(k){
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(
1/sig**2*(sum(Rki*(Zki-a0-a1*Rki)))
)
}))
}),paste0("eta_",1:nd))
# ddLLi[k1,k2] <- ddLLi[k2,k1] <- 1/Li*sum(mh.parm$Weights*ddfact*R.margDensity*S.margDensity)-prod(dLLi[c(k1,k2)])
res.aux <- lapply(1:nd,function(ei){
decomp_intermediaire = decomp(mh.parm$Weights[[ei]]*ddfact[[ei]]*R.margDensity[[ei]][["mantissa"]]*S.margDensity[[ei]])
mantissa_res = decomp_intermediaire["mantissa"]
exponent_res = R.margDensity[[ei]][["exponent"]] + decomp_intermediaire["exponent"]
return(c(exponent_res,mantissa_res))
})
max_exponent <- max(sapply(res.aux, function(x) x[["exponent"]]))
decomp.aux <- decomp(sum(sapply(res.aux,function(ri){
return(ri[["mantissa"]]*10**(ri[["exponent"]] - max_exponent))
})))
aux =c(exponent=(max_exponent+decomp.aux[["exponent"]]),mantissa=decomp.aux[["mantissa"]])
ddLLi[k1,k2] <- ddLLi[k2,k1] <- (invLi[["mantissa"]]*aux[["mantissa"]]*10**(invLi[["exponent"]]+aux[["exponent"]])) - prod(dLLi[c(k1,k2)])
}
}
}
}else{
dLLi <- ddLLi <- NULL
}
return(list(Li=Li[["mantissa"]]*10**(Li[["exponent"]]),LLi=LLi,dLLi=dLLi,ddLLi=ddLLi))
}
# One dimensional Adaptative Hermite Gauss ------------------------------------------
agauss.hermite <- function(n,mu=0,sd=1){
gh <- fastGHQuad::gaussHermiteData(n)
names(gh) <- c("Points","Weights")
ww <- gh$Weights * sqrt(2) * sd * exp(gh$Points**2)
xx <- mu + sd * sqrt(2) * gh$Points
return(data.frame(Points=xx,Weights=ww))
}
# Multi-dimensional Adaptative Hermite Gauss -----------------------------------------
amgauss.hermite <- function(n,mu=rep(0,ncol(Omega)),Omega=diag(rep(1,length(mu))),prune=NULL){
dm = length(mu)
gh <- as.data.frame(fastGHQuad::gaussHermiteData(n))
idx <- as.matrix(expand.grid(rep(list(1:n),dm)))
pts <- matrix(gh[idx,1],nrow(idx),dm)
wts <- apply(matrix(gh[idx,2],nrow(idx),dm), 1, prod)
# plot(pts, cex=-5/log(wts), pch=19)
if(!is.null(prune)) {
qwt <- quantile(wts, probs=prune)
pts <- pts[wts > qwt,]
wts <- wts[wts > qwt]
}
eig <- eigen(Omega)
S = eig$vectors
L = diag(sqrt(eig$values))
rot <- S %*% L
ppts <-t( sqrt(2) * rot %*% t(pts) + mu)
wwts <- wts * sqrt(2)**dm * abs(det(rot)) * exp(apply(pts,1,FUN=function(x){t(x)%*%x}))
return(list(Points=ppts, Weights=wwts))
}
# individual contribution to log-lik gradient in alpha1=0 ----------------------------
lambda.max.ind <- function(
dynFUN,
y,
mu_i=NULL,
Omega_i=NULL,
theta=NULL,
alpha1=NULL,
covariates_i=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs_i=NULL,
Robs_i=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
ind = NULL,
n = NULL,
prune=NULL){
mu <- Omega <- Sobs <- Robs <- covariates <- NULL
if(is.null(data)){
test <- sapply(c("mu_i","Omega_i","theta","alpha1","covariates_i","ParModel.transfo","ParModel.transfo.inv","Sobs_i","Robs_i","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu_i) || is.null(Omega_i)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
argsNONind = c("theta","alpha1","ParModel.transfo","ParModel.transfo.inv","Serr","Rerr","ObsModel.transfo")
for(d in 1:length(argsNONind)){
test <- eval(parse(text=paste0("is.null(",argsNONind[d],")")))
if(test){
eval(parse(text=paste0(argsNONind[d],"<- data[[argsNONind[d]]]")))
}
}
argsInd = c("mu","Omega","Robs","Sobs","covariates")
for(d in 1:length(argsInd)){
test <- eval(parse(text=paste0("is.null(",paste0(argsInd[d],"_i"),")")))
if(test){
eval(parse(text=paste0(argsInd[d],"<- data[[argsInd[d]]]")))
}
}
if((is.null(mu_i) || is.null(Omega_i) || is.null(Robs_i) || is.null(Sobs_i) || is.null(covariates_i)) & is.null(ind)){
stop("Please provide individual information (mu_i,Omega_i,Sobs_i,Robs_i,covariates_i), or the individual id ind.")
}
if(is.null(mu_i)){
mu_i = mu[[ind]]
}
if(is.null(Omega_i)){
Omega_i = Omega[[ind]]
}
if(is.null(Sobs_i)){
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==ind,]})
}
if(is.null(Robs_i)){
Robs_i = lapply(Robs,FUN=function(R){R[R$id==ind,]})
}
if(is.null(covariates_i)){
covariates_i = covariates[ind,,drop=F]
}
}
if(is.null(n)){
n <- floor(100**(1/length(theta$psi_pop)))
}
dm = length(mu_i)
if(dm!=1){
mh.parm <- amgauss.hermite(n,mu=mu_i,Omega=Omega_i,prune=prune)
detsq.omega = prod(theta$omega)
root.omega = diag(1/theta$omega**2)
}else{
mh.parm <- agauss.hermite(n,mu = mu_i,sd=sqrt(as.numeric(Omega_i)))
detsq.omega = as.numeric(theta$omega)
root.omega = 1/as.numeric(theta$omega)**2
}
nd = length(mh.parm$Weights)
R.sz = length(Rerr)
S.sz = length(Serr)
all.tobs = sort(union(unlist(lapply(Robs_i,FUN=function(x){x$time})),unlist(lapply(Sobs_i,FUN=function(x){x$time}))))
dyn <- setNames(lapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
PSI_i = indParm(theta[c("phi_pop","psi_pop","gamma","beta")],covariates_i,setNames(eta_i,colnames(Omega_i)),ParModel.transfo,ParModel.transfo.inv)
dyn_eta_i <- dynFUN(all.tobs,y,unlist(unname(PSI_i)))
return(dyn_eta_i)
}),paste0("eta_",1:nd))
# Compute marginal latent density for each individual and value of RE
R.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
res = prod(sapply(1:R.sz,FUN=function(k){
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Zki-a0-a1*Rki)/sig)**2))
}))
}),paste0("eta_",1:nd))
# A(Y|theta0)
if(S.sz!=0){
S.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
eta_i = split(mh.parm$Points,1:nd)[[ei]]
dyn.ei = dyn[[paste0("eta_",ei)]]
res = prod(sapply(1:S.sz,FUN=function(p){
Yp = ObsModel.transfo$linkS[p]
sig = Serr[[Yp]]
tpi = Sobs_i[[Yp]]$time
Ypi = Sobs_i[[Yp]][,Yp]
trs = ObsModel.transfo$S[which(ObsModel.transfo$linkS==Yp)]
Spi = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tpi,names(trs)])
prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Ypi-Spi)/sig)**2))
}))*(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}else{
S.margDensity = setNames(sapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
}),paste0("eta_",1:nd))
}
# s'k(Ri(krij))A(Y|theta0)
S2.margDensity <- lapply(1:R.sz,FUN=function(k){
setNames(sapply(1:nd,FUN=function(ei){
dyn.ei = dyn[[paste0("eta_",ei)]]
yGk = ObsModel.transfo$linkR[k]
sig = Rerr[[yGk]]
tki = Robs_i[[yGk]]$time
Zki = Robs_i[[yGk]][,yGk]
a0 = theta$alpha0[[yGk]]
a1 = alpha1[[yGk]]
trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
return(1/sig**2*sum((Zki-a0)*Rki)*S.margDensity[[ei]])
}),paste0("eta_",1:nd))})
indi.contrib <- sapply(1:R.sz,FUN=function(k){
sum(mh.parm$Weights*S2.margDensity[[k]])
})/sum(mh.parm$Weights*S.margDensity)
return(indi.contrib)
}
# log-lik gradient in alpha1=0 for maximum penalization parameter to reach ----------
lambda.max <- function(
dynFUN,
y,
mu=NULL,
Omega=NULL,
theta=NULL,
alpha1=NULL,
covariates=NULL,
ParModel.transfo=NULL,
ParModel.transfo.inv=NULL,
Sobs=NULL,
Robs=NULL,
Serr=NULL,
Rerr=NULL,
ObsModel.transfo=NULL,
data=NULL,
n = NULL,
prune=NULL,
parallel = TRUE,
ncores=NULL,
onlyLL=FALSE,
verbose=TRUE){
if(is.null(data)){
test <- sapply(c("mu","Omega","theta","alpha1","covariates","ParModel.transfo","ParModel.transfo.inv","Sobs","Robs","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
if(any(test))
stop("Please provide all necessary arguments.")
}else{
if((is.null(mu) || is.null(Omega)) & !all(c("mu","Omega") %in% names(data))){
stop("Please provide mu and Omega if these are missing from data.")
}
for(d in 1:length(data)){
test <- eval(parse(text=(paste0("is.null(",names(data)[d],")"))))
if(test){
eval(parse(text=paste0(names(data[d]),"<- data[[d]]")))
}
}
}
if(is.null(n)){
n <- floor(100**(1/length(theta$psi_pop)))
}
if(parallel){
if(!is.null(ncores)){
cluster <- parallel::makeCluster(ncores)
}else{
cluster <- parallel::makeCluster(parallel::detectCores())
}
doSNOW::registerDoSNOW(cluster)
}
N=length(mu)
i = 1
if(verbose){
ntasks <- N
pb <- txtProgressBar(max = ntasks, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)
}else{
opts <- NULL
}
res = foreach::foreach(i = 1:N,.packages = "REMixed",.export = c("lambda.max.ind","amgauss.hermite"),.options.snow=opts)%dopar%{
if(0 %in% diag(Omega[[i]])){
diag(Omega[[i]])[diag(Omega[[i]])==0] <- 10**(-5)
}
lambda.max.ind(mu_i = mu[[i]],
Omega_i = Omega[[i]],
theta = theta,
alpha1 = alpha1,
dynFUN = dynFUN,
y = y,
covariates_i = covariates[i,,drop=F],
ParModel.transfo = ParModel.transfo,
ParModel.transfo.inv = ParModel.transfo.inv,
Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==i,]}),
Robs_i = lapply(Robs,FUN=function(R){R[R$id==i,]}),
Serr = Serr,
Rerr = Rerr,
ObsModel.transfo = ObsModel.transfo,
n = n,
prune = prune)
}
if(verbose)
close(pb)
if(parallel)
snow::stopCluster(cluster)
return(max(Reduce("+",res)))
}
decomp <- function(x){
if(x==0){
return(c(exponent=0,mantissa=0))
}
exponent_x <- floor(log10(abs(x))) # Exponent
exponent_lim <- abs(floor(log10(.Machine$double.xmin)))
if(abs(exponent_x)>exponent_lim){
mantissa_x <- x*10**(-sign(exponent_x)*(exponent_lim-1))
mantissa_x <- mantissa_x*10**(-(exponent_x+(exponent_lim-1)))
}else{
mantissa_x <- x*10**(-exponent_x)
}
return(c(exponent=exponent_x,mantissa=mantissa_x))
}
#
# gh.LL.ind_OLD <- function(
# dynFUN,
# y,
# mu_i=NULL,
# Omega_i=NULL,
# theta=NULL,
# alpha1=NULL,
# covariates_i=NULL,
# ParModel.transfo=NULL,
# ParModel.transfo.inv=NULL,
# Sobs_i=NULL,
# Robs_i=NULL,
# Serr=NULL,
# Rerr=NULL,
# ObsModel.transfo=NULL,
# data=NULL,
# ind = NULL,
# n = NULL,
# prune=NULL,
# onlyLL=FALSE){
#
# mu <- Omega <- Sobs <- Robs <- covariates <- NULL
#
# if(is.null(data)){
# test <- sapply(c("mu_i","Omega_i","theta","alpha1","covariates_i","ParModel.transfo","ParModel.transfo.inv","Sobs_i","Robs_i","Serr","Rerr","ObsModel.transfo"),FUN=is.null)
# if(any(test))
# stop("Please provide all necessary arguments.")
# }else{
# if((is.null(mu_i) || is.null(Omega_i)) & !all(c("mu","Omega") %in% names(data))){
# stop("Please provide mu and Omega if these are missing from data.")
# }
# argsNONind = c("theta","alpha1","ParModel.transfo","ParModel.transfo.inv","Serr","Rerr","ObsModel.transfo")
# for(d in 1:length(argsNONind)){
# test <- eval(parse(text=paste0("is.null(",argsNONind[d],")")))
# if(test){
# eval(parse(text=paste0(argsNONind[d],"<- data[[argsNONind[d]]]")))
# }
# }
# argsInd = c("mu","Omega","Robs","Sobs","covariates")
# for(d in 1:length(argsInd)){
# test <- eval(parse(text=paste0("is.null(",paste0(argsInd[d],"_i"),")")))
# if(test){
# eval(parse(text=paste0(argsInd[d],"<- data[[argsInd[d]]]")))
#
# }
# }
# if((is.null(mu_i) || is.null(Omega_i) || is.null(Robs_i) || is.null(Sobs_i) || is.null(covariates_i)) & is.null(ind)){
# stop("Please provide individual information (mu_i,Omega_i,Sobs_i,Robs_i,covariates_i), or the individual id ind.")
# }
# if(is.null(mu_i)){
# mu_i = mu[[ind]]
# }
# if(is.null(Omega_i)){
# Omega_i = Omega[[ind]]
# }
# if(is.null(Sobs_i)){
# Sobs_i = lapply(Sobs,FUN=function(S){S[S$id==ind,]})
# }
# if(is.null(Robs_i)){
# Robs_i = lapply(Robs,FUN=function(R){R[R$id==ind,]})
# }
# if(is.null(covariates_i)){
# covariates_i = covariates[ind,,drop=F]
# }
# }
#
# if(is.null(n)){
# n <- floor(100**(1/length(theta$psi_pop)))
# }
#
# dm = length(mu_i)
#
# if(dm!=1){
# mh.parm <- amgauss.hermite(n,mu=mu_i,Omega=Omega_i,prune=prune)
# detsq.omega = prod(theta$omega)
# root.omega = diag(1/theta$omega**2)
# }else{
# mh.parm <- agauss.hermite(n,mu = mu_i,sd=sqrt(as.numeric(Omega_i)))
# detsq.omega = as.numeric(theta$omega)
# root.omega = 1/as.numeric(theta$omega)**2
# }
#
# nd = length(mh.parm$Weights)
# R.sz = length(Rerr)
# S.sz = length(Serr)
# all.tobs = sort(union(unlist(lapply(Robs_i,FUN=function(x){x$time})),unlist(lapply(Sobs_i,FUN=function(x){x$time}))))
#
# # Need to compute, for each eta_i x individual i, the dynamics of the model
# # eta_i = split(mh.parm$Points,1:nd)[[1]]
# dyn <- setNames(lapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
# PSI_i = indParm(theta[c("phi_pop","psi_pop","gamma","beta")],covariates_i,setNames(eta_i,colnames(Omega_i)),ParModel.transfo,ParModel.transfo.inv)
# dyn_eta_i <- dynFUN(all.tobs,y,unlist(unname(PSI_i)))
#
# return(dyn_eta_i)
# }),paste0("eta_",1:nd))
#
# # Compute marginal latent density for each individual and value of RE
# R.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
#
#
# res = prod(sapply(1:R.sz,FUN=function(k){
# yGk = ObsModel.transfo$linkR[k]
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Zki-a0-a1*Rki)/sig)**2))
# }))
# }),paste0("eta_",1:nd))
#
# # Compute marginal density for each individual and value of RE
# if(S.sz!=0){
# S.margDensity <- setNames(sapply(1:nd,FUN=function(ei){
# eta_i = split(mh.parm$Points,1:nd)[[ei]]
# dyn.ei = dyn[[paste0("eta_",ei)]]
#
# res = prod(sapply(1:S.sz,FUN=function(p){
# Yp = ObsModel.transfo$linkS[p]
# sig = Serr[[Yp]]
# tpi = Sobs_i[[Yp]]$time
# Ypi = Sobs_i[[Yp]][,Yp]
#
# trs = ObsModel.transfo$S[which(ObsModel.transfo$linkS==Yp)]
# Spi = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tpi,names(trs)])
#
# prod(1/(sig*sqrt(2*pi))*exp(-1/2*((Ypi-Spi)/sig)**2))
# }))*(1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
# }),paste0("eta_",1:nd))
# }else{
# S.margDensity = setNames(sapply(split(mh.parm$Points,1:nd),FUN=function(eta_i){
# (1/((2*pi)**(dm/2)*detsq.omega)*exp(-1/2*eta_i%*%root.omega%*%eta_i))
# }),paste0("eta_",1:nd))
# }
#
# # Compute individual log-Likelihood
# Li = sum(mh.parm$Weights*R.margDensity*S.margDensity)
# LLi = log(Li)
#
# if(!onlyLL){
# # Compute gradient of individual log-Likelihood
# # dim.alpha = length(alpha1) # = K = length(Robs_i )
# dfact = lapply(1:R.sz,FUN=function(k){
# setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
# yGk = ObsModel.transfo$linkR[k]
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# return(1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))
# }),paste0("eta_",1:nd))
# })
#
# dLLi = sapply(dfact,FUN=function(f){
# return(1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
# })
#
# # Compute hessienne of individual log-Likelihood
# ddLLi = matrix(0,ncol=R.sz,nrow=R.sz)
# # diagonal term
# ddfact = lapply(1:R.sz,FUN=function(k){
# setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
# yGk = ObsModel.transfo$linkR[k]
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# return(
# (1/(sig**2)*(sum(Rki*(Zki-a0-a1*Rki))))**2 - 1/(sig**2)*(sum(Rki**2))
# )
# }),paste0("eta_",1:nd))
# })
#
# diag(ddLLi) = sapply(ddfact,FUN=function(f){
# return(1/Li*sum(mh.parm$Weights*f*R.margDensity*S.margDensity))
# })-(dLLi)**2
#
# # non diagonal term
# for(k1 in 1:R.sz){
# for(k2 in 1:R.sz){
# if(k1<k2){
# ddfact = setNames(sapply(1:nd,FUN=function(ei){
# dyn.ei = dyn[[paste0("eta_",ei)]]
# prod(sapply(c(k1,k2),FUN=function(k){
# yGk = ObsModel.transfo$linkR[k]
#
# sig = Rerr[[yGk]]
# tki = Robs_i[[yGk]]$time
# Zki = Robs_i[[yGk]][,yGk]
#
# a0 = theta$alpha0[[yGk]]
# a1 = alpha1[[yGk]]
# trs = ObsModel.transfo$R[which(ObsModel.transfo$linkR==yGk)]
# Rki = trs[[1]](dyn.ei[dyn.ei[,"time"] %in% tki,names(trs)])
#
# return(
# 1/sig**2*(sum(Rki*(Zki-a0-a1*Rki)))
# )
# }))
# }),paste0("eta_",1:nd))
#
# ddLLi[k1,k2] <- ddLLi[k2,k1] <- 1/Li*sum(mh.parm$Weights*ddfact*R.margDensity*S.margDensity)-prod(dLLi[c(k1,k2)])
# }
# }
# }
# }else{
# dLLi <- ddLLi <- NULL
# }
# return(list(Li=Li,LLi=LLi,dLLi=dLLi,ddLLi=ddLLi))
# }
Try the REMixed package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.