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.
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:
To run the test with the default settings, simply use:
varCompTest(m1,m0)
#> Variance components testing in mixed effects models
#> Testing that the covariance matrix of lrc and c0 is equal to 0
#> Likelihood ratio test statistic:
#> LRT = 28.98229
#> bounds on p-value: lower 3.652152e-08 upper 2.909383e-07
#> 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")
#> Variance components testing in mixed effects models
#> Testing that the covariance matrix of lrc and c0 is equal to 0
#> Likelihood ratio test statistic:
#> LRT = 28.98229
#>
#> p-value from estimated weights: 1.647692e-07
#> bounds on p-value: lower 3.652152e-08 upper 2.909383e-07
#> 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))
#> Variance components testing in mixed effects models
#> Testing that the covariance matrix of lrc and c0 is equal to 0
#> Computing Fisher Information Matrix by bootstrap...
#> ...generating the B=100 bootstrap samples ...
#>
|
| | 0%
|
|. | 1%
|
|. | 2%
|
|.. | 3%
|
|... | 4%
|
|.... | 5%
|
|.... | 6%
|
|..... | 7%
|
|...... | 8%
|
|...... | 9%
|
|....... | 10%
|
|........ | 11%
|
|........ | 12%
|
|......... | 13%
|
|.......... | 14%
|
|........... | 15%
|
|........... | 16%
|
|............ | 17%
|
|............. | 18%
|
|............. | 19%
|
|.............. | 20%
|
|............... | 21%
|
|................ | 22%
|
|................ | 23%
|
|................. | 24%
|
|.................. | 25%
|
|.................. | 26%
|
|................... | 27%
|
|.................... | 28%
|
|..................... | 29%
|
|..................... | 30%
|
|...................... | 31%
|
|....................... | 32%
|
|....................... | 33%
|
|........................ | 34%
|
|......................... | 35%
|
|......................... | 36%
|
|.......................... | 37%
|
|........................... | 38%
|
|............................ | 39%
|
|............................ | 40%
|
|............................. | 41%
|
|.............................. | 42%
|
|.............................. | 43%
|
|............................... | 44%
|
|................................ | 45%
|
|................................. | 46%
|
|................................. | 47%
|
|.................................. | 48%
|
|................................... | 49%
|
|................................... | 51%
|
|.................................... | 52%
|
|..................................... | 53%
|
|..................................... | 54%
|
|...................................... | 55%
|
|....................................... | 56%
|
|........................................ | 57%
|
|........................................ | 58%
|
|......................................... | 59%
|
|.......................................... | 60%
|
|.......................................... | 61%
|
|........................................... | 62%
|
|............................................ | 63%
|
|............................................. | 64%
|
|............................................. | 65%
|
|.............................................. | 66%
|
|............................................... | 67%
|
|............................................... | 68%
|
|................................................ | 69%
|
|................................................. | 70%
|
|................................................. | 71%
|
|.................................................. | 72%
|
|................................................... | 73%
|
|.................................................... | 74%
|
|.................................................... | 75%
|
|..................................................... | 76%
|
|...................................................... | 77%
|
|...................................................... | 78%
|
|....................................................... | 79%
|
|........................................................ | 80%
|
|......................................................... | 81%
|
|......................................................... | 82%
|
|.......................................................... | 83%
|
|........................................................... | 84%
|
|........................................................... | 85%
|
|............................................................ | 86%
|
|............................................................. | 87%
|
|.............................................................. | 88%
|
|.............................................................. | 89%
|
|............................................................... | 90%
|
|................................................................ | 91%
|
|................................................................ | 92%
|
|................................................................. | 93%
|
|.................................................................. | 94%
|
|.................................................................. | 95%
|
|................................................................... | 96%
|
|.................................................................... | 97%
|
|..................................................................... | 98%
|
|..................................................................... | 99%
|
|......................................................................| 100%
#> Likelihood ratio test statistic:
#> LRT = 28.98229
#>
#> p-value from estimated weights: 1.831489e-07
#> bounds on p-value: lower 3.652152e-08 upper 2.909383e-07
#> 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.