R/func_FIM.R

In saemixdevelopment/saemix: Stochastic Approximation Expectation Maximization (SAEM) Algorithm

###########################  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 the 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 Comets  E, Lavenu A, Lavielle M. Parameter estimation in nonlinear mixed effect models using saemix, an R implementation of the SAEM algorithm. Journal of Statistical Software 80, 3 (2017), 1-41.
#' 
#' Kuhn E, Lavielle M. Maximum likelihood estimation in nonlinear mixed effects models. Computational Statistics and Data Analysis 49, 4 (2005), 1020-1038.
#' 
#' Comets E, Lavenu A, Lavielle M. SAEMIX, an R version of the SAEM algorithm.
#' 20th meeting of the Population Approach Group in Europe, Athens, Greece
#' (2011), 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)])
  nres<-length(saemix.res["indx.res"])
  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))
  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,]
    Vi<- DF[j1:j2,,drop=FALSE] %*% omega %*% t(DF[j1:j2,,drop=FALSE]) + mydiag((g0[j1:j2])^2, nrow=ni)
    #    invVi[[i]]<-solve(Vi[[i]])
    # Invert avoiding numerical problems
    Gi[[i]]<-round(Vi*1e10)/1e10
    VD<-try(eigen(Gi[[i]]))
    if(is(VD,"try-error")) {
      cat("Unable to compute the FIM by linearisation.\n")
      stop()
      #    return(saemixObject)
    }
    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"]
  j<-which(saemix.model["betaest.model"]>0)
  ind.fixed.est<-(covariate.estim[j]>0)
  npar<-length(ind.fixed.est)
  
  # 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]
    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()
    myidx.omega<-c()
    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
          if(iom==jom) myidx.omega<-c(myidx.omega,ipar)
          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,
          DV[[ipar]]<-DFi %*% domega %*% t(DFi)
        }
      }
    }
    # for(ipar in 1:nomega) {
    #   domega<-omega.null
    #   domega[ipar,ipar]<-sqrt(omega[ipar,ipar])*2
    #   DV[[ipar+npar]] <- DFi %*% t(DFi)
    # }
    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)
      }
    }
    #    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))
#    indMF[[i]]<-MFi
    MF <- MF+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(is(Cth,"try-error")) {
      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(is(CO,"try-error")) {
    CO<-NA*FO
    cat("Error computing the Fisher Information Matrix: singular system.\n")
  }
  sO<-sqrt(mydiag(CO))
  se.omega<-matrix(0,nphi,1)
  se.cov<-matrix(0,nphi,nphi)
  se.omega[saemix.model["indx.omega"]]<-sO[myidx.omega-npar]
  if(length(saemix.model["indx.cov"])>0) {
    ipar<-0
    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]
        }
      }
    }
  } else diag(se.cov)<-se.omega
  se.res<-matrix(0,2*nytype,1)
  se.res[saemix.res["indx.res"]]<-sO[(nomega+1):length(sO)]    
  
  saemix.res["se.fixed"]<-se.fixed
  saemix.res["se.omega"]<-c(se.omega)
  saemix.res["se.cov"]<-se.cov
  saemix.res["se.respar"]<-c(se.res)
  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"]
  
  ##################################
  #delete(hw)
  saemixObject["results"]<-saemix.res
  return(saemixObject)
#  return(list(ll.lin,fim,DFi, Dzi, invVi))
}