Homework_1130_2
yuwenchen
12/6/2020
Exercise 1
In a learning task known as a free-recall list, each subject was presented with a list of words, one at a time, for a few second each. After all the words had been presented, the subject was asked to recall as many of the words as possible. On each successive study trial the order of presenting the words differ.
Control: 2 trials on List A, followed by 10 trials on List B. Experimental: 2 trials on List B, followed by 10 trials on List B.
The mean number of correct responses given on each trial for the two group was the outcome variable.
Source: Schwartz, R.M., & Humphreys, M.S. (1973). List differentiation in part/whole free recall. American Journal of Psychology, 86, 79-88.
- Column 1: Group
- Column 2: Time in half years unit
- Column 3: Percent of details recalled
pacman::p_load(tidyverse, nlme)
dta1<-read.table("freeRecall.asc.txt", h=T)
names(dta1)<-c("Group", "Trial", "p_Correct")
head(dta1)## Group Trial p_Correct
## 1 C 1 7.9
## 2 C 2 10.9
## 3 C 3 11.9
## 4 C 4 13.0
## 5 C 5 14.2
## 6 C 6 14.2Different Groups’ correct rate in each Trial
ot<-theme_set(theme_bw())
ggplot(dta1, aes(factor(Trial), p_Correct, group=Group, color=Group))+
geom_point()+
geom_line()+
labs(x="Trial", y="Percent of details recalled")+
theme(legend.position = c(.9,.2))
ggplot(dta1, aes(Group, p_Correct, group=factor(Trial), color=factor(Trial)))+
geom_point()+
geom_line()+
labs(x="Group", y="Percent of details recalled")
Model
start1<- c(a = 5, b = log(10/5))
summary(m0 <- nls(p_Correct ~ a*exp(b*sqrt(Trial)), data = dta1,
start = start1))##
## Formula: p_Correct ~ a * exp(b * sqrt(Trial))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 8.92939 0.50253 17.769 7.34e-13 ***
## b 0.17714 0.02269 7.806 3.47e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.864 on 18 degrees of freedom
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 4.602e-07Fitted Plot and Residual Plot
ggplot(dta1, aes(Trial, p_Correct, color = Group)) +
stat_smooth(method = "nls", formula = y ~ a*exp(b*sqrt(x)),
method.args = list(start = start1),
se = F,
size = rel(.5),
linetype="dotted") +
geom_point(pch = 1, size = rel(2)) +
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)) +
labs(x = "Trials",
y = "Percent of correct responses")
plot(m0)
Exercise 2
During a 3-week period, subjects played a chess board game called, “Fox & Geese”. Each subject played up to 27 games. The number of moves made by the subject before commiting a catastrophic error (resulting in defeat) was the outcome of interest. The reading scores of these subjects were also collected as a potential covariate.
Source: Tivnan, T. (1980). Dissertation at Harvard Graduate School of Education. Reported in Singer, J.D., & Willet, J.B. (2003). Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence.
- Column 1: Subject ID
- Column 2: Game played
- Column 3: Number of moves made before defeat
- Column 4: Reading score
dta2<-read.table("fox_geese.txt", h=T)
str(dta2)## 'data.frame': 445 obs. of 4 variables:
## $ ID : chr "S101" "S101" "S101" "S101" ...
## $ Game : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Moves: int 4 7 8 3 3 3 7 6 3 7 ...
## $ Read : num 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 ...psych::describe(dta2)## Warning in psych::describe(dta2): 強制變更過程中產生了 NA## Warning in FUN(newX[, i], ...): min 中沒有無漏失的引數; 回傳 Inf## Warning in FUN(newX[, i], ...): max 中沒有無漏失的引數；回傳 -Inf## vars n mean sd median trimmed mad min max range skew kurtosis
## ID* 1 445 NaN NA NA NaN NA Inf -Inf -Inf NA NA
## Game 2 445 14.03 7.79 14.0 14.03 10.38 1 27 26 0.00 -1.21
## Moves 3 445 7.50 5.80 6.0 6.59 4.45 1 20 19 1.21 0.10
## Read 4 445 1.96 0.78 1.8 1.84 0.44 1 4 3 1.41 1.03
## se
## ID* NA
## Game 0.37
## Moves 0.27
## Read 0.04Individual profiles in separate panels
ggplot(dta2, aes(Game, Moves))+
geom_point()+
geom_line()+
facet_wrap(~ID)+
labs(x="Game", "Moves")
Individual profiles in one panel
ggplot(dta2, aes(Game, Moves, group=ID))+
geom_point(alpha=.3)+
geom_line(alpha=.5)+
labs(x="Game", "Moves")
dta2_g <- groupedData(Moves ~ Game | ID, outer = ~ Read,
data = as.data.frame(dta2),
labels = list(x = "Game", y = "Moves"),
units = list(x = "", y = "") )
plot(dta2_g, pch = 1, aspect = .9)
mfx <- function(x, A, B) { (1 + (19 /(1 + (A*exp(-B*x)))))}
m0 <- nlsList(Moves ~ mfx(Game, A, B),
start = c(A = 150, B = 0.1),
data = dta2_g)
coef(m0)## A B
## S113 3.444376 0.019393820
## S110 19.021257 0.083381612
## S101 3.085531 -0.034697912
## S102 28.499569 0.148899146
## S103 5.385205 0.005747677
## S104 6.257776 0.133481440
## S108 26.681234 0.250404162
## S116 7.093699 0.141653735
## S106 8.742240 0.017753536
## S105 16.633832 0.125658158
## S117 NA NA
## S115 16.674963 0.160175972
## S107 3.146240 0.047900002
## S114 101.147194 0.298797443
## S112 88.215308 0.176818526
## S111 56.277737 0.198327022
## S109 8.227554 0.155514017m1 <- nlme(m0, random=pdDiag(A + B ~ 1))
summary(m1)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Moves ~ mfx(Game, A, B)
## Data: dta2_g
## AIC BIC logLik
## 2500.012 2520.502 -1245.006
##
## Random effects:
## Formula: list(A ~ 1, B ~ 1)
## Level: ID
## Structure: Diagonal
## A B Residual
## StdDev: 2.62854 0.04992857 3.745157
##
## Fixed effects: list(A ~ 1, B ~ 1)
## Value Std.Error DF t-value p-value
## A 10.357095 1.592630 427 6.503138 0
## B 0.109389 0.014638 427 7.472905 0
## Correlation:
## A
## B 0.472
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.5723172 -0.5887838 -0.1064716 0.5565760 3.4574177
##
## Number of Observations: 445
## Number of Groups: 17#重新取得起始值
start2 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)
start2## A B
## 10.357 0.109#將ID的變異數加權
m2 <- update(m1, weights=varIdent(form = ~ 1 | ID),
control=nlmeControl(maxIter=1e8, msMaxIter=1e8, niterEM=500))
summary(m2)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Moves ~ mfx(Game, A, B)
## Data: dta2_g
## AIC BIC logLik
## 2456.174 2542.234 -1207.087
##
## Random effects:
## Formula: list(A ~ 1, B ~ 1)
## Level: ID
## Structure: Diagonal
## A B Residual
## StdDev: 3.468549 0.05936852 3.982058
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | ID
## Parameter estimates:
## S113 S110 S101 S102 S103 S104 S108 S116
## 1.0000000 0.8277692 0.4279192 1.1215884 0.5515779 1.1217123 0.7118693 1.2247191
## S106 S105 S117 S115 S107 S114 S112 S111
## 0.3392784 0.8255142 1.0328388 0.9281142 0.6649246 0.9527534 0.9995646 1.1471665
## S109
## 1.2466191
## Fixed effects: list(A ~ 1, B ~ 1)
## Value Std.Error DF t-value p-value
## A 8.579726 1.2977630 427 6.611165 0
## B 0.096326 0.0159375 427 6.044000 0
## Correlation:
## A
## B 0.287
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.1980904 -0.6733035 -0.1749345 0.5312257 4.0182558
##
## Number of Observations: 445
## Number of Groups: 17#新的起始值
m3 <- update(m2,
fixed = list(A + B ~ Read),
random = pdDiag(A + B ~ 1),
groups = ~ ID,
start = list(fixed = c(A = c(start2[1], 0),
B = c(start2[2], 0))))
summary(m3)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Moves ~ mfx(Game, A, B)
## Data: dta2_g
## AIC BIC logLik
## 2452.096 2546.352 -1203.048
##
## Random effects:
## Formula: list(A ~ 1, B ~ 1)
## Level: ID
## Structure: Diagonal
## A.(Intercept) B.(Intercept) Residual
## StdDev: 2.026809 0.05351836 3.894972
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | ID
## Parameter estimates:
## S113 S110 S101 S102 S103 S104 S108 S116
## 1.0000000 0.8633942 0.4407039 1.1670976 0.5603939 1.1435756 0.7424579 1.2491535
## S106 S105 S117 S115 S107 S114 S112 S111
## 0.3458038 0.8407412 1.0628410 0.9404114 0.8148872 0.9537728 0.9889875 1.1371785
## S109
## 1.3115072
## Fixed effects: list(A + B ~ Read)
## Value Std.Error DF t-value p-value
## A.(Intercept) -4.631554 4.837132 425 -0.9575002 0.3389
## A.Read 8.021220 3.152021 425 2.5447863 0.0113
## B.(Intercept) -0.001241 0.040061 425 -0.0309764 0.9753
## B.Read 0.054672 0.019629 425 2.7852689 0.0056
## Correlation:
## A.(In) A.Read B.(In)
## A.Read -0.975
## B.(Intercept) 0.306 -0.273
## B.Read -0.395 0.389 -0.928
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.2978838 -0.6328087 -0.1662573 0.5546342 4.0940870
##
## Number of Observations: 445
## Number of Groups: 17#將相關性考慮進去
m4 <- update(m3, correlation=corAR1(form = ~ 1 | ID))
summary(m4)$coef## $fixed
## A.(Intercept) A.Read B.(Intercept) B.Read
## -2.74734519 7.22295204 0.01193236 0.04959426
##
## $random
## $random$ID
## A.(Intercept) B.(Intercept)
## S113 -4.370989e-09 -0.0270401920
## S110 6.599619e-07 -0.0337840225
## S101 -6.034670e-07 -0.0648744033
## S102 2.580311e-07 -0.0009086408
## S103 -2.790867e-08 -0.0540501711
## S104 -3.402358e-07 0.0525318046
## S108 -2.708467e-08 0.0778790958
## S116 -3.300954e-07 0.0486220769
## S106 8.861005e-08 -0.0727530324
## S105 2.053069e-07 0.0035211144
## S117 2.956522e-07 0.0494358038
## S115 -1.688804e-08 0.0296559422
## S107 -3.239044e-07 0.0029401455
## S114 3.511984e-08 0.0487637634
## S112 1.399938e-07 -0.0476461186
## S111 5.596740e-08 -0.0267635663
## S109 -6.468817e-08 0.0144704006anova(m1, m2,m3,m4)## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 5 2500.012 2520.502 -1245.006
## m2 2 21 2456.174 2542.234 -1207.087 1 vs 2 75.83759 <.0001
## m3 3 23 2452.096 2546.352 -1203.048 2 vs 3 8.07781 0.0176
## m4 4 24 2423.659 2522.013 -1187.830 3 vs 4 30.43707 <.0001Plot Fitted values and Residuals
dta2$yhat <- fitted(m4)
ggplot(dta2, aes(Game, Moves, group=Read, color=Read )) +
geom_point(size = rel(.5), pch = 1) +
geom_line(aes(y = yhat, group = ID), col = "indianred") +
labs(x = "Game played",
y = "Number of moves made before defeat (steps)") +
facet_wrap(~ ID, ncol = 6)+
theme_bw()+
theme(legend.position = "NONE")
plot(m4,id = .05, cex = .6, adj=-.05)
# random effects
plot(ranef(m4, augFrame = T))