HW3
1 daya
## 'data.frame': 252 obs. of 3 variables:
## $ ID : Factor w/ 28 levels "S11","S12","S13",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ Trial: int 1 2 3 4 5 6 7 8 9 1 ...
## $ Plane: int 39 40 47 51 52 50 60 52 56 38 ...## 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 502 plot-descriptive
# individual profiles on one panel
library(ggplot2)
ggplot(dta3, aes(Trial, Plane, group = ID)) +
geom_line() +
geom_point() +
# facet_wrap( ~ ID) +
labs(x = "Trial", y = "Plane")
# individual profiles on separate panels
ggplot(dta3, aes(Trial, Plane)) +
geom_line() +
geom_point() +
facet_wrap( ~ ID) +
labs(x = "Trial", y = "Plane")
# coerce data into groupedData object
library(nlme)
dta3_g <- groupedData(Plane ~ Trial | ID, data = as.data.frame(dta3),
labels = list(x = "Trial", y = "number of planes land safely"),
units = list(x = "(order)", y = "(unit)") )
# lattice plot by subject
plot(dta3_g, pch = 1, aspect = .9)
3 nls
dta3_init <- c(a = 1, b = log(10/1))
summary(m0 <- nls(Plane ~ a*exp(b*sqrt(Trial)), data = dta3,
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-074 plot-model fit
ggplot(dta3, aes(Trial, Plane)) +
stat_smooth(method = "nls", formula = y ~ a*exp(b*sqrt(x)),
method.args = list(start = dta3_init),
se = F,
size = rel(.5)) +
geom_point(pch = 1, size = rel(2)) +
scale_x_continuous(limits = c(0, 10), breaks = seq(0, 10, by = 5)) +
scale_y_continuous(limits = c(0, 70), breaks = seq(0, 70, by = 20)) +
facet_wrap( ~ ID) +
labs(x = "Trial (order)",
y = "Number of planes land safely (unit)")
5 Self-Starting Nls Asymptotic Regression Model with an Offset
## 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## 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## 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## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Plane ~ SSmicmen(Trial, Vm, K)
## Data: dta3_g
## AIC BIC logLik
## 1632.314 1646.431 -812.1568
##
## Random effects:
## Formula: Vm ~ 1 | ID
## Vm Residual
## StdDev: 8.988347 5.192305
##
## Fixed effects: list(Vm ~ 1, K ~ 1)
## Value Std.Error DF t-value p-value
## Vm 42.64821 1.9577317 223 21.78450 0
## K 1.22035 0.1196625 223 10.19829 0
## Correlation:
## Vm
## K 0.439
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.79153678 -0.52887678 -0.03144453 0.53400462 3.16157119
##
## Number of Observations: 252
## Number of Groups: 28
6 individual profiles on separate panel
dta3_m1a <- data.frame(dta3, fit = fitted(m1a_g), fit1 = fitted(m1_g))
ggplot(dta3_m1a, aes(Trial, fit)) +
geom_point(aes(Trial, Plane), pch = 1)+
geom_line(col ="steelblue", lwd = .6)+
geom_line(aes(Trial, fit1), col = "indianred", lwd = .6) +
facet_wrap( ~ ID) +
labs(x = "Trial (order)", y = "Number of planes land safely (unit)")
7 individual profiles on one panel
ggplot(dta3_m1a, aes(Trial, fit, group = ID)) +
geom_point(aes(Trial, fit), pch = 1, alpha = .6)+
geom_line(alpha = .3, col="steelblue")+
geom_line(aes(Trial, fit1), col = "indianred", alpha = .3) +
labs(x = "Trial (order)", y = "Number of planes land safely (unit)")
比較兩個fit 的curve