# R/func_FIM.R In saemix: Stochastic Approximation Expectation Maximization (SAEM) Algorithm

#### Documented in fim.saemix

###########################  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))
}


## Try the saemix package in your browser

Any scripts or data that you put into this service are public.

saemix documentation built on Aug. 5, 2022, 5:25 p.m.