1 input data

dta2 <- read.table("fox_geese.txt", h=T)

dta2$ID <- as.factor(dta2$ID)

str(dta2)
## 'data.frame':    445 obs. of  4 variables:
##  $ ID   : Factor w/ 17 levels "S101","S102",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ 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 ...
head(dta2)
##     ID Game Moves Read
## 1 S101    1     4  1.4
## 2 S101    2     7  1.4
## 3 S101    3     8  1.4
## 4 S101    4     3  1.4
## 5 S101    5     3  1.4
## 6 S101    6     3  1.4

2 individual profiles on separate panels

library(ggplot2)
ggplot(dta2, aes(Game , Moves )) +
 geom_point(size = rel(.8), pch = 20)+
 geom_line()+
 facet_wrap( ~ ID) +
 labs(x = "Number of game played", y = "Number of moves made before defeat")

3 individual profiles on one panel

ggplot(dta2, aes(Game , Moves , group = ID)) +
 geom_point(pch = 20)+
 geom_line(alpha = .3) +
 labs(x = "Number of game played", y = "Number of moves made before defeat")

4 coerce data into groupedData object

library(nlme)
dta_g <- groupedData(Moves ~ Game | ID, data = as.data.frame(dta2),
                     labels = list(x = "Number of game played", 
                                    y = "Number of moves made before defeat"),
                      units = list(x = "(bouts)", y = "(steps)"))

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

5 define the 3-parameter curve fitting function

fx <- function(x, b0, b1) {
               1 + 19 /(1 + b0*exp(-b1*x))}

6 operating characteristic of b1?

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dd <- expand.grid(b1=1:10, x=seq(0, 1.5, by=0.1)) %>% 
  rowwise %>%
  mutate(y=fx(b0=20,
              b1=b1,  
              x=x))

ggplot(dd, 
       aes(x, 
           y, 
           group=b1, 
           color=factor(b1))) +
  geom_path() +   
  labs(x="Number of game played(log)", 
       y="Number of moves made before defeat ") +
  scale_color_grey() +
  theme_minimal()+
  theme(legend.position="NONE")

7 searching for starting values

#b0與b1的seq不知道如何設才是make sense
v0 <- expand.grid(list(b0=seq(1,150, by=1), b1=seq(0, 1, by=0.1)))
res_v0 <- nls2::nls2(Moves ~ fx(Game, b0, b1), data=subset(dta_g, ID=="S113"), start=v0,
                     algorithm="brute-force")
coef(res_v0)
## b0 b1 
##  3  0
mv0 <- nls(Moves ~ fx(Game, b0, b1), data=subset(dta_g, ID=="S113"), 
          start=coef(res_v0),control = list(maxiter = 500))
coef(mv0)
##         b0         b1 
## 3.44442216 0.01939469

8 nlsList

m0 <- nlsList(Moves ~ fx(Game, b0, b1), start=coef(mv0), data=dta_g)
## Warning: 1 error caught in nls(model, data = data, control = controlvals, start
## = start): number of iterations exceeded maximum of 50
coef(m0)
##              b0           b1
## S101   3.085535 -0.034697889
## S106   8.742300  0.017753960
## S103   5.385170  0.005747236
## S107   3.146234  0.047899854
## S102  28.499604  0.148899204
## S104   6.257682  0.133480447
## S105  16.633842  0.125658185
## S108  26.681191  0.250404054
## S109   8.227605  0.155514511
## S110  19.019834  0.083378090
## S111  56.278986  0.198328088
## S112  88.216608  0.176819183
## S113   3.444422  0.019394694
## S114 101.144893  0.298796006
## S115  16.675123  0.160176422
## S116   7.093723  0.141654034
## S117         NA           NA

S117的b0與b1和其他人的差太多了

9 add random term in model

m1 <- nlme(m0, random=pdDiag(b0 + b1 ~ 1))
prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)

prm0
##     b0     b1 
## 10.357  0.109
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 ~ fx(Game, b0, b1) 
##  Data: dta_g 
##        AIC      BIC    logLik
##   2456.174 2542.233 -1207.087
## 
## Random effects:
##  Formula: list(b0 ~ 1, b1 ~ 1)
##  Level: ID
##  Structure: Diagonal
##               b0         b1 Residual
## StdDev: 3.468563 0.05936851 1.704016
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | ID 
##  Parameter estimates:
##      S101      S106      S103      S107      S102      S104      S105      S108 
## 1.0000000 0.7928477 1.2889613 1.5538390 2.6209915 2.6213001 1.9291149 1.6635381 
##      S109      S110      S111      S112      S113      S114      S115      S116 
## 2.9131894 1.9343810 2.6807712 2.3358371 2.3368856 2.2264429 2.1688762 2.8620167 
##      S117 
## 2.4135995 
## Fixed effects: list(b0 ~ 1, b1 ~ 1) 
##       Value Std.Error  DF  t-value p-value
## b0 8.579780 1.2977712 427 6.611166       0
## b1 0.096326 0.0159375 427 6.044021       0
##  Correlation: 
##    b0   
## b1 0.287
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1981094 -0.6733040 -0.1749361  0.5312284  4.0182636 
## 
## Number of Observations: 445
## Number of Groups: 17
m3 <- update(m2, 
             fixed = list(b0 + b1 ~ Read),
             random = pdDiag(b0+ b1 ~ 1),
             groups = ~ ID, 
             start = list(fixed = c(A = c(prm0[1], 0), 
                                    B = c(prm0[2], 0))))
summary(m3)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: Moves ~ fx(Game, b0, b1) 
##  Data: dta_g 
##        AIC      BIC    logLik
##   2452.096 2546.352 -1203.048
## 
## Random effects:
##  Formula: list(b0 ~ 1, b1 ~ 1)
##  Level: ID
##  Structure: Diagonal
##         b0.(Intercept) b1.(Intercept) Residual
## StdDev:       2.026821     0.05351841 1.716477
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | ID 
##  Parameter estimates:
##      S101      S106      S103      S107      S102      S104      S105      S108 
## 1.0000000 0.7846877 1.2716289 1.8490997 2.6483453 2.5949663 1.9077857 1.6847646 
##      S109      S110      S111      S112      S113      S114      S115      S116 
## 2.9760207 1.9591906 2.5804538 2.2441820 2.2691635 2.1642777 2.1339569 2.8345455 
##      S117 
## 2.4117818 
## Fixed effects: list(b0 + b1 ~ Read) 
##                    Value Std.Error  DF    t-value p-value
## b0.(Intercept) -4.631539  4.837153 425 -0.9574929  0.3389
## b0.Read         8.021226  3.152035 425  2.5447771  0.0113
## b1.(Intercept) -0.001241  0.040061 425 -0.0309762  0.9753
## b1.Read         0.054672  0.019629 425  2.7852686  0.0056
##  Correlation: 
##                b0.(I) b0.Red b1.(I)
## b0.Read        -0.975              
## b1.(Intercept)  0.306 -0.273       
## b1.Read        -0.395  0.389 -0.928
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2978743 -0.6328282 -0.1662579  0.5546352  4.0940724 
## 
## Number of Observations: 445
## Number of Groups: 17
m4 <- update(m3, correlation=corAR1(form = ~ 1 | ID))

summary(m4)$coef
## $fixed
## b0.(Intercept)        b0.Read b1.(Intercept)        b1.Read 
##    -2.74769356     7.22310620     0.01193027     0.04959508 
## 
## $random
## $random$ID
##      b0.(Intercept) b1.(Intercept)
## S101  -5.649274e-07  -0.0648747333
## S106   8.294107e-08  -0.0727527961
## S103  -2.611597e-08  -0.0540500659
## S107  -3.031676e-07   0.0029403976
## S102   2.415260e-07  -0.0009085568
## S104  -3.184562e-07   0.0525319134
## S105   1.921668e-07   0.0035213222
## S108  -2.534033e-08   0.0778792407
## S109  -6.054672e-08   0.0144701654
## S110   6.177424e-07  -0.0337843473
## S111   5.238442e-08  -0.0267638159
## S112   1.310308e-07  -0.0476462762
## S113  -4.073225e-09  -0.0270411817
## S114   3.287551e-08   0.0487641053
## S115  -1.580596e-08   0.0296561981
## S116  -3.089659e-07   0.0486223934
## S117   2.767322e-07   0.0494360370
m5 <- update(m4, random=pdDiag(b0 ~ 1))

10 model comparison

anova(m1, m2, m3, m4, m5)
##    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.233 -1207.087 1 vs 2 75.83771  <.0001
## m3     3 23 2452.096 2546.352 -1203.048 2 vs 3  8.07770  0.0176
## m4     4 24 2423.659 2522.013 -1187.830 3 vs 4 30.43710  <.0001
## m5     5 23 2494.401 2588.656 -1224.200 4 vs 5 72.74152  <.0001

比較各模型，都有顯著，選擇AIC最小

11 plot-model fit

dta2$yhat <- fitted(m4)

ggplot(dta2, aes(Game, Moves), group = Read, coloer = Read ) +
  geom_point(size = rel(.5), pch = 1, col = "skyblue") +
  geom_line(aes(y = yhat, group = ID), col = "indianred") +
  labs(x = "Number of game played (bouts)", 
      y = "Number of moves made before defeat (steps)") +
  facet_wrap(~ ID, ncol = 6)+
  theme(legend.position = "NONE")