Estimation of the Fisher Information Matrix

In varTestnlme: Variance Components Testing for Linear and Nonlinear Mixed Effects Models

LOCAL <- identical(Sys.getenv("LOCAL"), "true")
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  purl = LOCAL)

library(varTestnlme)

When testing that the variance of at least one random effect is equal to 0, the limiting distribution of the test statistic is a chi-bar-square distribution whose weights depend on the Fisher Information Matrix (FIM) of the model.

varCompTestnlme provides different ways to handle the FIM. Depending on the package used to fit the model and on your preferences, it is possible to specify different options with the argument fim:

• fim="extract": (the default option) extract the FIM computed by the package used to fit the models. Please note that it might not be available in some cases. For example, for nonlinear models fitted with lme4. It might also rely on some linearization of the model.

• fim="compute": to estimate the FIM using parametric bootstrap. Note that it can take some time to run. In this case, one must provide the bootstrap sample size with the control argument

• fim=I: with I a symmetric positive definite matrix, to provide your own FIM

We provide some examples below.

Fitting the models

We use a dataset on high-flux hemodialyzers, available in the nlme package, in which ultrafiltration rates of 20 dialyzers were measured at 7 different dates. A nonlinear model is considered, with two random effects on the asymptote and on the rate.

$$y_{ij} = A_i(1 - \exp(-e^{l_i}(t_{ij} - c_i))) + \varepsilon_{ij}, \quad \varepsilon_{ij} \sim \mathcal{N}(0,\sigma^2)$$ $$(A_i,l_i,c_i)^t \sim \mathcal{N}(\beta,\Gamma) $$

Suppose we wish to test if the asymptote is indeed random, i.e.: $$H_0: \Gamma = \begin{pmatrix} \gamma_1^2 & 0 & 0\ 0 & 0 & 0 \ 0 & 0 & 0 \end{pmatrix} \quad \text{versus} \quad H_1 : \Gamma = \begin{pmatrix} \gamma_1^2 & 0 & 0 \ 0 & \gamma_2^2 & 0 \ 0 & 0 & \gamma_3^2 \end{pmatrix}$$

Using nlme, the two models under $H_0$ and $H_1$ can be specified in the following way:

library(nlme)
fm1.lis <- nlsList(rate ~ SSasympOff(pressure, Asym, lrc, c0), data=Dialyzer)
m1 <- nlme(fm1.lis, random = pdDiag(Asym + lrc + c0 ~ 1))
m0 <- nlme(fm1.lis, random = pdDiag(Asym ~ 1))

Run the test

To run the test with the default settings, simply use:

varCompTest(m1,m0)

With the default settings, only bounds on the p-value are computed. In this case it is largely enough to reject the null hypothesis that the asymptote is the only random effect. However, if more precision is needed, one should specify it using pval.comp = "both" or pval.comp="approx". By default, the fim option is set to "extract", which means that varCompTestnlme will extract the FIM computed by the package used to fit the model (i.e. nlme, lme4 or saemix).

varCompTest(m1,m0,pval.comp = "both",fim="extract")

In this case we get an error message from nlme package stating that the FIM is non-positive definite. We will then use the fim="compute" option to approximate the FIM using parametric bootstrap, and we set the bootstrap sample size to 100 with the control argument.

varCompTest(m1,m0,pval.comp="both",fim="compute",control=list(B=100))

The weights of the chi-bar-square distribution are provided, along with their standard deviations associated to the Monte Carlo approximation of the weights. Note that exact formulas exist for the weights when the number of chi-bar-square components is less than or equal to 3, in which cases no Monte Carlo approximation is needed and standard deviations are thus null.

varTestnlme documentation built on Sept. 22, 2023, 5:13 p.m.