1 Data Management

dtaHW3 <- read.table("C:/Users/ASUS/Desktop/data/planes.txt", h=T)
head(dtaHW3)
##    ID Trial Plane
## 1 S11     1    39
## 2 S11     2    40
## 3 S11     3    47
## 4 S11     4    51
## 5 S11     5    52
## 6 S11     6    50

2 Background info. Graph

2.1 All subjects

with the increase in trials, more planes were brought in to land safely

library(ggplot2)
ggplot(dtaHW3, aes(Trial, Plane, group = ID)) + 
  geom_line() +
  geom_point() +
  # facet_wrap( ~ ID) +
  labs(x = "Trial", y = "Plane")

2.2 each subject

ggplot(dtaHW3, aes(Trial, Plane)) + 
  geom_line() +
  geom_point() +
  facet_wrap( ~ ID) +
  labs(x = "Trial", y = "Plane")

3 non-linear model

3.1 nls

給定起始值的非線性模型以nls()做分析

dta3_init <- c(a = 1, b = log(10/1))

summary(m0 <- nls(Plane ~ a*exp(b*sqrt(Trial)), data = dtaHW3,
                  start = dta3_init))
## 
## Formula: Plane ~ a * exp(b * sqrt(Trial))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 16.35533    1.17486   13.92   <2e-16 ***
## b  0.30581    0.02953   10.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.637 on 250 degrees of freedom
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 4.348e-07

3.2 nlme

自動給定起始值的非線性模型以nlme()做分析

library (nlme)
dta3_g <- groupedData(Plane ~ Trial | ID, data = as.data.frame(dtaHW3),
                      labels = list(x = "Trial", y = "number of planes land safely"),
                      units = list(x = "(order)", y = "(unit)") )
m0_g <- nlsList(SSmicmen, data = dta3_g)
## Warning: 2 times caught the same error in nls(y ~ x/(K + x), data = xy,
## start = list(K = abs(pars[2L]/pars[1L])), algorithm = "plinear"): step factor
## 0.000488281 reduced below 'minFactor' of 0.000976562
m0_g
## Call:
##   Model: Plane ~ SSmicmen(Trial, Vm, K) | ID 
##    Data: dta3_g 
## 
## Coefficients:
##            Vm          K
## S37        NA         NA
## S32  38.42522  2.6382919
## S21  35.58824  0.6117632
## S26  33.49397  0.8342775
## S31  38.25314  1.4534320
## S29  46.66226  2.9971659
## S35        NA         NA
## S23  35.37661  0.6846419
## S27  39.04922  0.8537614
## S28  46.16958  2.2907550
## S33  58.51689  4.8733793
## S34  58.73462  5.3899245
## S36  72.57835  7.9155112
## S14  30.51662 -0.1808489
## S30  48.34655  1.5879277
## S19  41.70994  0.9235549
## S25  46.88423  1.9489698
## S20  45.17192  1.2972892
## S38 102.88531 12.2902544
## S18  45.74841  0.8423952
## S22  45.58569  0.9419266
## S13  43.09619  0.1979081
## S24  70.15462  5.2446174
## S15  46.70072  0.4936488
## S17  52.95904  1.5766528
## S12  53.66601  0.7619322
## S16  54.07300  0.7409107
## S11  58.77729  0.6510748
## 
## Degrees of freedom: 234 total; 182 residual
## Residual standard error: 3.575345
m1_g <- nlme(m0_g, control=nlmeControl(maxIter=1e6, msMaxIter=1e6, niterEM=500))

summary(m1_g)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: Plane ~ SSmicmen(Trial, Vm, K) 
##  Data: dta3_g 
##        AIC     BIC    logLik
##   1554.784 1575.96 -771.3919
## 
## Random effects:
##  Formula: list(Vm ~ 1, K ~ 1)
##  Level: ID
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev   Corr 
## Vm       7.582077 Vm   
## K        1.681533 0.441
## Residual 3.861442      
## 
## Fixed effects: list(Vm ~ 1, K ~ 1) 
##       Value Std.Error  DF   t-value p-value
## Vm 45.83635 1.7120608 223 26.772621       0
## K   1.93060 0.3545793 223  5.444764       0
##  Correlation: 
##   Vm   
## K 0.546
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.72172107 -0.66635039 -0.03338046  0.57852813  2.08310019 
## 
## Number of Observations: 252
## Number of Groups: 28