This tutorial explains some of the output from the lme function and related methods. The same dataset is analyzed as in the related tutorial on grouped data analysis. Here, purely technical details are explained and not all other important issues of a full statistical analysis!
library("nlme")
library("lattice")
data(Orthodont, package = "nlme")
This data frame contains the following columns:
distance A numeric vector of distances (mm) from the pituitary (hypophysis gland) to the pterygomaxillary fissure (see the picture below). These distances are measured on x-ray images of the skull.age A numeric vector of ages of the subject (yr).Subject An ordered factor indicating the subject on which the measurement was made. The levels are labelled M01 to M16 for the males and F01 to F13 for the females. The ordering is by increasing average distance within sex.Sex A factor with levels Male and Female.We create few derived variables that will be used later on. Variable Sex3 is an artificial factor variable with levels Male1 (the first four boys in the dataset), Male2 (remaining boys in the dataset) and Female (all girls). It is included here purely technically to have a factor with more than two levels. In practice, such split of boys into two groups would be sensibe if there is a supposition that one group of boys may behave differently from the second group of boys.
Ort <- as.data.frame(Orthodont)
head(Ort)
## distance age Subject Sex
## 1 26.0 8 M01 Male
## 2 25.0 10 M01 Male
## 3 29.0 12 M01 Male
## 4 31.0 14 M01 Male
## 5 21.5 8 M02 Male
## 6 22.5 10 M02 Male
Ort <- transform(Ort, age8 = age - 8,
Sex3 = factor(2*(Sex == "Male" & (Subject %in% c("M01", "M02", "M03", "M04")))
+ 3*(Sex == "Male" & !(Subject %in% c("M01", "M02", "M03", "M04")))
+ 1 * (Sex == "Female"), labels = c("Female", "Male1", "Male2")))
Ort8 <- groupedData(formula = distance ~ age8 | Subject,
data = Ort, outer = ~Sex3,
labels = list(x = "Age", y = "Distance from pituitary to pterygomaxillary fissure"),
units = list(x = "[yr]", y = "[mm]"))
head(Ort8)
## Grouped Data: distance ~ age8 | Subject
## <environment: 0xd1e1e54>
## distance age Subject Sex age8 Sex3
## 1 26.0 8 M01 Male 0 Male1
## 2 25.0 10 M01 Male 2 Male1
## 3 29.0 12 M01 Male 4 Male1
## 4 31.0 14 M01 Male 6 Male1
## 5 21.5 8 M02 Male 0 Male1
## 6 22.5 10 M02 Male 2 Male1
The following model will be fitted: random intercept and random slope in age starting at 8 years with the mean of both the intercept and slope depending on gender.
\[
Y_{ij}=\beta_0 + \beta_1\,I[Male1] + \beta_2\,I[Male2] + (t_j - 8)\,(\beta_3 + \beta_4\,I[Male1] + \beta_5\,I[Male2]) + b_{i0}+(t_j-8)b_{i1}+\varepsilon_{ij},
\]
where \(Y_{ij}\) is the \(j\)-th distance measurement on the \(i\)-th subject, \(t_j\) is the age at the \(j\)-th measurement (the same for all subjects), \(\beta_0\), \(\beta_0 + \beta_1\), and \(\beta_0 + \beta_2\) respectively, is the population mean distance at 8 years for Female, Male1 and Male2, respectively, \(\beta_3\), \(\beta_3 + \beta_4\) and \(\beta_3 + \beta_5\), respectively, is the mean growth in distance in one year for Female, Male1 and Male2, respectively, \(b_{i0}\) is the zero-mean random intercept for the \(i\)-th subject (the deviation of \(Y_{ij}-\varepsilon_{ij}\) from \(\beta_0\) or \(\beta_0 + \beta_1\) or \(\beta_0 + \beta_2\) at the age 8), \(b_{i1}\) is the zero-mean random slope for the \(i\)-th subject (the deviation of the growth rate from \(\beta_3\) or \(\beta_3 + \beta_4\) or \(\beta_3 + \beta_5\)), and \(\varepsilon_{ij}\) are uncorrelated random errors.
We calculate both REML (default) and ML estimates.
fit <- lme(distance ~ Sex3 + age8 + Sex3:age8, random = ~ age8 | Subject, data = Ort8)
fit.ML <- lme(distance ~ Sex3 + age8 + Sex3:age8, random = ~ age8 | Subject, data = Ort8, method = "ML")
The most important results as provided by the summary function:
summary(fit)
## Linear mixed-effects model fit by REML
## Data: Ort8
## AIC BIC logLik
## 450.712 476.9618 -215.356
##
## Random effects:
## Formula: ~age8 | Subject
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.8000018 (Intr)
## age8 0.1848086 -0.051
## Residual 1.3100400
##
## Fixed effects: distance ~ Sex3 + age8 + Sex3:age8
## Value Std.Error DF t-value p-value
## (Intercept) 21.209091 0.6354204 78 33.37805 0.0000
## Sex3Male1 2.303409 1.2304863 24 1.87195 0.0734
## Sex3Male2 1.107576 0.8796997 24 1.25904 0.2201
## age8 0.479545 0.1044311 78 4.59198 0.0000
## Sex3Male1:age8 0.182955 0.2022300 78 0.90469 0.3684
## Sex3Male2:age8 0.345455 0.1445784 78 2.38939 0.0193
## Correlation:
## (Intr) Sx3Ml1 Sx3Ml2 age8 S3M1:8
## Sex3Male1 -0.516
## Sex3Male2 -0.722 0.373
## age8 -0.376 0.194 0.272
## Sex3Male1:age8 0.194 -0.376 -0.140 -0.516
## Sex3Male2:age8 0.272 -0.140 -0.376 -0.722 0.373
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.075968859 -0.385289821 -0.004697224 0.429112299 3.838937259
##
## Number of Observations: 108
## Number of Groups: 27
summary(fit.ML)
## Linear mixed-effects model fit by maximum likelihood
## Data: Ort8
## AIC BIC logLik
## 446.4945 473.3158 -213.2473
##
## Random effects:
## Formula: ~age8 | Subject
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.6572682 (Intr)
## age8 0.1443073 0.056
## Residual 1.3100398
##
## Fixed effects: distance ~ Sex3 + age8 + Sex3:age8
## Value Std.Error DF t-value p-value
## (Intercept) 21.209091 0.6164494 78 34.40524 0.0000
## Sex3Male1 2.303409 1.1937492 24 1.92956 0.0656
## Sex3Male2 1.107576 0.8534356 24 1.29778 0.2067
## age8 0.479545 0.1013130 78 4.73331 0.0000
## Sex3Male1:age8 0.182955 0.1961917 78 0.93253 0.3539
## Sex3Male2:age8 0.345455 0.1402614 78 2.46293 0.0160
## Correlation:
## (Intr) Sx3Ml1 Sx3Ml2 age8 S3M1:8
## Sex3Male1 -0.516
## Sex3Male2 -0.722 0.373
## age8 -0.376 0.194 0.272
## Sex3Male1:age8 0.194 -0.376 -0.140 -0.516
## Sex3Male2:age8 0.272 -0.140 -0.376 -0.722 0.373
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.326883708 -0.423752992 -0.003944058 0.445384713 3.851262361
##
## Number of Observations: 108
## Number of Groups: 27
Issues that are common for models fitted by both REML and ML method will only be shown on the object fit (REML results). Only those issues that make only sense for the ML fits will be shown on the fit.ML object (ML results).
Estimated fixed effects which were also seen in the output from summary(fit):
fixef(fit)
## (Intercept) Sex3Male1 Sex3Male2 age8 Sex3Male1:age8
## 21.2090909 2.3034091 1.1075758 0.4795455 0.1829545
## Sex3Male2:age8
## 0.3454545
Approximate covariance matrix for the fixed effects (as also shown in the lecture):
vcov(fit)
## (Intercept) Sex3Male1 Sex3Male2 age8
## (Intercept) 0.40375909 -0.40375909 -0.40375909 -0.02495666
## Sex3Male1 -0.40375909 1.51409658 0.40375909 0.02495666
## Sex3Male2 -0.40375909 0.40375909 0.77387159 0.02495666
## age8 -0.02495666 0.02495666 0.02495666 0.01090586
## Sex3Male1:age8 0.02495666 -0.09358747 -0.02495666 -0.01090586
## Sex3Male2:age8 0.02495666 -0.02495666 -0.04783360 -0.01090586
## Sex3Male1:age8 Sex3Male2:age8
## (Intercept) 0.02495666 0.02495666
## Sex3Male1 -0.09358747 -0.02495666
## Sex3Male2 -0.02495666 -0.04783360
## age8 -0.01090586 -0.01090586
## Sex3Male1:age8 0.04089698 0.01090586
## Sex3Male2:age8 0.01090586 0.02090290
fit$varFix
## (Intercept) Sex3Male1 Sex3Male2 age8
## (Intercept) 0.40375909 -0.40375909 -0.40375909 -0.02495666
## Sex3Male1 -0.40375909 1.51409658 0.40375909 0.02495666
## Sex3Male2 -0.40375909 0.40375909 0.77387159 0.02495666
## age8 -0.02495666 0.02495666 0.02495666 0.01090586
## Sex3Male1:age8 0.02495666 -0.09358747 -0.02495666 -0.01090586
## Sex3Male2:age8 0.02495666 -0.02495666 -0.04783360 -0.01090586
## Sex3Male1:age8 Sex3Male2:age8
## (Intercept) 0.02495666 0.02495666
## Sex3Male1 -0.09358747 -0.02495666
## Sex3Male2 -0.02495666 -0.04783360
## age8 -0.01090586 -0.01090586
## Sex3Male1:age8 0.04089698 0.01090586
## Sex3Male2:age8 0.01090586 0.02090290
Approximate standard errors for the fixed effects which were also seen in the output from summary(fit):
sqrt(diag(vcov(fit)))
## (Intercept) Sex3Male1 Sex3Male2 age8 Sex3Male1:age8
## 0.6354204 1.2304863 0.8796997 0.1044311 0.2022300
## Sex3Male2:age8
## 0.1445784
In case of ML estimation, standard errors printed by default by summary(fit.ML) are corrected by a factor
\[
\sqrt{\frac{N}{N - p}},
\]
where \(N\) is the total number of observations and \(p\) is the number of fixed effects and hence output from sqrt(diag(vcov(fit.ML))) does not exactly match the standard errors seen in the output from summary(fit.ML). The code below provides the standard errors seen before in the output from summary(fit.ML).
print(N <- fit.ML$dims$N)
## [1] 108
print(p <- length(fixef(fit.ML)))
## [1] 6
print(corr.SE.ML <- sqrt(N/(N - p)))
## [1] 1.028992
corr.SE.ML * sqrt(diag(vcov(fit.ML)))
## (Intercept) Sex3Male1 Sex3Male2 age8 Sex3Male1:age8
## 0.6164494 1.1937492 0.8534356 0.1013130 0.1961917
## Sex3Male2:age8
## 0.1402614
Method summary can also be asked to provide uncorrected standard errors:
sqrt(diag(vcov(fit.ML)))
## (Intercept) Sex3Male1 Sex3Male2 age8 Sex3Male1:age8
## 0.59908118 1.16011571 0.82939033 0.09845849 0.19066405
## Sex3Male2:age8
## 0.13630961
summary(fit.ML, adjustSigma = FALSE)
## Linear mixed-effects model fit by maximum likelihood
## Data: Ort8
## AIC BIC logLik
## 446.4945 473.3158 -213.2473
##
## Random effects:
## Formula: ~age8 | Subject
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.6572682 (Intr)
## age8 0.1443073 0.056
## Residual 1.3100398
##
## Fixed effects: distance ~ Sex3 + age8 + Sex3:age8
## Value Std.Error DF t-value p-value
## (Intercept) 21.209091 0.5990812 78 35.40270 0.0000
## Sex3Male1 2.303409 1.1601157 24 1.98550 0.0586
## Sex3Male2 1.107576 0.8293903 24 1.33541 0.1943
## age8 0.479545 0.0984585 78 4.87053 0.0000
## Sex3Male1:age8 0.182955 0.1906641 78 0.95956 0.3402
## Sex3Male2:age8 0.345455 0.1363096 78 2.53434 0.0133
## Correlation:
## (Intr) Sx3Ml1 Sx3Ml2 age8 S3M1:8
## Sex3Male1 -0.516
## Sex3Male2 -0.722 0.373
## age8 -0.376 0.194 0.272
## Sex3Male1:age8 0.194 -0.376 -0.140 -0.516
## Sex3Male2:age8 0.272 -0.140 -0.376 -0.722 0.373
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.326883708 -0.423752992 -0.003944058 0.445384713 3.851262361
##
## Number of Observations: 108
## Number of Groups: 27
Approximate degrees of freedom for the fixed effects shown in the output from summary(fit):
fit$fixDF$X
## (Intercept) Sex3Male1 Sex3Male2 age8 Sex3Male1:age8
## 78 24 24 78 78
## Sex3Male2:age8
## 78
In case that a factor with more than two levels is included in the model, the approximate degrees of freedom for each term in the model (group of $\beta$s that relate to all reparameterizing (pseudo)contrasts) are obtained as
fit$fixDF$terms
## (Intercept) Sex3 age8 Sex3:age8
## 78 24 78 78
Approximate t-tests on fixed effects which were also seen in the output from summary(fit):
tstat <- fixef(fit)/sqrt(diag(vcov(fit)))
pval <- 2*pt(-abs(tstat), df = fit$fixDF$X)
Tresult <- data.frame(t = tstat, p = round(pval, 4))
print(Tresult)
## t p
## (Intercept) 33.3780451 0.0000
## Sex3Male1 1.8719502 0.0734
## Sex3Male2 1.2590384 0.2201
## age8 4.5919778 0.0000
## Sex3Male1:age8 0.9046854 0.3684
## Sex3Male2:age8 2.3893932 0.0193
Also function anova is supplied by the respective L argument provides above t-tests (the test statistic is reported as F-statistic, i.e., is equal to the squared t-statistic):
anova(fit, L = c(0, 1, 0, 0, 0, 0)) ## Sex3Male1
## F-test for linear combination(s)
## [1] 1
## numDF denDF F-value p-value
## 1 1 24 3.504197 0.0734
anova(fit, L = c(0, 0, 1, 0, 0, 0)) ## Sex3Male2
## F-test for linear combination(s)
## [1] 1
## numDF denDF F-value p-value
## 1 1 24 1.585178 0.2201
anova(fit, L = c(0, 0, 0, 1, 0, 0)) ## age8
## F-test for linear combination(s)
## [1] 1
## numDF denDF F-value p-value
## 1 1 78 21.08626 <.0001
Approximate confidence intervals for the fixed effects (based on the above t-statistics) are calculated by the intervals function:
intervals(fit)
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 19.94406606 21.2090909 22.4741158
## Sex3Male1 -0.23618985 2.3034091 4.8430080
## Sex3Male2 -0.70803522 1.1075758 2.9231867
## age8 0.27163904 0.4795455 0.6874519
## Sex3Male1:age8 -0.21965450 0.1829545 0.5855636
## Sex3Male2:age8 0.05762114 0.3454545 0.6332879
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: Subject
## lower est. upper
## sd((Intercept)) 1.20349767 1.80000182 2.6921586
## sd(age8) 0.05986179 0.18480865 0.5705515
## cor((Intercept),age8) -0.79122253 -0.05138195 0.7495108
##
## Within-group standard error:
## lower est. upper
## 1.083906 1.310040 1.583352
And now by hand:
print(tquants <- qt(0.975, df = fit$fixDF$X))
## (Intercept) Sex3Male1 Sex3Male2 age8 Sex3Male1:age8
## 1.990847 2.063899 2.063899 1.990847 1.990847
## Sex3Male2:age8
## 1.990847
Low <- fixef(fit) - tquants * sqrt(diag(vcov(fit)))
Upp <- fixef(fit) + tquants * sqrt(diag(vcov(fit)))
EstInt <- data.frame(Lower = Low, Estimate = fixef(fit), Upper = Upp)
print(EstInt)
## Lower Estimate Upper
## (Intercept) 19.94406606 21.2090909 22.4741158
## Sex3Male1 -0.23618985 2.3034091 4.8430080
## Sex3Male2 -0.70803522 1.1075758 2.9231867
## age8 0.27163904 0.4795455 0.6874519
## Sex3Male1:age8 -0.21965450 0.1829545 0.5855636
## Sex3Male2:age8 0.05762114 0.3454545 0.6332879
Approximate F-tests on fixed effects but corresponding to type I (sequential) tests:
anova(fit)
## numDF denDF F-value p-value
## (Intercept) 1 78 3928.133 <.0001
## Sex3 2 24 4.051 0.0305
## age8 1 78 98.094 <.0001
## Sex3:age8 2 78 2.855 0.0636
And now by hand for the Sex3:age8 interaction:
efNames <- c("Sex3Male1:age8", "Sex3Male2:age8")
beta <- fixef(fit)[efNames]
Vbeta <- vcov(fit)[efNames, efNames]
numDF <- length(beta)
denDF <- fit$fixDF$terms["Sex3:age8"]
F <- as.numeric(crossprod(beta, solve(Vbeta, beta))) / numDF
pval <- pf(F, df1 = numDF, df2 = denDF, lower.tail = FALSE)
Fresult <- data.frame(F = F, p = round(pval, 4), row.names = "Sex3:age8")
print(Fresult)
## F p
## Sex3:age8 2.854705 0.0636
The Wald test based on F-approximation for arbitrary linear combination(s) of the fixed effects is obtained by invoking the anova command with a relevant value of its L argument. For example, the above F-test on Sex3:age8 interaction is obtained as follows:
anova(fit, L = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 2, byrow = TRUE))
## F-test for linear combination(s)
## Sex3Male1:age8 Sex3Male2:age8
## 1 1 0
## 2 0 1
## numDF denDF F-value p-value
## 1 2 78 2.854705 0.0636
Similarly, the F-test for the main effect of Sex3 (what is it the interpretation of the test?) is obtained using
anova(fit, L = matrix(c(0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), nrow = 2, byrow = TRUE))
## F-test for linear combination(s)
## Sex3Male1 Sex3Male2
## 1 1 0
## 2 0 1
## numDF denDF F-value p-value
## 1 2 24 1.934759 0.1663
Once again, note that output on row Sex3 in anova(fit) is telling something else (because it is doing something else):
anova(fit)
## numDF denDF F-value p-value
## (Intercept) 1 78 3928.133 <.0001
## Sex3 2 24 4.051 0.0305
## age8 1 78 98.094 <.0001
## Sex3:age8 2 78 2.855 0.0636
Wald tests on fixed effects based on asymptotic normal approximations (the test statistic are indeed the same as for the Wald tests based on the t-approximation):
tstat <- fixef(fit)/sqrt(diag(vcov(fit)))
pval <- 2*pnorm(-abs(tstat))
Wresult <- data.frame(Z = tstat, p = round(pval, 4))
print(Wresult)
## Z p
## (Intercept) 33.3780451 0.0000
## Sex3Male1 1.8719502 0.0612
## Sex3Male2 1.2590384 0.2080
## age8 4.5919778 0.0000
## Sex3Male1:age8 0.9046854 0.3656
## Sex3Male2:age8 2.3893932 0.0169
Wald test on a group of fixed effects based on the asymptotic \(\chi^2\) approximation will be shown on a test on the significance of the Sex3:age8 interaction.
efNames <- c("Sex3Male1:age8", "Sex3Male2:age8")
beta <- fixef(fit)[efNames]
Vbeta <- vcov(fit)[efNames, efNames]
Q <- as.numeric(crossprod(beta, solve(Vbeta, beta)))
pval <- pchisq(Q, df = length(beta), lower.tail = FALSE)
Waldresult <- data.frame(Q = Q, p = round(pval, 4), row.names = "Sex3:age8")
print(Waldresult)
## Q p
## Sex3:age8 5.70941 0.0576
Likelihood ratio tests on fixed effects can only be performed if a method of maximum-likelihood (ML) is used to fit the model. Below, we show the likelihood-ratio test to evaluate the siginificance of the Sex3:age8 interaction:
fit0.ML <- lme(distance ~ Sex3 + age8, random = ~ age8 | Subject, data = Ort8, method = "ML")
anova(fit0.ML, fit.ML)
## Model df AIC BIC logLik Test L.Ratio p-value
## fit0.ML 1 8 448.2566 469.7136 -216.1283
## fit.ML 2 10 446.4945 473.3158 -213.2473 1 vs 2 5.762069 0.0561
And now by hand:
logLik <- logLik(fit.ML)
logLik0 <- logLik(fit0.ML)
LR <- 2 * (logLik - logLik0)
pval <- pchisq(LR, df = 2, lower.tail = FALSE)
LRresult <- data.frame(LR = LR, p = round(pval, 4), row.names = "Sex3:age8")
print(LRresult)
## LR p
## Sex3:age8 5.762069 0.0561
Let us now assume three nested models: (1) random intercept and random slope model with general positive definite covariance matrix for random effects, (2) random intercept and random slope model with a diagonal covariance matrix for random effects, (3) random intercept, fixed slope model. In all models, no other covariates are considered.
fit.1 <- lme(distance ~ age8, random = pdLogChol(~ age8), data = Ort8)
fit.1.uncor <- lme(distance ~ age8, random = pdDiag(~ age8), data = Ort8)
fit.1.fixSlp <- lme(distance ~ age8, random = pdLogChol(~ 1), data = Ort8)
summary(fit.1)
## Linear mixed-effects model fit by REML
## Data: Ort8
## AIC BIC logLik
## 454.6367 470.6173 -221.3183
##
## Random effects:
## Formula: ~age8 | Subject
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.886629 (Intr)
## age8 0.226428 0.209
## Residual 1.310039
##
## Fixed effects: distance ~ age8
## Value Std.Error DF t-value p-value
## (Intercept) 22.042593 0.4199079 80 52.49387 0
## age8 0.660185 0.0712533 80 9.26533 0
## Correlation:
## (Intr)
## age8 -0.208
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.223105126 -0.493761342 0.007316612 0.472151020 3.916033672
##
## Number of Observations: 108
## Number of Groups: 27
summary(fit.1.uncor)
## Linear mixed-effects model fit by REML
## Data: Ort8
## AIC BIC logLik
## 452.8916 466.2088 -221.4458
##
## Random effects:
## Formula: ~age8 | Subject
## Structure: Diagonal
## (Intercept) age8 Residual
## StdDev: 1.962957 0.2469767 1.293767
##
## Fixed effects: distance ~ age8
## Value Std.Error DF t-value p-value
## (Intercept) 22.042593 0.4314009 80 51.09538 0
## age8 0.660185 0.0732042 80 9.01841 0
## Correlation:
## (Intr)
## age8 -0.294
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.98617508 -0.49418192 0.01618754 0.42884417 3.96519081
##
## Number of Observations: 108
## Number of Groups: 27
summary(fit.1.fixSlp)
## Linear mixed-effects model fit by REML
## Data: Ort8
## AIC BIC logLik
## 455.0025 465.6563 -223.5013
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 2.114724 1.431592
##
## Fixed effects: distance ~ age8
## Value Std.Error DF t-value p-value
## (Intercept) 22.042593 0.4677240 80 47.12735 0
## age8 0.660185 0.0616059 80 10.71626 0
## Correlation:
## (Intr)
## age8 -0.395
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.66453932 -0.53507984 -0.01289591 0.48742859 3.72178465
##
## Number of Observations: 108
## Number of Groups: 27
We will now calculate likelihood-ratio tests comparing the above models. Note that models being compared must be fitted using the same method (REML or ML). But still, classical asyptotic theory in many situations fails due to the fact that under the null hypothesis (smaller model), parameters of the bigger model reach a boundary of the parameter space. Classical likelihood-ratio test is then conservative (giving larger p-values than the p-values it should provide). Hence be careful when deciding on whether a simpler model is really suitable!
anova(fit.1, fit.1.uncor)
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.1 1 6 454.6367 470.6173 -221.3183
## fit.1.uncor 2 5 452.8916 466.2088 -221.4458 1 vs 2 0.2548889 0.6137
anova(fit.1.uncor, fit.1.fixSlp)
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.1.uncor 1 5 452.8916 466.2088 -221.4458
## fit.1.fixSlp 2 4 455.0025 465.6563 -223.5013 1 vs 2 4.110941 0.0426
anova(fit.1, fit.1.fixSlp)
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.1 1 6 454.6367 470.6173 -221.3183
## fit.1.fixSlp 2 4 455.0025 465.6563 -223.5013 1 vs 2 4.36583 0.1127
Nevertheless, a bootstrap method (see NMST434 Modern Statistical Methods class) can be used to calculate a correct p-value. To this end, function simulate.lme is useful to simulate from a model under the null hypothesis and then to fit both null and alternative model to generated data. The code below provides the bootstrap simulation while generating data according to the estimated parameters of model fit.1.fixSlp (estimates based on our data) and then fitting both models fit.1 and fit.1.fixSlp to generated data.
set.seed(221913282)
msim <- simulate.lme(fit.1.fixSlp, m2 = fit.1, nsim = 1000, method = "REML")
Maximized log-likelihoods under the null model (first 10 replicated datasets):
print(msim$null$REML[1:10,])
## info logLik
## 1 0 -224.6887
## 2 0 -234.0964
## 3 0 -221.3357
## 4 0 -227.8287
## 5 0 -215.2709
## 6 0 -226.2425
## 7 0 -218.3154
## 8 0 -215.2921
## 9 0 -222.4872
## 10 0 -232.3065
Maximized log-likelihoods under the alternative model (first 10 replicated datasets):
print(msim$alt$REML[1:10,])
## info logLik
## 1 0 -224.2911
## 2 0 -233.1589
## 3 0 -221.3393
## 4 0 -227.4392
## 5 0 -212.0168
## 6 0 -226.2454
## 7 0 -217.1792
## 8 0 -215.2995
## 9 0 -221.4039
## 10 0 -231.8556
Bootstrap sample of the values of the likelihood-ratio test statistic. Some values are slightly negative. This is likely caused by not enough precise convergence of an algorithm used to fit the models. We change negative values to zeros.
LRboot <- -2 * (msim$null$REML[, "logLik"] - msim$alt$REML[, "logLik"])
summary(LRboot)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.01709 0.22403 0.87788 1.40596 1.92214 12.60870
summary(LRboot[LRboot < 0])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.017091 -0.011050 -0.007684 -0.007861 -0.003432 -0.000352
LRboot[LRboot < 0] <- 0
Histogram of the values of the likelihood-ratio statistics (bootstrap estimate of the distribution of the likelihood-ratio statistic under the null hypothesis) compared to assumed asymptotic \(\chi^2_2\) density:
x <- seq(0, 13, length = 500)
dx <- dchisq(x, df = 2)
px <- pchisq(x, df = 2)
#
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1) + 0.1)
hist(LRboot, prob = TRUE, main = "", xlab = "Likelihood-ratio statistic", ylab = "Density", col = "palegreen")
lines(x, dx, col = "red3", lwd = 2)
#
plot(ecdf(LRboot), pch = NA, col = "darkgreen", lwd = 2, xlab = "Likelihood-ratio statistic", ylab = "CDF", main = "")
lines(x, px, col = "red3", lwd = 2)
legend(4, 0.3, legend = c("True", "Chi^2_2"), col = c("darkgreen", "red3"), lwd = 2, lty = 1)
par(mfrow = c(1, 1), mar = c(5, 4, 4, 1) + 0.1)
Can you see that the true distribution of the likelihood-ratio test statistic is stochastically smaller than the usual \(\chi^2_2\) distribution?
And finally, original (conservative) p-value and a bootstrap based p-value of the likelihood ratio test:
anova(fit.1, fit.1.fixSlp)
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.1 1 6 454.6367 470.6173 -221.3183
## fit.1.fixSlp 2 4 455.0025 465.6563 -223.5013 1 vs 2 4.36583 0.1127
LRdata <- anova(fit.1, fit.1.fixSlp)[2, "L.Ratio"]
(p.boot <- sum(LRboot > LRdata) / length(LRboot))
## [1] 0.066