1 Description

MNREAD Acuity Charts (Optelec, US Inc.) were used to measure reading speed as a function of print size. MNREAD data were collected from 42 adult observers with normal vision. Reading speed was measured in log word per minute. Print size was measured in in logMAR (related to degrees in visual angle).

Source: Cheung, C., Kallie, S., Legge, G.E., & Cheong, A.M.Y. (2008). Nonlinear Mixed-Effects Modeling of MNREAD Data. Investigative Ophthalmology & Visual Science, 49(2), 828-835.

Column 1: Subject ID
Column 2: Reading time in logWPM
Column 3: Print size measured in logMAR

# set digit and stars options
options(digits = 4, show.signif.stars = FALSE)

2 Data management

# set path to working directory
setwd("C:/Users/HANK/Desktop/HOMEWORK/")
# regular contrast polarity (i.e., black on white)
dta <- read.table("C:/Users/HANK/Desktop/HOMEWORK//mnread0.txt", header = T)
# load packages for later use
pacman::p_load(tidyverse, nlme)
# set plot theme to black-n-white
ot <- theme_set(theme_bw())
# individual profiles on separate panels
ggplot(dta, aes(pSize, logRS)) +
 geom_point(size = rel(.8), pch = 20)+
 geom_line()+
 facet_wrap( ~ ID) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

#  individual profiles on one panel
ggplot(dta, aes(pSize, logRS, group = ID)) +
 geom_point(pch = 20)+
 geom_line(alpha = .3) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

# coerce data into groupedData object
dta_g <- groupedData(logRS ~ pSize | ID, data = as.data.frame(dta),
                     labels = list(x = "Print size", y = "Reading speed"),
                     units = list(x = "(logMAR)", y = "(logWPM)") )
# lattice plot by subject
plot(dta_g, pch = 20, aspect = .89)

# try nls for all subjects
m0_g <- nlsList(SSasympOff, data = dta_g)
## Warning: 2 times caught the same error in nls(y ~ cbind(1, exp(-exp(lrc) * x)),
## data = xy, algorithm = "plinear", start = list(lrc = lrc)): singular gradient
# results of parameter estimates
coef(m0_g)
##       Asym   lrc       c0
## S128 2.131 1.590 -0.30739
## S132 2.186 2.388 -0.40812
## S102 2.154 2.350 -0.04337
## S131 2.193 1.551 -0.60031
## S122 2.166 2.270 -0.34818
## S109 2.202 3.196 -0.01809
## S108 2.171 2.353 -0.23285
## S130 2.187 2.332 -0.31713
## S106 2.184 2.372 -0.38127
## S117 2.212 1.819 -0.22989
## S136 2.275 2.769 -0.08696
## S103    NA    NA       NA
## S113 2.240 2.177 -0.21904
## S101 2.251 2.758 -0.03504
## S123 2.274 1.853 -0.30819
## S133 2.225 2.362 -0.44646
## S112 2.270 2.075 -0.23090
## S142 2.258 2.169 -0.33552
## S118 2.267 2.607 -0.32027
## S104 2.236 2.431 -0.24719
## S120 2.268 1.521 -0.31138
## S105 2.289 2.420 -0.25797
## S124 2.290 1.788 -0.16797
## S129 2.258 2.130 -0.32087
## S134 2.290 2.236 -0.37175
## S135 2.287 2.071 -0.29440
## S125 2.272 1.961 -0.28334
## S126 2.317 1.685 -0.47288
## S141 2.323 2.082 -0.29600
## S111 2.329 2.233 -0.42614
## S115 2.306 1.623 -0.56310
## S110 2.288 2.401 -0.14867
## S107 2.321 1.720 -0.37909
## S114 2.334 2.406 -0.24746
## S139    NA    NA       NA
## S137 2.334 2.622 -0.26045
## S140 2.352 2.190 -0.22763
## S138 2.325 2.358 -0.37132
## S127 2.316 1.982 -0.32574
## S119 2.386 1.902 -0.30892
## S116 2.412 1.710 -0.41221
## S121 2.456 2.476 -0.22895
# random effects model equivalent
m1_g <- nlme(m0_g)
# show summary
summary(m1_g)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: logRS ~ SSasympOff(pSize, Asym, lrc, c0) 
##  Data: dta_g 
##     AIC   BIC logLik
##   -1369 -1324  694.7
## 
## Random effects:
##  Formula: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
##  Level: ID
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev  Corr         
## Asym     0.06717 Asym   lrc   
## lrc      0.33175 -0.032       
## c0       0.10966 -0.103  0.492
## Residual 0.06186              
## 
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1) 
##        Value Std.Error  DF t-value p-value
## Asym  2.2700   0.01079 626  210.47       0
## lrc   2.2028   0.05630 626   39.13       0
## c0   -0.2856   0.01739 626  -16.43       0
##  Correlation: 
##     Asym   lrc   
## lrc -0.072       
## c0  -0.111  0.514
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -6.62871 -0.39007  0.08947  0.56091  2.61865 
## 
## Number of Observations: 670
## Number of Groups: 42
# fix c0
m1a_g <- update(m1_g, random = Asym + lrc ~ 1)
# model without random effects for the covariance terms involving Asym
m2_g <- update(m1_g, random = pdBlocked(list(Asym ~ 1, lrc + c0 ~ 1)))

3 Model comparison

# model comparison
anova(m1_g, m2_g, m1a_g)
##       Model df     AIC     BIC logLik   Test L.Ratio p-value
## m1_g      1 10 -1369.5 -1324.4  694.7                       
## m2_g      2  8 -1373.0 -1337.0  694.5 1 vs 2     0.4  0.8178
## m1a_g     3  7  -551.2  -519.7  282.6 2 vs 3   823.8  <.0001

4 Plot Fitted values and Residuals

# residuals
plot(m1a_g, id = .05, cex = .6, adj=-.05)

# random effects
plot(ranef(m1a_g, augFrame = T))

# parameter estimates CIs
intervals(m1a_g)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##        lower    est.   upper
## Asym  2.2820  2.3055  2.3289
## lrc   1.4610  1.5918  1.7226
## c0   -0.4176 -0.4096 -0.4016
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: ID 
##                  lower     est.    upper
## sd(Asym)       0.04262  0.06013  0.08484
## sd(lrc)        0.31114  0.39712  0.50684
## cor(Asym,lrc) -0.69649 -0.42048 -0.03608
## 
##  Within-group standard error:
##  lower   est.  upper 
## 0.1336 0.1415 0.1498
# plot model fits
dta_m1a <- data.frame(dta, fit = fitted(m1a_g), fit1 = fitted(m1_g))
#
ggplot(dta_m1a, aes(pSize, fit)) +
 geom_point(aes(pSize, logRS), pch = 1)+
 geom_line(col ="steelblue", lwd = 1)+
 geom_line(aes(pSize, fit1), col = "indianred", lwd = 1) +
 facet_wrap( ~ ID) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

#  individual profiles on one panel
ggplot(dta_m1a, aes(pSize, fit, group = ID)) +
 geom_point(aes(pSize, logRS), pch = 1, alpha = .6)+
 geom_line(alpha = .3, col="steelblue")+
 geom_line(aes(pSize, fit1), col = "indianred", alpha = .3) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

4.0.1