## Load dataset
dental <- read.table("http://isites.harvard.edu/fs/docs/icb.topic1100055.files/dental.txt",
col.names = c("id","gender","age8","age10","age12","age14"))
dental
id gender age8 age10 age12 age14
1 1 F 21.0 20.0 21.5 23.0
2 2 F 21.0 21.5 24.0 25.5
3 3 F 20.5 24.0 24.5 26.0
4 4 F 23.5 24.5 25.0 26.5
5 5 F 21.5 23.0 22.5 23.5
6 6 F 20.0 21.0 21.0 22.5
7 7 F 21.5 22.5 23.0 25.0
8 8 F 23.0 23.0 23.5 24.0
9 9 F 20.0 21.0 22.0 21.5
10 10 F 16.5 19.0 19.0 19.5
11 11 F 24.5 25.0 28.0 28.0
12 12 M 26.0 25.0 29.0 31.0
13 13 M 21.5 22.5 23.0 26.5
14 14 M 23.0 22.5 24.0 27.5
15 15 M 25.5 27.5 26.5 27.0
16 16 M 20.0 23.5 22.5 26.0
17 17 M 24.5 25.5 27.0 28.5
18 18 M 22.0 22.0 24.5 26.5
19 19 M 24.0 21.5 24.5 25.5
20 20 M 23.0 20.5 31.0 26.0
21 21 M 27.5 28.0 31.0 31.5
22 22 M 23.0 23.0 23.5 25.0
23 23 M 21.5 23.5 24.0 28.0
24 24 M 17.0 24.5 26.0 29.5
25 25 M 22.5 25.5 25.5 26.0
26 26 M 23.0 24.5 26.0 30.0
27 27 M 22.0 21.5 23.5 25.0
## Load reshape2
library(reshape2)
dental.long <- melt(dental, id.vars = c("id","gender"),
value.name = "distance", variable.name = "age")
## delete age from age variable values
dental.long$age <- as.numeric(gsub("age", "", dental.long$age))
dental.long$id <- factor(dental.long$id)
## Reordering
library(doBy)
dental.long <- orderBy( ~ id + age, dental.long)
rownames(dental.long) <- NULL
dental.long
id gender age distance
1 1 F 8 21.0
2 1 F 10 20.0
3 1 F 12 21.5
4 1 F 14 23.0
5 2 F 8 21.0
6 2 F 10 21.5
7 2 F 12 24.0
8 2 F 14 25.5
9 3 F 8 20.5
10 3 F 10 24.0
11 3 F 12 24.5
12 3 F 14 26.0
13 4 F 8 23.5
14 4 F 10 24.5
15 4 F 12 25.0
16 4 F 14 26.5
17 5 F 8 21.5
18 5 F 10 23.0
19 5 F 12 22.5
20 5 F 14 23.5
21 6 F 8 20.0
22 6 F 10 21.0
23 6 F 12 21.0
24 6 F 14 22.5
25 7 F 8 21.5
26 7 F 10 22.5
27 7 F 12 23.0
28 7 F 14 25.0
29 8 F 8 23.0
30 8 F 10 23.0
31 8 F 12 23.5
32 8 F 14 24.0
33 9 F 8 20.0
34 9 F 10 21.0
35 9 F 12 22.0
36 9 F 14 21.5
37 10 F 8 16.5
38 10 F 10 19.0
39 10 F 12 19.0
40 10 F 14 19.5
41 11 F 8 24.5
42 11 F 10 25.0
43 11 F 12 28.0
44 11 F 14 28.0
45 12 M 8 26.0
46 12 M 10 25.0
47 12 M 12 29.0
48 12 M 14 31.0
49 13 M 8 21.5
50 13 M 10 22.5
51 13 M 12 23.0
52 13 M 14 26.5
53 14 M 8 23.0
54 14 M 10 22.5
55 14 M 12 24.0
56 14 M 14 27.5
57 15 M 8 25.5
58 15 M 10 27.5
59 15 M 12 26.5
60 15 M 14 27.0
61 16 M 8 20.0
62 16 M 10 23.5
63 16 M 12 22.5
64 16 M 14 26.0
65 17 M 8 24.5
66 17 M 10 25.5
67 17 M 12 27.0
68 17 M 14 28.5
69 18 M 8 22.0
70 18 M 10 22.0
71 18 M 12 24.5
72 18 M 14 26.5
73 19 M 8 24.0
74 19 M 10 21.5
75 19 M 12 24.5
76 19 M 14 25.5
77 20 M 8 23.0
78 20 M 10 20.5
79 20 M 12 31.0
80 20 M 14 26.0
81 21 M 8 27.5
82 21 M 10 28.0
83 21 M 12 31.0
84 21 M 14 31.5
85 22 M 8 23.0
86 22 M 10 23.0
87 22 M 12 23.5
88 22 M 14 25.0
89 23 M 8 21.5
90 23 M 10 23.5
91 23 M 12 24.0
92 23 M 14 28.0
93 24 M 8 17.0
94 24 M 10 24.5
95 24 M 12 26.0
96 24 M 14 29.5
97 25 M 8 22.5
98 25 M 10 25.5
99 25 M 12 25.5
100 25 M 14 26.0
101 26 M 8 23.0
102 26 M 10 24.5
103 26 M 12 26.0
104 26 M 14 30.0
105 27 M 8 22.0
106 27 M 10 21.5
107 27 M 12 23.5
108 27 M 14 25.0
## Load ggplot2
library(ggplot2)
## plot
plot1 <- ggplot(data = dental.long, mapping = aes(x = age, y = distance, color = id)) +
layer(geom = "line") +
facet_grid(. ~ gender) +
theme_bw() +
theme(legend.key = element_blank())
## Directly label
library(directlabels)
direct.label(plot1)
## Load plyr package
library(plyr)
dental.summary <- ddply(dental.long,
"age",
summarize,
mean = mean(distance),
sd = sd(distance))
dental.summary
age mean sd
1 8 22.18519 2.434322
2 10 23.16667 2.157277
3 12 24.64815 2.817578
4 14 26.09259 2.766687
## Plot mean and +/- 1 SD
plot2 <- ggplot(data = dental.long, mapping = aes(x = age, y = distance)) +
layer(geom = "line", mapping = aes(color = id), alpha = 1/2) +
layer(geom = "line", mapping = aes(y = mean), size = 2, data = dental.summary) +
layer(geom = "errorbar", mapping = aes(ymax = mean + sd, ymin = mean - sd, y = NULL),
width = 0.5, data = dental.summary) +
facet_grid(. ~ gender) +
theme_bw() +
theme(legend.key = element_blank())
## Directly label
direct.label(plot2)
## Load plyr
library(plyr)
## Fit individual models
res.stage1 <- ddply(.data = dental.long,
.variables = c("id","gender"),
.fun = function(DF) {
res <- lm(distance ~ age, data = DF)
coef(res)
})
res.stage1
id gender (Intercept) age
1 1 F 17.25 0.375
2 2 F 14.20 0.800
3 3 F 14.40 0.850
4 4 F 19.65 0.475
5 5 F 19.60 0.275
6 6 F 17.00 0.375
7 7 F 16.95 0.550
8 8 F 21.45 0.175
9 9 F 18.10 0.275
10 10 F 13.55 0.450
11 11 F 18.95 0.675
12 12 M 17.30 0.950
13 13 M 14.85 0.775
14 14 M 16.00 0.750
15 15 M 24.70 0.175
16 16 M 13.65 0.850
17 17 M 18.95 0.675
18 18 M 14.95 0.800
19 19 M 19.75 0.375
20 20 M 14.40 0.975
21 21 M 21.25 0.750
22 22 M 20.05 0.325
23 23 M 13.25 1.000
24 24 M 2.80 1.950
25 25 M 19.10 0.525
26 26 M 13.50 1.125
27 27 M 16.95 0.550
## Means by gender
res.means <- ddply(.data = res.stage1,
.variables = c("gender"),
.fun = function(DF) {
colMeans(DF[,3:4])
})
res.means
gender (Intercept) age
1 F 17.37273 0.4795455
2 M 16.34063 0.7843750
## Intercept model
intercept.model <- lm(`(Intercept)` ~ gender, res.stage1)
summary(intercept.model)
Call:
lm(formula = `(Intercept)` ~ gender, data = res.stage1)
Residuals:
Min 1Q Median 3Q Max
-13.5406 -2.3156 -0.1227 2.4433 8.3594
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.373 1.228 14.143 1.97e-13 ***
genderM -1.032 1.596 -0.647 0.524
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.074 on 25 degrees of freedom
Multiple R-squared: 0.01646, Adjusted R-squared: -0.02288
F-statistic: 0.4183 on 1 and 25 DF, p-value: 0.5237
## Slope model
slope.model <- lm(age ~ gender, res.stage1)
summary(slope.model)
Call:
lm(formula = age ~ gender, data = res.stage1)
Residuals:
Min 1Q Median 3Q Max
-0.60937 -0.20455 -0.02955 0.17812 1.16563
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4795 0.1037 4.623 0.0000989 ***
genderM 0.3048 0.1347 2.262 0.0326 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.344 on 25 degrees of freedom
Multiple R-squared: 0.1699, Adjusted R-squared: 0.1367
F-statistic: 5.119 on 1 and 25 DF, p-value: 0.03261
## SD of residuals from slope model
sd(residuals(slope.model))
[1] 0.3373178
## Load nlme
library(nlme)
## Fit a random intercept and slope model
res.lme <- lme(fixed = distance ~ 1 + gender + age + gender:age,
data = dental.long,
random = ~ 1 + age| id,
method = "REML"
)
## Results
## logLik compatible with SAS. AIC by different method.
## Random effects part: compatible with SAS although shown in SD and Corr rather than Var and Cov
## Fixed effects part: DFs different, thus different t-values, and p-values
summary(res.lme)
Linear mixed-effects model fit by REML
Data: dental.long
AIC BIC logLik
448.5817 469.7368 -216.2908
Random effects:
Formula: ~1 + age | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 2.4055009 (Intr)
age 0.1803455 -0.668
Residual 1.3100396
Fixed effects: distance ~ 1 + gender + age + gender:age
Value Std.Error DF t-value p-value
(Intercept) 17.372727 1.2283958 79 14.142614 0.0000
genderM -1.032102 1.5957329 25 -0.646789 0.5237
age 0.479545 0.1037193 79 4.623492 0.0000
genderM:age 0.304830 0.1347353 79 2.262432 0.0264
Correlation:
(Intr) gendrM age
genderM -0.770
age -0.880 0.678
genderM:age 0.678 -0.880 -0.770
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.168077730 -0.385939008 0.007104087 0.445154545 3.849463576
Number of Observations: 108
Number of Groups: 27
## Interval estimates
intervals(res.lme)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 14.92766665 17.3727273 19.8177879
genderM -4.31857573 -1.0321023 2.2543712
age 0.27309731 0.4795455 0.6859936
genderM:age 0.03664554 0.3048295 0.5730136
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: id
lower est. upper
sd((Intercept)) 1.00632033 2.4055009 5.7500923
sd(age) 0.05846112 0.1803455 0.5563439
cor((Intercept),age) -0.96063873 -0.6676191 0.3285937
Within-group standard error:
lower est. upper
1.084766 1.310040 1.582095
## Estimated G matrix for variances and covariances of random effects (SAS-compatible)
getVarCov(res.lme)
Random effects variance covariance matrix
(Intercept) age
(Intercept) 5.78640 -0.289630
age -0.28963 0.032524
Standard Deviations: 2.4055 0.18035
## Estimated G correlation matrix (SAS-compatible)
cov2cor(getVarCov(res.lme))
Random effects variance covariance matrix
(Intercept) age
(Intercept) 1.00000 -0.66762
age -0.66762 1.00000
Standard Deviations: 1 1
## fixed effect
fixef(res.lme)
(Intercept) genderM age genderM:age
17.3727273 -1.0321023 0.4795455 0.3048295
## random effect
ranef(res.lme)
(Intercept) age
1 -0.64132024 -0.044754845
2 -0.66020223 0.090293750
3 -0.24892690 0.113565208
4 1.66111350 0.028212826
5 0.57096833 -0.054963722
6 -0.82630562 -0.048057940
7 0.05820188 0.023482879
8 1.41328613 -0.071778787
9 -0.53894398 -0.074782288
10 -2.98417340 -0.062697170
11 2.19630252 0.101480090
12 1.58201968 0.081009128
13 -1.15234167 -0.023562635
14 -0.43305242 -0.018682891
15 2.97663820 -0.140968498
16 -1.64534099 -0.008474015
17 1.35484459 -0.010649850
18 -0.94670401 -0.011926906
19 0.36707568 -0.123853840
20 -0.43216727 0.053007723
21 3.45164067 0.050682093
22 0.32577111 -0.140519109
23 -1.15145653 0.048127980
24 -3.88139216 0.302009291
25 0.67597475 -0.070554939
26 -0.30825358 0.103003530
27 -0.78325605 -0.088647061
## ## V matrix: Cov(Yi) = Zi*G*Zi’ + sigma^2*Ini
## G <- as.matrix(getVarCov(res.lme))
## Z <- as.matrix(ranef(res.lme))
## i <- 1
## Z[i,] %*% G %*% t(Z[i,]) +
## Load lme4
library(lme4)
## Fit a random intercept and slope model
res.lmer <- lmer(formula = distance ~ 1 + gender + age + gender:age + (1 + age | id),
data = dental.long,
REML = TRUE
)
## Result
res.lmer
Linear mixed model fit by REML
Formula: distance ~ 1 + gender + age + gender:age + (1 + age | id)
Data: dental.long
AIC BIC logLik deviance REMLdev
448.6 470 -216.3 427.9 432.6
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 5.786427 2.40550
age 0.032524 0.18035 -0.668
Residual 1.716205 1.31004
Number of obs: 108, groups: id, 27
Fixed effects:
Estimate Std. Error t value
(Intercept) 17.3727 1.2284 14.143
genderM -1.0321 1.5957 -0.647
age 0.4795 0.1037 4.623
genderM:age 0.3048 0.1347 2.262
Correlation of Fixed Effects:
(Intr) gendrM age
genderM -0.770
age -0.880 0.678
genderM:age 0.678 -0.880 -0.770
## Fit statistics
logLik(res.lmer)
'log Lik.' -216.2908 (df=8)
AIC(res.lmer)
[1] 448.5817
BIC(res.lmer)
[1] 470.0387
## ANOVA table
anova(res.lmer)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
gender 1 13.758 13.758 8.0164
age 1 171.189 171.189 99.7487
gender:age 1 8.790 8.790 5.1219
## Estimated G matrix and G correlation matrix. SD of residuals (sc). Compatible with SAS
## VarCorr(res.lmer)$id[1:2,1:2] for G matrix only
VarCorr(res.lmer)
$id
(Intercept) age
(Intercept) 5.7864268 -0.28962636
age -0.2896264 0.03252433
attr(,"stddev")
(Intercept) age
2.405499 0.180345
attr(,"correlation")
(Intercept) age
(Intercept) 1.000000 -0.667619
age -0.667619 1.000000
attr(,"sc")
[1] 1.31004
## Extract residual standard error
sigma(res.lmer)
[1] 1.31004
## Fixed effects (population average effects in linear mixed models)
fixef(res.lmer)
(Intercept) genderM age genderM:age
17.3727273 -1.0321023 0.4795455 0.3048295
## Random effects
ranef(res.lmer)
$id
(Intercept) age
1 -0.64126894 -0.044761999
2 -0.66027220 0.090296261
3 -0.24902488 0.113571913
4 1.66105960 0.028225002
5 0.57100802 -0.054964218
6 -0.82624817 -0.048066461
7 0.05817974 0.023484934
8 1.41332651 -0.071775363
9 -0.53886740 -0.074790991
10 -2.98406574 -0.062719915
11 2.19617346 0.101500837
12 1.58191960 0.081024646
13 -1.15230069 -0.023571371
14 -0.43302816 -0.018686886
15 2.97671405 -0.140960547
16 -1.64530512 -0.008484666
17 1.35483093 -0.010642356
18 -0.94667703 -0.011933545
19 0.36718090 -0.123860570
20 -0.43220758 0.053008926
21 3.45153581 0.050706819
22 0.32589205 -0.140527297
23 -1.15148011 0.048124441
24 -3.88159749 0.302007482
25 0.67602668 -0.070555925
26 -0.30834104 0.103009108
27 -0.78316279 -0.088658258
## Individual effects (fixed + random effects)
coef(res.lmer)
$id
(Intercept) genderM age genderM:age
1 16.73146 -1.032102 0.4347835 0.3048295
2 16.71246 -1.032102 0.5698417 0.3048295
3 17.12370 -1.032102 0.5931174 0.3048295
4 19.03379 -1.032102 0.5077705 0.3048295
5 17.94374 -1.032102 0.4245812 0.3048295
6 16.54648 -1.032102 0.4314790 0.3048295
7 17.43091 -1.032102 0.5030304 0.3048295
8 18.78605 -1.032102 0.4077701 0.3048295
9 16.83386 -1.032102 0.4047545 0.3048295
10 14.38866 -1.032102 0.4168255 0.3048295
11 19.56890 -1.032102 0.5810463 0.3048295
12 18.95465 -1.032102 0.5605701 0.3048295
13 16.22043 -1.032102 0.4559741 0.3048295
14 16.93970 -1.032102 0.4608586 0.3048295
15 20.34944 -1.032102 0.3385849 0.3048295
16 15.72742 -1.032102 0.4710608 0.3048295
17 18.72756 -1.032102 0.4689031 0.3048295
18 16.42605 -1.032102 0.4676119 0.3048295
19 17.73991 -1.032102 0.3556849 0.3048295
20 16.94052 -1.032102 0.5325544 0.3048295
21 20.82426 -1.032102 0.5302523 0.3048295
22 17.69862 -1.032102 0.3390182 0.3048295
23 16.22125 -1.032102 0.5276699 0.3048295
24 13.49113 -1.032102 0.7815529 0.3048295
25 18.04875 -1.032102 0.4089895 0.3048295
26 17.06439 -1.032102 0.5825546 0.3048295
27 16.58956 -1.032102 0.3908872 0.3048295
## Define a function to get cov(Yi)
## V matrix: Cov(Yi) = Zi*G*Zi’ + sigma^2*Ini
GetCovYi <- function(RES.LMER, i = 1, RANEF = c(1)) {
## G matrix
G <- VarCorr(RES.LMER)$id
## Select columns of the design matrix
Z <- model.matrix(RES.LMER)[, RANEF, drop = F]
## Design matrix for random effect for subject i
Zi <- Z[RES.LMER@frame$id == i,]
## Within-individual variability
sigma <- sigma(RES.LMER)
Ini <- diag(nrow(Zi))
## Get Cov(Yi)
Vi <- Zi %*% G %*% t(Zi) + sigma^2 * Ini
Vi
}
## Get cov(Yi) for i = 1
GetCovYi(res.lmer, i = 1, RANEF = c(1,3))
[,1] [,2] [,3] [,4]
[1,] 4.950167 3.175099 3.116235 3.057372
[2,] 3.175099 4.962538 3.317566 3.388800
[3,] 3.116235 3.317566 5.235103 3.720229
[4,] 3.057372 3.388800 3.720229 5.767863
## Load MCMCglmm
library(MCMCglmm)
## Perform analysis
res.MCMCglmm <- MCMCglmm(fixed = distance ~ gender + age + gender:age,
random = ~ us(age):id,
family = "gaussian",
mev = NULL,
data = dental.long,
start = NULL, prior = NULL, tune = NULL, pedigree = NULL, nodes = "ALL",
scale = TRUE, nitt = 13000, thin = 10, burnin = 3000, pr = FALSE,
pl = FALSE, verbose = TRUE, DIC = TRUE, singular.ok = FALSE, saveX = TRUE,
saveZ = TRUE, saveXL = TRUE, slice = FALSE, ginverse = NULL)
MCMC iteration = 0
MCMC iteration = 1000
MCMC iteration = 2000
MCMC iteration = 3000
MCMC iteration = 4000
MCMC iteration = 5000
MCMC iteration = 6000
MCMC iteration = 7000
MCMC iteration = 8000
MCMC iteration = 9000
MCMC iteration = 10000
MCMC iteration = 11000
MCMC iteration = 12000
MCMC iteration = 13000
## Show results
summary(res.MCMCglmm)
Iterations = 3001:12991
Thinning interval = 10
Sample size = 1000
DIC: 409.8364
G-structure: ~us(age):id
post.mean l-95% CI u-95% CI eff.samp
age:age.id 0.02746 0.01235 0.04665 1000
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 2.073 1.487 2.802 1106
Location effects: distance ~ gender + age + gender:age
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 17.35360 15.41198 19.44500 1000 <0.001 ***
genderM -1.04659 -3.75004 1.49614 1000 0.404
age 0.48072 0.27944 0.69142 1000 0.002 **
genderM:age 0.30618 0.02945 0.56566 1000 0.020 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Plot diagnostics
plot(res.MCMCglmm)
## Fit model with unstructured covariance structure
res.gls <- gls(model = distance ~ gender + age + gender:age,
data = dental.long,
## Covariance structure: symmetry is the only restriction
correlation = corSymm(form = ~ as.numeric(factor(age)) | id),
## Variance structure: heterogeneous
weights = varIdent(form = ~ 1 | factor(age)),
method = "REML"
)
summary(res.gls)
Generalized least squares fit by REML
Model: distance ~ gender + age + gender:age
Data: dental.long
AIC BIC logLik
452.5468 489.5683 -212.2734
Correlation Structure: General
Formula: ~as.numeric(factor(age)) | id
Parameter estimate(s):
Correlation:
1 2 3
2 0.568
3 0.659 0.581
4 0.522 0.725 0.740
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | factor(age)
Parameter estimates:
8 10 12 14
1.0000000 0.8788793 1.0744596 0.9586878
Coefficients:
Value Std.Error t-value p-value
(Intercept) 17.425369 1.1726425 14.859916 0.0000
genderM -1.583086 1.5233073 -1.039243 0.3011
age 0.476365 0.0991583 4.804083 0.0000
genderM:age 0.350439 0.1288104 2.720580 0.0076
Correlation:
(Intr) gendrM age
genderM -0.770
age -0.875 0.674
genderM:age 0.674 -0.875 -0.770
Standardized residuals:
Min Q1 Med Q3 Max
-2.34272828 -0.63481507 -0.07904148 0.63772265 2.16523297
Residual standard error: 2.329213
Degrees of freedom: 108 total; 104 residual
intervals(res.gls)
Approximate 95% confidence intervals
Coefficients:
lower est. upper
(Intercept) 15.09997483 17.4253690 19.7507631
genderM -4.60386190 -1.5830863 1.4376893
age 0.27973007 0.4763647 0.6729993
genderM:age 0.09500314 0.3504390 0.6058749
attr(,"label")
[1] "Coefficients:"
Correlation structure:
lower est. upper
cor(1,2) 0.2530506 0.5681970 0.7743263
cor(1,3) 0.3843126 0.6589493 0.8264396
cor(1,4) 0.1851099 0.5220393 0.7491473
cor(2,3) 0.2731591 0.5806064 0.7804339
cor(2,4) 0.4823896 0.7249210 0.8642209
cor(3,4) 0.5135371 0.7396207 0.8696791
attr(,"label")
[1] "Correlation structure:"
Variance function:
lower est. upper
10 0.6328974 0.8788793 1.220464
12 0.7998338 1.0744596 1.443379
14 0.6830031 0.9586878 1.345649
attr(,"label")
[1] "Variance function:"
Residual standard error:
lower est. upper
1.760799 2.329213 3.081118
## Variance-Covariance structure compatible with SAS
getVarCov(res.gls)
Marginal variance covariance matrix
[,1] [,2] [,3] [,4]
[1,] 5.4252 2.7092 3.8411 2.7152
[2,] 2.7092 4.1906 2.9745 3.3137
[3,] 3.8411 2.9745 6.2632 4.1333
[4,] 2.7152 3.3137 4.1333 4.9862
Standard Deviations: 2.3292 2.0471 2.5026 2.233
See larger standard errors due to overestimation of variance in difference.
\( Var(a - b) = Var(a) + Var(b) - 2Cov(a,b) \)
The subtraction of the positive covariance (positive correlation in longitudinal data) at the end is ignored is data are assumed to be independent, giving a falsely high variance.
res.lm <- lm(formula = distance ~ 1 + gender + age + gender:age, data = dental.long)
summary(res.lm)
Call:
lm(formula = distance ~ 1 + gender + age + gender:age, data = dental.long)
Residuals:
Min 1Q Median 3Q Max
-5.6156 -1.3219 -0.1682 1.3299 5.2469
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.3727 1.7080 10.171 < 2e-16 ***
genderM -1.0321 2.2188 -0.465 0.64279
age 0.4795 0.1522 3.152 0.00212 **
genderM:age 0.3048 0.1977 1.542 0.12608
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.257 on 104 degrees of freedom
Multiple R-squared: 0.4227, Adjusted R-squared: 0.4061
F-statistic: 25.39 on 3 and 104 DF, p-value: 2.108e-12