## Tutorial on the inference on the parameters of the linear mixed model using the R package nlme

Back to course page

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!

## Data

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.

### Data transformations

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)

##   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]"))

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


## Model

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


## Inference on fixed effects

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).

### Estimates

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


### Covariance matrix of estimates and standard errors

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)  ##  108  print(p <- length(fixef(fit.ML)))  ##  6  print(corr.SE.ML <- sqrt(N/(N - p)))  ##  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  ### Wald tests and confidence intervals using the t- and F-approximations 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 ## 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 ## 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 ## 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") ##  "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 and confidence intervals using the asymptotic normal and $$\chi^2$$-approximations 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 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  ### Inference on variance components 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))

##  0.066


Back to course page