Skip to content

Commit 9de1fb5

Browse files
julianfarawayjulianfaraway
committed
vision freq
1 parent 6683bcd commit 9de1fb5

File tree

3 files changed

+561
-1
lines changed

3 files changed

+561
-1
lines changed

‎README.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -50,4 +50,5 @@ we explore below using the following packages:
5050
- [Nested Effects](mixed/eggsfreq.md) - the `eggs` data
5151
- [Crossed Effects](mixed/abrafreq.md) - the `abrasion` data
5252
- [Multilevel Models](mixed/jspmlevelfreq.md) - the `jsp` data
53-
- [Longitudinal Models](mixed/longitfreq.md) - the `psid` data
53+
- [Longitudinal Models](mixed/longitfreq.md) - the `psid` data
54+
- [Repeated Measures](mixed/visionfreq.md) - the `vision` data

‎mixed/visionfreq.md‎

Lines changed: 355 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,355 @@
1+
# Repeated Measures with Vision Data using Frequentist Methods
2+
[Julian Faraway](https://julianfaraway.github.io/)
3+
2024-10-09
4+
5+
- [Data](#data)
6+
- [Mixed Effect Model](#mixed-effect-model)
7+
- [LME4](#lme4)
8+
- [NLME](#nlme)
9+
- [MMRM](#mmrm)
10+
- [GLMMTMB](#glmmtmb)
11+
- [Discussion](#discussion)
12+
- [Package version info](#package-version-info)
13+
14+
See the [introduction](index.md) for an overview.
15+
16+
See a [mostly Bayesian analysis](vision.md) analysis of the same data.
17+
18+
This example is discussed in more detail in my book [Extending the
19+
Linear Model with R](https://julianfaraway.github.io/faraway/ELM/)
20+
21+
Required libraries:
22+
23+
``` r
24+
library(faraway)
25+
library(ggplot2)
26+
```
27+
28+
# Data
29+
30+
The acuity of vision for seven subjects was tested. The response is the
31+
lag in milliseconds between a light flash and a response in the cortex
32+
of the eye. Each eye is tested at four different powers of lens. An
33+
object at the distance of the second number appears to be at distance of
34+
the first number.
35+
36+
Load in and look at the data:
37+
38+
``` r
39+
data(vision, package="faraway")
40+
ftable(xtabs(acuity ~ eye + subject + power, data=vision))
41+
```
42+
43+
power 6/6 6/18 6/36 6/60
44+
eye subject
45+
left 1 116 119 116 124
46+
2 110 110 114 115
47+
3 117 118 120 120
48+
4 112 116 115 113
49+
5 113 114 114 118
50+
6 119 115 94 116
51+
7 110 110 105 118
52+
right 1 120 117 114 122
53+
2 106 112 110 110
54+
3 120 120 120 124
55+
4 115 116 116 119
56+
5 114 117 116 112
57+
6 100 99 94 97
58+
7 105 105 115 115
59+
60+
We create a numeric version of the power to make a plot:
61+
62+
``` r
63+
vision$npower <- rep(1:4,14)
64+
ggplot(vision, aes(y=acuity, x=npower, linetype=eye)) + geom_line() + facet_wrap(~ subject, ncol=4) + scale_x_continuous("Power",breaks=1:4,labels=c("6/6","6/18","6/36","6/60"))
65+
```
66+
67+
![](figs/visplot-1..svg)
68+
69+
# Mixed Effect Model
70+
71+
The power is a fixed effect. In the model below, we have treated it as a
72+
nominal factor, but we could try fitting it in a quantitative manner.
73+
The subjects should be treated as random effects. Since we do not
74+
believe there is any consistent right-left eye difference between
75+
individuals, we should treat the eye factor as nested within subjects.
76+
77+
The model can be written as:
78+
79+
``` math
80+
y_{ijk} = \mu + p_j + s_i + e_{ik} + \epsilon_{ijk}
81+
```
82+
83+
where $i=1, \dots ,7$ runs over individuals, $j=1, \dots ,4$ runs over
84+
power and $k=1,2$ runs over eyes. The $p_j$ term is a fixed effect, but
85+
the remaining terms are random. Let $s_i \sim N(0,\sigma^2_s)$,
86+
$e_{ik} \sim N(0,\sigma^2_e)$ and
87+
$\epsilon_{ijk} \sim N(0,\sigma^2\Sigma)$ where we take $\Sigma=I$.
88+
89+
# LME4
90+
91+
``` r
92+
library(lme4)
93+
library(pbkrtest)
94+
```
95+
96+
We can fit the model with:
97+
98+
``` r
99+
mmod <- lmer(acuity~power + (1|subject) + (1|subject:eye),vision)
100+
faraway::sumary(mmod)
101+
```
102+
103+
Fixed Effects:
104+
coef.est coef.se
105+
(Intercept) 112.64 2.23
106+
power6/18 0.79 1.54
107+
power6/36 -1.00 1.54
108+
power6/60 3.29 1.54
109+
110+
Random Effects:
111+
Groups Name Std.Dev.
112+
subject:eye (Intercept) 3.21
113+
subject (Intercept) 4.64
114+
Residual 4.07
115+
---
116+
number of obs: 56, groups: subject:eye, 14; subject, 7
117+
AIC = 342.7, DIC = 349.6
118+
deviance = 339.2
119+
120+
We can check for a power effect using a Kenward-Roger adjusted $F$-test:
121+
122+
``` r
123+
mmod <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision,REML=FALSE)
124+
nmod <- lmer(acuity~1+(1|subject)+(1|subject:eye),vision,REML=FALSE)
125+
KRmodcomp(mmod, nmod)
126+
```
127+
128+
large : acuity ~ power + (1 | subject) + (1 | subject:eye)
129+
small : acuity ~ 1 + (1 | subject) + (1 | subject:eye)
130+
stat ndf ddf F.scaling p.value
131+
Ftest 2.83 3.00 39.00 1 0.051
132+
133+
The power just fails to meet the 5% level of significance (although note
134+
that there is a clear outlier in the data which deserves some
135+
consideration). We can also compute bootstrap confidence intervals:
136+
137+
``` r
138+
set.seed(123)
139+
print(confint(mmod, method="boot", oldNames=FALSE, nsim=1000),digits=3)
140+
```
141+
142+
2.5 % 97.5 %
143+
sd_(Intercept)|subject:eye 0.172 5.27
144+
sd_(Intercept)|subject 0.000 6.88
145+
sigma 2.943 4.70
146+
(Intercept) 108.623 116.56
147+
power6/18 -1.950 3.94
148+
power6/36 -3.868 1.88
149+
power6/60 0.388 6.22
150+
151+
We see that lower ends of the CIs for random effect SDs are zero or
152+
close to it.
153+
154+
# NLME
155+
156+
See the discussion for the [single random effect
157+
example](pulpfreq.md#NLME) for some introduction.
158+
159+
The syntax for specifying the nested/heirarchical model is different
160+
from `lme4`:
161+
162+
``` r
163+
library(nlme)
164+
nlmod = lme(acuity ~ power,
165+
vision,
166+
~ 1 | subject/eye)
167+
summary(nlmod)
168+
```
169+
170+
Linear mixed-effects model fit by REML
171+
Data: vision
172+
AIC BIC logLik
173+
342.71 356.37 -164.35
174+
175+
Random effects:
176+
Formula: ~1 | subject
177+
(Intercept)
178+
StdDev: 4.6396
179+
180+
Formula: ~1 | eye %in% subject
181+
(Intercept) Residual
182+
StdDev: 3.2052 4.0746
183+
184+
Fixed effects: acuity ~ power
185+
Value Std.Error DF t-value p-value
186+
(Intercept) 112.643 2.2349 39 50.401 0.0000
187+
power6/18 0.786 1.5400 39 0.510 0.6128
188+
power6/36 -1.000 1.5400 39 -0.649 0.5199
189+
power6/60 3.286 1.5400 39 2.134 0.0392
190+
Correlation:
191+
(Intr) pw6/18 pw6/36
192+
power6/18 -0.345
193+
power6/36 -0.345 0.500
194+
power6/60 -0.345 0.500 0.500
195+
196+
Standardized Within-Group Residuals:
197+
Min Q1 Med Q3 Max
198+
-3.424009 -0.323213 0.010949 0.440812 2.466186
199+
200+
Number of Observations: 56
201+
Number of Groups:
202+
subject eye %in% subject
203+
7 14
204+
205+
The results are presented somewhat differently but match those presented
206+
by `lme4` earlier. We do get p-values for the fixed effects but these
207+
are not so useful.
208+
209+
We can get tests on the fixed effects with:
210+
211+
``` r
212+
anova(nlmod)
213+
```
214+
215+
numDF denDF F-value p-value
216+
(Intercept) 1 39 3132.91 <.0001
217+
power 3 39 2.83 0.0511
218+
219+
In this case, the results are the same as with the `pbkrtest` output
220+
because the degrees of freedom adjustment has no effect. This is not
221+
always the case.
222+
223+
# MMRM
224+
225+
This package is not designed to fit models with a hierarchical
226+
structure. In principle, it should be possible to specify a
227+
parameterised covariance structure corresponding to this design but the
228+
would reduce to the previous computations.
229+
230+
# GLMMTMB
231+
232+
See the discussion for the [single random effect
233+
example](pulpfreq.md#GLMMTMB) for some introduction.
234+
235+
``` r
236+
library(glmmTMB)
237+
```
238+
239+
The default fit uses ML (not REML)
240+
241+
``` r
242+
gtmod <- glmmTMB(acuity~power+(1|subject)+(1|subject:eye),data=vision)
243+
summary(gtmod)
244+
```
245+
246+
Family: gaussian ( identity )
247+
Formula: acuity ~ power + (1 | subject) + (1 | subject:eye)
248+
Data: vision
249+
250+
AIC BIC logLik deviance df.resid
251+
353.2 367.4 -169.6 339.2 49
252+
253+
Random effects:
254+
255+
Conditional model:
256+
Groups Name Variance Std.Dev.
257+
subject (Intercept) 17.4 4.17
258+
subject:eye (Intercept) 10.6 3.25
259+
Residual 15.4 3.93
260+
Number of obs: 56, groups: subject, 7; subject:eye, 14
261+
262+
Dispersion estimate for gaussian family (sigma^2): 15.4
263+
264+
Conditional model:
265+
Estimate Std. Error z value Pr(>|z|)
266+
(Intercept) 112.643 2.084 54.0 <2e-16
267+
power6/18 0.786 1.484 0.5 0.596
268+
power6/36 -1.000 1.484 -0.7 0.500
269+
power6/60 3.286 1.484 2.2 0.027
270+
271+
Another option is to use REML with:
272+
273+
``` r
274+
gtmodr = glmmTMB(acuity~power+(1|subject)+(1|subject:eye),data=vision,
275+
REML=TRUE)
276+
summary(gtmodr)
277+
```
278+
279+
Family: gaussian ( identity )
280+
Formula: acuity ~ power + (1 | subject) + (1 | subject:eye)
281+
Data: vision
282+
283+
AIC BIC logLik deviance df.resid
284+
342.7 356.9 -164.4 328.7 53
285+
286+
Random effects:
287+
288+
Conditional model:
289+
Groups Name Variance Std.Dev.
290+
subject (Intercept) 21.5 4.64
291+
subject:eye (Intercept) 10.3 3.21
292+
Residual 16.6 4.07
293+
Number of obs: 56, groups: subject, 7; subject:eye, 14
294+
295+
Dispersion estimate for gaussian family (sigma^2): 16.6
296+
297+
Conditional model:
298+
Estimate Std. Error z value Pr(>|z|)
299+
(Intercept) 112.643 2.235 50.4 <2e-16
300+
power6/18 0.786 1.540 0.5 0.610
301+
power6/36 -1.000 1.540 -0.6 0.516
302+
power6/60 3.286 1.540 2.1 0.033
303+
304+
The result is appears identical with the previous REML fits.
305+
306+
If we want to test the significance of the `power` fixed effect, we have
307+
the same methods available as for the `lme4` fit.
308+
309+
# Discussion
310+
311+
See the [Discussion of the single random effect
312+
model](pulp.md#Discussion) for general comments.
313+
314+
`lme4`, `nlme` and `glmmTMB` were all able to fit this model. `mmrm` was
315+
not designed for this type of model.
316+
317+
# Package version info
318+
319+
``` r
320+
sessionInfo()
321+
```
322+
323+
R version 4.4.1 (2024-06-14)
324+
Platform: x86_64-apple-darwin20
325+
Running under: macOS Sonoma 14.7
326+
327+
Matrix products: default
328+
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
329+
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
330+
331+
locale:
332+
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
333+
334+
time zone: Europe/London
335+
tzcode source: internal
336+
337+
attached base packages:
338+
[1] stats graphics grDevices utils datasets methods base
339+
340+
other attached packages:
341+
[1] glmmTMB_1.1.9 nlme_3.1-166 pbkrtest_0.5.3 lme4_1.1-35.5 Matrix_1.7-0 ggplot2_3.5.1 faraway_1.0.8
342+
343+
loaded via a namespace (and not attached):
344+
[1] utf8_1.2.4 generics_0.1.3 tidyr_1.3.1 lattice_0.22-6 digest_0.6.37
345+
[6] magrittr_2.0.3 estimability_1.5.1 evaluate_0.24.0 grid_4.4.1 mvtnorm_1.2-6
346+
[11] fastmap_1.2.0 jsonlite_1.8.8 backports_1.5.0 mgcv_1.9-1 purrr_1.0.2
347+
[16] fansi_1.0.6 scales_1.3.0 numDeriv_2016.8-1.1 cli_3.6.3 rlang_1.1.4
348+
[21] munsell_0.5.1 splines_4.4.1 withr_3.0.1 yaml_2.3.10 tools_4.4.1
349+
[26] parallel_4.4.1 coda_0.19-4.1 nloptr_2.1.1 minqa_1.2.8 dplyr_1.1.4
350+
[31] colorspace_2.1-1 boot_1.3-31 broom_1.0.6 vctrs_0.6.5 R6_2.5.1
351+
[36] emmeans_1.10.4 lifecycle_1.0.4 MASS_7.3-61 pkgconfig_2.0.3 pillar_1.9.0
352+
[41] gtable_0.3.5 glue_1.7.0 Rcpp_1.0.13 systemfonts_1.1.0 xfun_0.47
353+
[46] tibble_3.2.1 tidyselect_1.2.1 rstudioapi_0.16.0 knitr_1.48 xtable_1.8-4
354+
[51] farver_2.1.2 htmltools_0.5.8.1 rmarkdown_2.28 svglite_2.1.3 labeling_0.4.3
355+
[56] TMB_1.9.14 compiler_4.4.1

0 commit comments

Comments
 (0)