1 setting

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

2 data input

# regular contrast polarity (i.e., black on white)
dta <- read.table("C:/Users/user/Desktop/multilevel/1130/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())

3 individual profiles plot

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

42個人，每個人print size對 reading speed的變化curve

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

把42個人的curve畫在一起，可以看得出有個體間的差異(subject random effect)，尤其是print size越小的時候

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

groupedData object看起來與ggplot畫起來的curve差不多，不過它是按照斜率排列，而非subject number

# try nls for all subjects
# use nlsList function，data must inherit from class "groupedData".
# SSasympOff means Self-Starting Nls Asymptotic Regression Model with an Offset
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

Asym = a numeric parameter representing the horizontal asymptote on the right side

lrc = a numeric parameter representing the natural logarithm of the rate constant.

c0 = a numeric parameter representing the input for which the response is zero.

# 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)
summary(m1a_g)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: logRS ~ SSasympOff(pSize, Asym, lrc, c0) 
##  Data: dta_g 
##      AIC    BIC logLik
##   -551.2 -519.7  282.6
## 
## Random effects:
##  Formula: list(Asym ~ 1, lrc ~ 1)
##  Level: ID
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev  Corr 
## Asym     0.06013 Asym 
## lrc      0.39712 -0.42
## Residual 0.14145      
## 
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1) 
##        Value Std.Error  DF t-value p-value
## Asym  2.3055   0.01197 626  192.57       0
## lrc   1.5918   0.06675 626   23.85       0
## c0   -0.4096   0.00410 626 -100.00       0
##  Correlation: 
##     Asym   lrc   
## lrc -0.430       
## c0  -0.103  0.266
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -5.67772 -0.33577  0.02076  0.39501  4.07914 
## 
## Number of Observations: 670
## Number of Groups: 42
# model without random effects for the covariance terms involving Asym
m2_g <- update(m1_g, random = pdBlocked(list(Asym ~ 1, lrc + c0 ~ 1)))
summary(m2_g)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: logRS ~ SSasympOff(pSize, Asym, lrc, c0) 
##  Data: dta_g 
##     AIC   BIC logLik
##   -1373 -1337  694.5
## 
## Random effects:
##  Composite Structure: Blocked
## 
##  Block 1: Asym
##  Formula: Asym ~ 1 | ID
##            Asym
## StdDev: 0.06714
## 
##  Block 2: lrc, c0
##  Formula: list(lrc ~ 1, c0 ~ 1)
##  Level: ID
##  Structure: General positive-definite
##          StdDev  Corr
## lrc      0.33153 lrc 
## c0       0.10968 0.49
## Residual 0.06185     
## 
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1) 
##        Value Std.Error  DF t-value p-value
## Asym  2.2701   0.01078 626  210.58       0
## lrc   2.2026   0.05626 626   39.15       0
## c0   -0.2857   0.01739 626  -16.43       0
##  Correlation: 
##     Asym   lrc   
## lrc -0.044       
## c0  -0.014  0.512
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -6.63101 -0.38734  0.09036  0.55052  2.61625 
## 
## Number of Observations: 670
## Number of Groups: 42
# 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

m1a_g顯著比m2_g來得好

m1_g和m2_g差不多

# 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

標準化殘差圖結果顯示，看起來fit不是很好。

4 plot model fits

dta_m1a <- data.frame(dta, fit = fitted(m1a_g), fit1 = fitted(m1_g))

#individual profiles on separate panel
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)")

5 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)")

###

比較兩個fit 的curve