Exercise

1 Exercise 1

1.1 load and see

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

# regular contrast polarity (i.e., black on white)
dta1 <- read.table("mnread0.txt", header = T)
str(dta1)

## 'data.frame':    670 obs. of  3 variables:
##  $ ID   : chr  "S101" "S101" "S101" "S101" ...
##  $ pSize: num  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
##  $ logRS: num  0.951 2.028 2.092 2.219 2.292 ...

head(dta1)

##     ID pSize  logRS
## 1 S101   0.0 0.9512
## 2 S101   0.1 2.0276
## 3 S101   0.2 2.0924
## 4 S101   0.3 2.2194
## 5 S101   0.4 2.2924
## 6 S101   0.5 2.2206

# load packages for later use
pacman::p_load(tidyverse, nlme)

1.2 individual profiles on separate panels

# set plot theme to black-n-white
ot <- theme_set(theme_bw())

# individual profiles on separate panels
ggplot(dta1, 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)")

1.3 individual profiles on one panel

#  individual profiles on one panel
ggplot(dta1, 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(dta1),
                     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)))

1.4 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 is better

1.5 residuals

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

1.6 random effects

plot(ranef(m1a_g, augFrame = T))

1.7 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

1.8 plot model fits

dta_m1a <- data.frame(dta1, 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)")

###

2 Exercise 2

2.1 load and see

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

# load packages for later use
pacman::p_load(bdots, tidyverse, nlme)

# load data
data(ci)
str(ci)

## 'data.frame':    108216 obs. of  5 variables:
##  $ protocol : Factor w/ 2 levels "CI","NH": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Subject  : int  36 36 36 36 36 36 36 36 36 36 ...
##  $ Time     : int  0 4 8 12 16 20 24 28 32 36 ...
##  $ Fixations: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ LookType : Factor w/ 4 levels "Cohort","Rhyme",..: 1 1 1 1 1 1 1 1 1 1 ...

head(ci)

##   protocol Subject Time Fixations LookType
## 1       NH      36    0         0   Cohort
## 2       NH      36    4         0   Cohort
## 3       NH      36    8         0   Cohort
## 4       NH      36   12         0   Cohort
## 5       NH      36   16         0   Cohort
## 6       NH      36   20         0   Cohort

2.2 plot

# set plot themes to black-n-white
ot <- theme_set(theme_bw())

#ci$Subject <- as.factor(ci$Subject)
ggplot(ci, aes(Time, Fixations)) +
 geom_point(size=rel(.1), pch = 1, alpha = .2) +
 facet_grid(protocol ~ LookType)

2.3 change form

# select target looking time course data
dta_t <- subset(ci, LookType == "Target")

# order data by subject , time
dta_t <- dta_t[order(dta_t$Subject, dta_t$Time, decreasing = F), ]

# have a look
head(dta_t)

##       protocol Subject Time Fixations LookType
## 39079       CI       2    0         0   Target
## 39080       CI       2    4         0   Target
## 39081       CI       2    8         0   Target
## 39082       CI       2   12         0   Target
## 39083       CI       2   16         0   Target
## 39084       CI       2   20         0   Target

2.4 convert to grouped Data object for nlme processing

dtag <- groupedData(Fixations ~ Time | Subject, outer = ~ protocol,
                    data = as.data.frame(dta_t),
                    labels = list(x = "Time", y = "Fixations"),
                    units = list(x = "(ms)", y = "(proportion)") )

# lattice plot by subject
plot(dtag, pch = 20, aspect = .89)

# plot by subject given group = protocol
p <- plot(dtag, outer = TRUE, cex = .1, layout = c(2,1), aspect = .89)
update(p, legend = NULL)

2.5 fitting function

Model: f(x) = P / (1 + exp(S*(C - x)/P))

# define the 3-parameter curve fitting function
fx <- function(x, P, S, C) {
               P / (1 + exp(S*(C - x)/P))
}
#

2.6 m0

m0 <- nlsList(Fixations ~ fx(Time, P, S, C),
               start = c(P = 0.87, S = 0.004, C = 757.8),
               data = dtag)

# results of parameter estimates
coef(m0)

##          P        S     C
## 24  0.5571 0.004108 708.9
## 81  0.7501 0.004582 804.5
## 6   0.7480 0.004294 801.2
## 90  0.7795 0.004376 780.4
## 2   0.7770 0.003588 811.9
## 41  0.8007 0.006199 728.2
## 84  0.8190 0.004794 792.4
## 30  0.8259 0.005005 819.4
## 16  0.8409 0.004715 804.8
## 10  0.8438 0.004999 788.8
## 23  0.8913 0.005303 924.2
## 75  0.8929 0.005409 805.3
## 79  0.9074 0.006026 698.8
## 48  0.8881 0.004349 792.9
## 76  0.8908 0.005846 761.5
## 4   0.9082 0.006417 736.7
## 74  0.9020 0.004945 762.7
## 72  0.8887 0.007259 822.3
## 17  0.9086 0.005659 788.8
## 29  0.9401 0.006854 777.9
## 19  0.9224 0.007384 709.1
## 82  0.9483 0.007750 734.1
## 83  0.9332 0.009380 761.1
## 94  0.9353 0.008019 735.2
## 22  0.9474 0.005805 814.3
## 86  0.9674 0.007940 690.5
## 99  0.9681 0.007187 817.4
## 3   0.9716 0.007568 677.7
## 68  0.6759 0.004333 745.9
## 64  0.8652 0.007813 761.6
## 55  0.8759 0.007716 711.8
## 98  0.8861 0.005607 744.8
## 100 0.8969 0.007230 795.8
## 51  0.9269 0.006000 705.6
## 105 0.9254 0.007172 685.0
## 38  0.9375 0.008724 696.7
## 40  0.9477 0.007433 721.8
## 65  0.9450 0.009253 671.1
## 42  0.9622 0.008466 622.0
## 71  0.9539 0.008714 658.7
## 63  0.9633 0.007947 671.7
## 95  0.9625 0.007309 725.4
## 66  0.9669 0.008559 777.7
## 106 0.9630 0.005560 743.0
## 56  0.9674 0.007671 681.7
## 60  0.9673 0.008303 680.0
## 85  0.9713 0.007033 690.3
## 102 0.9636 0.008226 722.6
## 36  0.9795 0.007957 698.8
## 61  0.9684 0.007889 683.7
## 73  0.9773 0.008403 673.0
## 47  0.9812 0.010171 704.2
## 53  0.9811 0.009324 722.4
## 58  0.9768 0.007508 682.5

2.7 m1

# covert fitting separately by subject to subjects as random effects
m1 <- nlme(m0)

## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: false
## convergence (8)

## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 2, LME step: nlminb() did not converge (code = 1). PORT message: false
## convergence (8)

2.8 mean parameters values

# extract mean parameters values to use as initial values
prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)

2.9 m2

# update to let parameters dependent on protocol
# and simplify the covariance structure to a simple separate
# variances for each one of the parameters
#

m2 <- update(m1, 
             fixed = list(P + S + C ~ protocol),
             random = pdDiag(P + S + C ~ 1),
             groups = ~ Subject, 
             start = list(fixed = c(P = c(prm0[1], 0), 
                                    S = c(prm0[2], 0), 
                                    C = c(prm0[3], 0))))

2.10 m3

m3 <- update(m2, corr = corAR1(0.95, form = ~ 1 | Subject),
             control = lmeControl(opt = 'optim', 
                                  niterEM = 500,
                                  maxIter = 1e6,
                                  msMaxIter = 1e3))

## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value

## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value

## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value

## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value

## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value

## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 1, LME step: nlm() did not converge (code = 4). Do increase
## 'msMaxIter'!

## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 2, LME step: nlm() did not converge (code = 4). Do increase
## 'msMaxIter'!

## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 4, LME step: nlm() did not converge (code = 3).

## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 5, LME step: nlm() did not converge (code = 3).

# show summary
summary(m3)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: Fixations ~ fx(Time, P, S, C) 
##  Data: dtag 
##       AIC     BIC logLik
##   -244799 -244709 122411
## 
## Random effects:
##  Formula: list(P ~ 1, S ~ 1, C ~ 1)
##  Level: Subject
##  Structure: Diagonal
##         P.(Intercept) S.(Intercept) C.(Intercept) Residual
## StdDev:       0.07106      0.001399         46.43  0.02976
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | Subject 
##  Parameter estimate(s):
##    Phi 
## 0.9962 
## Fixed effects: list(P + S + C ~ protocol) 
##               Value Std.Error    DF t-value p-value
## P.(Intercept)   0.9     0.014 26995   61.37  0.0000
## P.protocolNH    0.1     0.021 26995    3.17  0.0015
## S.(Intercept)   0.0     0.000 26995   21.23  0.0000
## S.protocolNH    0.0     0.000 26995    4.49  0.0000
## C.(Intercept) 754.6     9.312 26995   81.04  0.0000
## C.protocolNH  -69.0    13.205 26995   -5.23  0.0000
##  Correlation: 
##               P.(In) P.prNH S.(In) S.prNH C.(In)
## P.protocolNH  -0.696                            
## S.(Intercept)  0.007 -0.005                     
## S.protocolNH  -0.005  0.008 -0.692              
## C.(Intercept)  0.012 -0.008  0.001  0.000       
## C.protocolNH  -0.008  0.008 -0.001  0.001 -0.705
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -3.30299 -1.01280 -0.47437 -0.06979  1.58923 
## 
## Number of Observations: 27054
## Number of Groups: 54

2.11 auto-correlation plot

#
plot(ACF(m3, resType = "n"), alpha = 0.05)

2.12 residual plot

# residual plot to identify outlier
plot(m3, id = 0.05, pch = '.', adj = -0.1, cex = .4)

# residual plot by protocol
plot(m3, resid(., type = "p") ~ fitted(.) | protocol, abline = 0, pch = '.')

# residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | Subject, abline = 0, pch = '.')

# residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | Subject*protocol, 
     abline = 0, pch = '.')

2.13 model fit

# plot fitted lines over observations by subject
plot(augPred(m3, level = 0:1))

# fortify data with fitted values
dta_t$yhat <- fitted(m3)

# Normal group
ggplot(subset(dta_t, protocol == "NH"), aes(Time, Fixations)) +
 geom_point(size = rel(.5), pch = 1, col = "skyblue") +
 geom_line(aes(y = yhat, group = Subject)) +
 coord_cartesian(ylim = c(0, 1)) +
 facet_wrap(~ Subject, ncol = 7) +
 labs(x = "Time (ms)", y = "Proportion of fixations") +
 theme(legend.position = "NONE")

# implant group
ggplot(subset(dta_t, protocol == "CI"), aes(Time, Fixations)) +
 geom_point(size = rel(.5), pch = 1, col = "pink") +
 geom_line(aes(y = yhat, group = Subject)) +
 coord_cartesian(ylim = c(0, 1)) +
 facet_wrap(~ Subject, ncol = 7) +
 labs(x = "Time (ms)", y = "Proportion of fixations") +
 theme(legend.position = "NONE")

3 Exercise 3

3.1 load and see

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

# load package for later use - must have pacman installed first
pacman::p_load(tidyverse, nlme, car)

# check the data set
?Wong

## starting httpd help server ... done

# save a copy
dta <- Wong

# check data structure
str(dta)

## 'data.frame':    331 obs. of  7 variables:
##  $ id      : int  3358 3535 3547 3592 3728 3790 3807 3808 4253 4356 ...
##  $ days    : num  30 16 40 13 19 13 37 31 40 31 ...
##  $ duration: int  4 17 1 10 6 3 5 7 3 7 ...
##  $ sex     : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
##  $ age     : num  20.7 55.3 55.9 61.7 30.1 ...
##  $ piq     : int  87 95 95 59 67 76 74 91 115 86 ...
##  $ viq     : int  89 77 116 73 73 69 77 110 110 83 ...

# see the first 6 lines
head(dta)

##     id days duration  sex   age piq viq
## 1 3358   30        4 Male 20.67  87  89
## 2 3535   16       17 Male 55.29  95  77
## 3 3547   40        1 Male 55.92  95 116
## 4 3592   13       10 Male 61.66  59  73
## 5 3728   19        6 Male 30.13  67  73
## 6 3790   13        3 Male 57.06  76  69

3.2 Plot

# plot individuals
ggplot(dta, aes(sqrt(days), piq, group = id)) +
 geom_point(pch = 20) +
 geom_line() +
 stat_smooth(aes(group=1), method = "gam", formula = y ~ s(x)) +
 labs(x = "Days after recovering from coma (square-root)",
      y = "Performance IQ score")

# examine the piq against coma duratin (sqrt)
ggplot(dta, aes(sqrt(duration), piq)) +
 geom_point(pch = 20) +
 stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) +
 labs(x = "Duration of coma (square-root)", y = "Performance IQ score")

隨著duration of coma增加，performannce IQ score下降

3.3 fit model

# use this to guess the initial values for parameters
summary(lm(piq ~ sqrt(duration), data = dta))

## 
## Call:
## lm(formula = piq ~ sqrt(duration), data = dta)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -40.30 -10.89  -1.73   9.33  45.34 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)      91.304      1.275   71.61  < 2e-16
## sqrt(duration)   -1.289      0.337   -3.82  0.00016
## 
## Residual standard error: 14.8 on 329 degrees of freedom
## Multiple R-squared:  0.0425, Adjusted R-squared:  0.0396 
## F-statistic: 14.6 on 1 and 329 DF,  p-value: 0.000158

# fit the model
m1c <- nlme(piq ~ b0 + b1*exp(-b2*days), data = dta,
           fixed = list(b0 ~ 1 + sqrt(duration),
                        b1 ~ 1 + sqrt(duration),
                        b2 ~ 1),
           random = list(id = list(b0 ~ 1, b1 ~ 1)),
           start = list(fixed = c(100, -1.5, -10, 0, 0.007)))

## Warning in nlme.formula(piq ~ b0 + b1 * exp(-b2 * days), data = dta, fixed =
## list(b0 ~ : Iteration 1, LME step: nlminb() did not converge (code = 1). Do
## increase 'msMaxIter'!

## Warning in nlme.formula(piq ~ b0 + b1 * exp(-b2 * days), data = dta, fixed =
## list(b0 ~ : Iteration 2, LME step: nlminb() did not converge (code = 1). Do
## increase 'msMaxIter'!

# show output summary
summary(m1c)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: piq ~ b0 + b1 * exp(-b2 * days) 
##  Data: dta 
##    AIC  BIC logLik
##   2593 2628  -1288
## 
## Random effects:
##  Formula: list(b0 ~ 1, b1 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##                StdDev Corr  
## b0.(Intercept) 13.769 b0.(I)
## b1.(Intercept)  2.600 -0.996
## Residual        6.736       
## 
## Fixed effects: list(b0 ~ 1 + sqrt(duration), b1 ~ 1 + sqrt(duration), b2 ~ 1) 
##                    Value Std.Error  DF t-value p-value
## b0.(Intercept)     97.09     2.037 127   47.68  0.0000
## b0.sqrt(duration)  -1.25     0.480 127   -2.59  0.0107
## b1.(Intercept)    -11.15     3.208 127   -3.47  0.0007
## b1.sqrt(duration)  -3.25     1.077 127   -3.02  0.0031
## b2                  0.01     0.002 127    5.00  0.0000
##  Correlation: 
##                   b0.(I) b0.s() b1.(I) b1.s()
## b0.sqrt(duration) -0.724                     
## b1.(Intercept)    -0.596  0.463              
## b1.sqrt(duration)  0.463 -0.455 -0.789       
## b2                -0.309  0.013  0.092 -0.380
## 
## Standardized Within-Group Residuals:
##       Min        Q1       Med        Q3       Max 
## -3.332390 -0.365703  0.008973  0.382757  2.303088 
## 
## Number of Observations: 331
## Number of Groups: 200

3.4 plot

# augment data with fitted values
dta_m1 <- data.frame(dta, fit = predict(m1c))

# fitted lines
ggplot(dta_m1, aes(sqrt(days), piq, group = id)) +
 geom_point(pch = 20, alpha = .5) +
 geom_line(aes(sqrt(days), fit, group = id)) +
 labs(x = "Days after recovering from coma (square-root)",
      y = "Performance IQ score")

# plot mode fits - basic graphics - detail contol
#
newdata <- expand.grid(duration = c(1, 10, 20, 50, 100, 200),
                       days = seq(0, 1000, 20))
newdata$piq <- predict(m1c, newdata, level = 0)

plot(piq ~ days, type="n", data = newdata,
     xlab = "Days Post Coma", ylab = "Average Performance IQ",
     ylim = c(15, 105), xlim = c(-100, 1000), axes = FALSE, frame=TRUE)
axis(2)
axis(4)
axis(1, at = seq(0, 1000, by=100))
for (dc in c(1, 10, 20, 50, 100, 200)){
 with(newdata, {
      lines(spline(seq(0, 1000, 20), piq[duration == dc]), lwd = 2)
      text(-25, piq[duration == dc][1], labels = dc, adj = 0.9)
 }
)}
grid()
text(-100, 95, labels="Length of\nComa", adj = 0, cex = .7)

Conclusion: 昏迷時間越長，表現智商得分越低。

4 Hw1

4.1 load and see

library(lattice)
library(nlme)
library(nlstools)

## 
## 'nlstools' has been loaded.

## IMPORTANT NOTICE: Most nonlinear regression models and data set examples

## related to predictive microbiolgy have been moved to the package 'nlsMicrobio'

dtah1 <-  read.table("freeRecall.asc.txt", header = T)

names(dtah1) <- c("Group","Trial","CorrectN")

str(dtah1)

## 'data.frame':    20 obs. of  3 variables:
##  $ Group   : chr  "C" "C" "C" "C" ...
##  $ Trial   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ CorrectN: num  7.9 10.9 11.9 13 14.2 14.2 14.7 15.1 14.8 15.2 ...

4.2 Different Groups’ correct responses in each Trial

ggplot(dtah1, aes(Trial, CorrectN, color = Group)) +
 geom_point(size = rel(2)) +
 geom_line() +
 labs(x = "Trial on list B", y = "Mean number of correct responses")+
 scale_x_continuous(limits = c(0, 10), breaks = seq(0, 10, by = 2)) +
 scale_y_continuous(limits = c(5, 20), breaks = seq(5, 20, by = 5)) +
 guides(color = guide_legend(reverse = T)) +
 theme(legend.position = c(.9,.2))  

4.3 model

Model: CRi = α × exp(β sqrt(triali))

initial <- c(a = 5, b = log(10/5))

summary(m0h <- nls(CorrectN ~ a*exp(b*sqrt(Trial)), data = dtah1,
                  start = initial))

## 
## Formula: CorrectN ~ a * exp(b * sqrt(Trial))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)
## a   8.9294     0.5025   17.77  7.3e-13
## b   0.1771     0.0227    7.81  3.5e-07
## 
## Residual standard error: 0.864 on 18 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 4.6e-07

4.4 observed and fitted

dtah1$yhat1 <- fitted(m0h)

ggplot(dtah1, aes(Trial, CorrectN, color = Group)) +
 geom_point(pch = 1, size = 3) +
geom_smooth(method = "nls", formula = y ~ a*exp(b*sqrt(x)), 
              method.args = list(start = initial), 
              se = F,
              size = rel(.5))+
   scale_x_continuous(limits = c(0, 10), breaks = seq(0, 10, by = 2)) +
  scale_y_continuous(limits = c(5, 20), breaks = seq(5, 20, by = 5)) +
 guides(color = guide_legend(reverse = T)) + 
 labs(x = "Trial on list B", y = "Mean number of correct responses")+
 theme(legend.position = c(.9,.9)) 

4.5 residual plot

plot(m0h)