Nothing
########################### Fisher Information Matrix & LL by linearisation #############################
#' Computes the Fisher Information Matrix by linearisation
#'
#' Estimate by linearisation the Fisher Information Matrix and the standard
#' error of the estimated parameters.
#'
#' The inverse of the Fisher Information Matrix provides an estimate of the
#' variance of the estimated parameters theta. This matrix cannot be computed
#' in closed-form for nonlinear mixed-effect models; instead, an approximation
#' is obtained as the Fisher Information Matrix of the Gaussian model deduced
#' from the nonlinear mixed effects model after linearisation of the function f
#' around the conditional expectation of the individual Gaussian parameters.
#' This matrix is a block matrix (no correlations between the estimated fixed
#' effects and the estimated variances).
#'
#' @param saemixObject an object returned by the \code{\link{saemix}} function
#' @return The function returns an updated version of the object saemix.fit in
#' which the following elements have been added: \describe{
#' \item{se.fixed:}{standard error of fixed effects, obtained as part of the
#' diagonal of the inverse of the Fisher Information Matrix (only when
#' fim.saemix has been run, or when saemix.options$algorithms\[2\] is 1)}
#' \item{se.omega:}{standard error of the variance of random effects, obtained
#' as part of the diagonal of the inverse of the Fisher Information Matrix
#' (only when fim.saemix has been run, or when the saemix.options$algorithms\[2\]
#' is 1)}
#' \item{se.res:}{standard error of the parameters of the residual error
#' model, obtained as part of the diagonal of the inverse of the Fisher
#' Information Matrix (only when fim.saemix has been run, or when the
#' saemix.options$algorithms\[2\] is 1)}
#' \item{fim:}{Fisher Information Matrix}
#' \item{ll.lin:}{ likelihood calculated by linearisation} }
#' @author Emmanuelle Comets <emmanuelle.comets@@inserm.fr>, Audrey Lavenu,
#' Marc Lavielle.
#' @seealso \code{\link{SaemixObject}},\code{\link{saemix}}
#' @references E Comets, A Lavenu, M Lavielle M (2017). Parameter estimation in nonlinear mixed effect models using saemix,
#' an R implementation of the SAEM algorithm. Journal of Statistical Software, 80(3):1-41.
#'
#' E Kuhn, M Lavielle (2005). Maximum likelihood estimation in nonlinear mixed effects models.
#' Computational Statistics and Data Analysis, 49(4):1020-1038.
#'
#' E Comets, A Lavenu, M Lavielle (2011). SAEMIX, an R version of the SAEM algorithm. 20th meeting of the
#' Population Approach Group in Europe, Athens, Greece, Abstr 2173.
#'
#' @keywords models
#' @examples
#'
#' # Running the main algorithm to estimate the population parameters
#' data(theo.saemix)
#'
#' saemix.data<-saemixData(name.data=theo.saemix,header=TRUE,sep=" ",na=NA,
#' name.group=c("Id"),name.predictors=c("Dose","Time"),
#' name.response=c("Concentration"),name.covariates=c("Weight","Sex"),
#' units=list(x="hr",y="mg/L",covariates=c("kg","-")), name.X="Time")
#'
#' model1cpt<-function(psi,id,xidep) {
#' dose<-xidep[,1]
#' tim<-xidep[,2]
#' ka<-psi[id,1]
#' V<-psi[id,2]
#' CL<-psi[id,3]
#' k<-CL/V
#' ypred<-dose*ka/(V*(ka-k))*(exp(-k*tim)-exp(-ka*tim))
#' return(ypred)
#' }
#'
#' saemix.model<-saemixModel(model=model1cpt,
#' description="One-compartment model with first-order absorption",
#' psi0=matrix(c(1.,20,0.5,0.1,0,-0.01),ncol=3, byrow=TRUE,
#' dimnames=list(NULL, c("ka","V","CL"))),transform.par=c(1,1,1),
#' covariate.model=matrix(c(0,1,0,0,0,0),ncol=3,byrow=TRUE),fixed.estim=c(1,1,1),
#' covariance.model=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE),
#' omega.init=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE), error.model="constant")
#'
#' saemix.options<-list(algorithm=c(1,0,0),seed=632545,save=FALSE,save.graphs=FALSE)
#'
#' # Not run (strict time constraints for CRAN)
#' # saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
#'
#' # Estimating the Fisher Information Matrix using the result of saemix
#' # & returning the result in the same object
#' # fim.saemix(saemix.fit)
#'
#'
#' @export fim.saemix
fim.saemix<-function(saemixObject) {
# Estimate the Fisher Information Matrix and the s.e. of the estimated parameters
saemix.model<-saemixObject["model"]
saemix.data<-saemixObject["data"]
saemix.res<-saemixObject["results"]
xind<-saemix.data["data"][,saemix.data["name.predictors"],drop=FALSE]
yobs<-saemix.data["data"][,saemix.data["name.response"]]
# covariance.model<-0*saemix.model["covariance.model"]
covariance.model<-saemix.model["covariance.model"]
omega<-saemix.res["omega"]
omega.null<-0*omega
# diag(covariance.model)<-mydiag(saemix.model["covariance.model"])
# omega<-0*saemix.res["omega"] # Why use only diag(omega) ???
# diag(omega)<-mydiag(saemix.res["omega"])
hat.phi<-saemix.res["cond.mean.phi"]
nphi<-dim(hat.phi)[2]
nomega<-sum(covariance.model[lower.tri(covariance.model,diag=TRUE)])
if (saemixObject["model"]["modeltype"]=="structural"){
nres<-length(saemix.res["indx.res"])
} else{
nres <- 0
}
nytype<-length(unique(saemix.data["data"]["ytype"]))
dphi<-cutoff(abs(colMeans(hat.phi))*1e-4,1e-10)
coefphi<-c(0,-1,1)
F<-array(data=0,dim=c(saemix.data["ntot.obs"],nphi,length(coefphi)))
gs<-matrix(0,saemix.data["ntot.obs"],4)
etype.exp<-which(saemix.model["error.model"]=='exponential')
for (l in 1:length(coefphi)) {
for (j in 1:nphi) {
phi<-hat.phi
phi[,j]<-phi[,j]+coefphi[l]*dphi[j]
psi<-transphi(phi,saemix.model["transform.par"])
f <- saemix.model["model"](psi, saemix.data["data"][,"index"],xind)
for(ityp in etype.exp) f[saemix.data["data"][,saemix.data["name.ytype"]]==ityp]<-log(cutoff(f[saemix.data["data"][,saemix.data["name.ytype"]]==ityp]))
F[,j,l]<-f
}
}
ind.covariates<-which(saemix.model["betaest.model"]>0)
f0<-F[,1,1]
# g0<-cutoff(saemix.res["respar"][1]+saemix.res["respar"][2]*abs(f0))
if (saemixObject["model"]["modeltype"]=="structural"){
g0<-error(f0,saemix.res@respar,saemix.data["data"][["ytype"]])
}
# DF<-(F[,,3]-F[,,2])/matrix(rep(dphi,each=saemix.data["ntot.obs"]), ncol=length(dphi))/2
DF<-(F[,,3]-F[,,1])/matrix(rep(dphi,each=saemix.data["ntot.obs"]), ncol=length(dphi)) #gradient of f (changed from F[,,2] to F[,,1])
z<-matrix(0,saemix.data["ntot.obs"],1)
invVi<-Gi<-list() # Individual variance matrices
j2<-0
for (i in 1:saemix.data["N"]) {
ni<-saemix.data["nind.obs"][i]
j1<-j2+1
j2<-j2+ni
z[j1:j2]<-yobs[j1:j2] - f0[j1:j2] + DF[j1:j2,,drop=FALSE]%*%hat.phi[i,]
if (saemixObject["model"]["modeltype"]=="structural"){
Vi<- DF[j1:j2,,drop=FALSE] %*% omega %*% t(DF[j1:j2,,drop=FALSE]) + mydiag((g0[j1:j2])^2, nrow=ni)
} else{
Vi<- DF[j1:j2,,drop=FALSE] %*% t(DF[j1:j2,,drop=FALSE])+ mydiag(1, nrow=ni)
}
# invVi[[i]]<-solve(Vi[[i]])
# Invert avoiding numerical problems
# Invert avoiding numerical problems
Gi[[i]]<-round(Vi*1e12)/1e12
VD<-try(eigen(Gi[[i]]))
if(inherits(VD,"try-error") || det(Gi[[i]])==0) {
cat("Unable to compute the FIM by linearisation.\n") # si matrice de variance non inversible
stop()
}
D<-Re(VD$values)
V<-Re(VD$vectors)
invVi[[i]] <- V%*%mydiag(1/D,nrow=length(D))%*%t(V)
}
# ECO ici modifie car role de covariate.estim pas clair
# covariate.estim=si un parametre (et ses covariables associees) sont estimees ou non
# covariate.estim<-matrix(rep(saemix.model["fixed.estim"], dim(saemix.model["betaest.model"])[1]),byrow=TRUE, ncol=length(saemix.model["fixed.estim"]))*saemix.model["betaest.model"] # 29/05/20
covariate.estim<-saemix.model["betaest.model"]
covariate.estim[1,]<-saemix.model["fixed.estim"]
j<-which(saemix.model["betaest.model"]>0)
ind.fixed.est<-(covariate.estim[j]>0)
npar<-sum(ind.fixed.est)
# Tracking indices for covariances
myidx.omega<-c()
myidx.cor<-c()
name.rand1<-name.rand2<-c()
myidx.track<-NULL
ipar<-npar
for(iom in 1:dim(covariance.model)[1]) {
for(jom in iom:dim(covariance.model)[1]) {
if(covariance.model[iom,jom]==1) {
ipar<-ipar+1
myidx.track<-rbind(myidx.track,c(ipar,iom,jom))
if(iom==jom) {
myidx.omega<-c(myidx.omega,ipar)
name.rand1<-c(name.rand1,paste("Var",saemixObject@model@name.modpar[iom],sep="."))
name.rand2<-c(name.rand2,paste("SD",saemixObject@model@name.modpar[iom],sep="."))
}
else {
myidx.cor<-c(myidx.cor,ipar)
name.rand1<-c(name.rand1,paste("Cov",saemixObject@model@name.modpar[iom],saemixObject@model@name.modpar[jom],sep="."))
name.rand2<-c(name.rand2,paste("Corr",saemixObject@model@name.modpar[iom],saemixObject@model@name.modpar[jom],sep="."))
}
}
}
}
if(length(myidx.cor)>0) {
track.var<-myidx.track[myidx.track[,1] %in% myidx.omega,]
for(i in myidx.cor) {
ij<-which(myidx.track[,1]==i)
myidx.track[ij,2]<-track.var[track.var[,2]==myidx.track[ij,2],1]
myidx.track[ij,3]<-track.var[track.var[,2]==myidx.track[ij,3],1]
}
}
if(saemixObject@model@modeltype=="structural")
namallpar<-c(saemixObject@results@name.fixed,name.rand1, saemixObject@results@name.sigma[saemixObject@results@indx.res], name.rand2) else
namallpar<-c(saemixObject@results@name.fixed,name.rand1, name.rand2)
# hw=waitbar(1,'Estimating the population parameters (SAEM). Wait...')
ll.lin<- -0.5*saemix.data["ntot.obs"]*log(2*pi)
j2<-0
# indMF<-list() # Individual FIM
MF<-matrix(0,nrow=(npar+nomega+nres),ncol=(npar+nomega+nres))
for (i in 1:saemix.data["N"]) {
#waitbar(i/N,hw)
ni<-saemix.data["nind.obs"][i]
j1<-j2+1
j2<-j2+ni
yi<-yobs[j1:j2]
DFi<-DF[j1:j2,,drop=FALSE]
f0i<-f0[j1:j2]
if (saemixObject["model"]["modeltype"]=="structural"){
g0i<-g0[j1:j2]
}
zi<-z[j1:j2]
Ai<-kronecker(diag(nphi),as.matrix(saemix.model["Mcovariates"][i,]))
Ai<-Ai[,ind.covariates,drop=FALSE]
DFAi<-DFi%*%Ai
Dzi<-zi-DFAi%*%saemix.res["betas"]
# Derivatives of Vi=var(yi) for subject i, w/r to lambda (FO approximation, neglecting dVi/dmu)
DV<-list()
for(ipar in 1:npar) {
DV[[ipar]]<-matrix(0,ncol=ni,nrow=ni)
}
for(iom in 1:dim(covariance.model)[1]) {
for(jom in iom:dim(covariance.model)[1]) {
if(covariance.model[iom,jom]==1) {
ipar<-ipar+1
domega<-omega.null
domega[iom,jom]<-domega[jom,iom]<-1
# if(iom==jom) domega[iom,jom]<-1*sqrt(omega[iom,jom]) else domega[iom,jom]<-1 # if parameterised in omega and not omega2,
if (saemixObject["model"]["modeltype"]=="structural"){
DV[[ipar]]<-DFi %*% domega %*% t(DFi)
} else {
DV[[ipar]]<-DFi %*% t(DFi)
}
}
}
}
# for(ipar in 1:nomega) {
# domega<-omega.null
# domega[ipar,ipar]<-sqrt(omega[ipar,ipar])*2
# DV[[ipar+npar]] <- DFi %*% t(DFi)
# }
if (saemixObject["model"]["modeltype"]=="structural"){
for(ipar.res in 1:(2*nytype)) {
if(!is.na(match(ipar.res,saemix.res@indx.res))) {
ipar<-ipar+1
if(ipar.res%%2 == 1) DV[[ipar]]<-mydiag(2*g0i, nrow=ni) else DV[[ipar]]<-mydiag(2*g0i*f0i, nrow=ni)
}
}
}
# for(ipar.res in 1:(2*nytype)) {
# if(!is.na(match(ipar.res,saemix.res@indx.res))) {
# ipar<-ipar+1
# if (saemixObject["model"]["modeltype"]=="structural"){
# if(ipar.res%%2 == 1) DV[[ipar]]<-mydiag(2*g0i, nrow=ni) else DV[[ipar]]<-mydiag(2*g0i*f0i, nrow=ni)
# } else{
# DV[[ipar]]<-mydiag(0, nrow=ni)
# }
# }
# }
# blocA <- t(DFAi) %*% invVi[[i]] %*% DFAi
if (sum(ind.fixed.est)>0) {
DFAiest<-DFAi[,ind.fixed.est,drop=FALSE]
blocA<-t(DFAiest)%*% invVi[[i]] %*%DFAiest
} else blocA<-NULL
# blocAbis<-matrix(0,ncol=(npar),nrow=(npar))
# for(ii in 1:npar) {
# for(ij in 1:npar) {
# blocAbis[ii,ij]<-DFi[,ii] %*% invVi[[i]] %*% DFi[,ij]
# }
# }
blocB<-matrix(0,ncol=(nomega+nres),nrow=(nomega+nres))
for(ij in 1:(nomega+nres)) { # columns
for(ii in 1:(nomega+nres)) { # lines, so that blocB is ordered according to c(covariance.model)
blocB[ii,ij]<-sum(diag(DV[[ii+npar]] %*% invVi[[i]] %*% DV[[ij+npar]] %*% invVi[[i]] ))/2
}
}
blocC<-matrix(0,ncol=(npar),nrow=(nomega+nres))
MFi <-rbind( cbind(blocA,t(blocC)), cbind(blocC, blocB))
MF <- MF+MFi
# FIMi[[i]]<-MFi
ll.lin <- ll.lin - 0.5*log(det(Gi[[i]])) - 0.5*t(Dzi)%*% invVi[[i]] %*%Dzi
}
for(ityp in etype.exp) ll.lin<-ll.lin-sum(yobs[saemix.data["data"][,saemix.data["name.ytype"]]==ityp])
if (sum(ind.fixed.est)>0) {
Mparam<-matrix(0,dim(saemix.model["betaest.model"])[1], dim(saemix.model["betaest.model"])[2])
Mparam[1,]<-saemix.model["transform.par"]
Mtp<-Mparam[saemix.model["betaest.model"]>0]
Mtp<-Mtp[ind.fixed.est]
dbetas <- dtransphi(saemix.res["betas"][ind.fixed.est],Mtp)
Mupth<-mydiag(1/dbetas,nrow=length(dbetas))
Fmu<-MF[1:npar,1:npar]
Fth<-t(Mupth)%*%Fmu%*%Mupth
MF[1:npar,1:npar]<-Fth
# Individual FIM
# for(i in 1:saemix.data["N"]) {
# Fmui<-t(Mupth) %*% indMF[[i]][1:npar,1:npar] %*% Mupth
# indMF[[i]][1:npar,1:npar]<-Fmui
# }
Cth<-try(solve(Fth))
if(inherits(Cth,"try-error")) {
if(saemixObject@options$warnings) cat("Error computing the Fisher Information Matrix: singular system.\n")
Cth<-NA*Fth
}
} else {
Cth<-NULL
}
fim<-MF
sTHest<-sqrt(mydiag(Cth))
#sTH<-matrix(0,1,length(saemix.res["betas"]))
sTH<-rep(0,length(saemix.res["betas"]))
sTH[ind.fixed.est]<-sTHest
se.fixed<-sTH
FO<-MF[-c(1:npar),-c(1:npar)]
CO<-try(solve(FO))
if(inherits(CO,"try-error")) {
CO<-NA*FO
if(saemixObject@options$warnings) cat("Error computing the Fisher Information Matrix: singular system.\n")
}
sO<-sqrt(mydiag(CO))
se.omega<-matrix(0,nphi,1)
se.sdcor<-se.cov<-matrix(0,nphi,nphi)
se.omega[saemix.model["indx.omega"]]<-sO[myidx.omega-npar]
se.res<-matrix(0,2*nytype,1)
if(saemixObject@model@modeltype=="structural") se.res[saemix.res["indx.res"]]<-sO[(nomega+1):length(sO)]
# Table with SE, CV and confidence intervals
estpar<-c(saemixObject@results@fixed.effects)
estSE<-c(se.fixed)
est1<-est2<-se1<-se2<-c()
if(length(myidx.cor)>0) {
ipar<-npar
for(iom in 1:nphi) {
for(jom in iom:nphi) {
if(covariance.model[iom,jom]==1) {
ipar<-ipar+1
se.cov[iom,jom]<-se.cov[jom,iom]<-sO[(ipar-npar)]
est1<-c(est1,omega[iom,jom])
se1<-c(se1,sO[ipar-npar])
if(iom==jom) {
se.sdcor[iom,iom]<-sqrt(CO[iom,iom])/2/sqrt(omega[iom,iom])
est2<-c(est2,sqrt(omega[iom,iom]))
se2<-c(se2,se.sdcor[iom,iom])
} else { # compute correlation and SE on correlation using the delta-method
ebet<-c(omega[iom,jom],omega[iom,iom],omega[jom,jom])
varbet<-CO[(myidx.track[myidx.track[,1]==ipar,]-npar),(myidx.track[myidx.track[,1]==ipar,]-npar)]
rho<-ebet[1]/sqrt(ebet[2]*ebet[3])
debet<-c(1/sqrt(ebet[2]*ebet[3]), -ebet[1]/(ebet[2]**(3/2))/sqrt(ebet[3])/2, -ebet[1]/(ebet[3]**(3/2))/sqrt(ebet[2])/2)
se.sdcor[iom,jom]<-se.sdcor[jom,iom]<-t(debet) %*% varbet %*% debet
est2<-c(est2,rho)
se2<-c(se2,se.sdcor[iom,jom])
}
}
}
}
if(saemixObject@model@modeltype=="structural") estpar<-c(estpar,est1,saemixObject@results@respar[saemixObject@results@indx.res],est2) else
estpar<-c(estpar,est1,est2)
if(saemixObject@model@modeltype=="structural") estSE<-c(estSE,se1,se.res[saemixObject@results@indx.res],se2) else estSE<-c(estSE,se1,se2)
} else {
diag(se.cov)<-se.omega
if(saemixObject@model@modeltype=="structural")
estpar<-c(estpar,diag(omega)[saemixObject@results@indx.omega],saemixObject@results@respar[saemixObject@results@indx.res], sqrt(diag(omega)[saemixObject@results@indx.omega])) else estpar<-c(estpar,diag(omega)[saemixObject@results@indx.omega], sqrt(diag(omega)[saemixObject@results@indx.omega]))
if(saemixObject@model@modeltype=="structural")
estSE<-c(estSE,se.omega[saemixObject@results@indx.omega],se.res[saemixObject@results@indx.res],se.omega[saemixObject@results@indx.omega]/2/sqrt(diag(omega)[saemixObject@results@indx.omega])) else estSE<-c(estSE,se.omega[saemixObject@results@indx.omega],se.omega[saemixObject@results@indx.omega]/2/sqrt(diag(omega)[saemixObject@results@indx.omega]))
}
conf.int<-data.frame(name=namallpar, estimate=estpar, se=estSE)
conf.int$cv<-100*conf.int$se/conf.int$estimate
conf.int$lower<-conf.int$estimate - 1.96*conf.int$se
conf.int$upper<-conf.int$estimate + 1.96*conf.int$se
saemix.res["se.fixed"]<-se.fixed
saemix.res["se.omega"]<-c(se.omega)
saemix.res["se.cov"]<-se.cov
if(saemixObject@model@modeltype=="structural") saemix.res["se.respar"]<-c(se.res)
saemix.res["conf.int"]<-conf.int
saemix.res["ll.lin"]<-c(ll.lin )
saemix.res["fim"]<-fim
saemix.res["aic.lin"]<-(-2)*saemix.res["ll.lin"]+ 2*saemix.res["npar.est"]
saemix.res["bic.lin"]<-(-2)*saemix.res["ll.lin"]+ log(saemix.data["N"])*saemix.res["npar.est"]
saemix.res["bic.covariate.lin"]<-(-2)*saemix.res["ll.lin"]+ log(saemix.data["N"])*saemix.res["nbeta.random"]+log(sum(saemix.data["nind.obs"]))*saemix.res["nbeta.fixed"]
##################################
#delete(hw)
saemixObject["results"]<-saemix.res
return(saemixObject)
# return(list(ll.lin,fim,DFi, Dzi, invVi))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.