1019_In Class Exercise 2
Pei-Yun,Li
Sun Oct 25 20:18:53 2020
1 Description
In the Netherlands, residents in the rural area of Vlagtwedde in the north-east and residents in the urban, industrial area of Vlaardingen in the south-west, were recruited to participate in a study on chronic obstructive lung diseases. The participants, initially aged 15-44, were surveyed approximately every 3 years for up to 21 years. At each survey, information on respiratory symptoms and smoking status was
collected by questionnaire and spirometry was performed. Pulmonary function was determined by spirometry and a measure of forced expiratory volume (FEV1) was obtained every three years for the first 15 years of the study, and also at year 19. The number of repeated measurements of FEV1 on each subject varied from 1 to 7.
Here the dataset is comprised of a sub-sample of 133 residents aged 36 or older at their entry into the study and whose smoking status did
not change over the 19 years of follow-up. Each participant was either a current or former smoker. Current smoking was defined as smoking at least one cigarette per day.
Column 1: Subject ID
Column 2: Smoker status
Column 3: Time from baseline (years), 0, 3, 6, 9, 12, 15, 19
Column 4: Forced expiratory volume (FEV1)
2 Load packages
pacman::p_load(tidyverse, nlme, car)3 Input data
dta <- read.csv("C:/Users/HANK/Desktop/HOMEWORK/fev_smoker.csv", h = T)str(dta)## 'data.frame': 771 obs. of 4 variables:
## $ ID : chr "P1" "P1" "P1" "P1" ...
## $ Smoker: chr "Former" "Former" "Former" "Former" ...
## $ Time : int 0 3 6 9 15 19 0 3 6 9 ...
## $ FEV1 : num 3.4 3.4 3.45 3.2 2.95 2.4 3.1 3.15 3.5 2.95 ...head(dta)## ID Smoker Time FEV1
## 1 P1 Former 0 3.40
## 2 P1 Former 3 3.40
## 3 P1 Former 6 3.45
## 4 P1 Former 9 3.20
## 5 P1 Former 15 2.95
## 6 P1 Former 19 2.403.1 save a copy of data in wide format
dtaW <- spread(dta, key = Time, value = FEV1)
names(dtaW)[3:9] <- paste0("Year", names(dtaW)[3:9])
head(dtaW)## ID Smoker Year0 Year3 Year6 Year9 Year12 Year15 Year19
## 1 P1 Former 3.40 3.40 3.45 3.20 NA 2.95 2.40
## 2 P10 Current 2.65 NA 3.30 2.40 2.15 2.15 2.25
## 3 P100 Current 3.00 2.75 2.45 NA NA 2.20 2.18
## 4 P101 Current 2.70 3.00 NA NA 2.25 1.85 1.60
## 5 P102 Current 2.80 3.35 NA 2.95 2.50 2.45 NA
## 6 P103 Current NA 3.10 NA 2.75 2.50 1.95 2.654 make ellipses in pairwise plot (by year,by group)
car::scatterplotMatrix(~ Year0 + Year3 + Year6 + Year9 + Year12 + Year15 + Year19,
data = dtaW,
smooth = F,
by.group = TRUE,
regLine = F,
ellipse = T,
diagonal = T,
pch = c(1, 20))
5 means & variances
aggregate(FEV1 ~ Time + Smoker, data = dta, mean, na.rm=T)## Time Smoker FEV1
## 1 0 Current 3.229294
## 2 3 Current 3.119474
## 3 6 Current 3.091461
## 4 9 Current 2.870706
## 5 12 Current 2.797901
## 6 15 Current 2.680959
## 7 19 Current 2.497703
## 8 0 Former 3.519130
## 9 3 Former 3.577778
## 10 6 Former 3.264643
## 11 9 Former 3.171667
## 12 12 Former 3.136207
## 13 15 Former 2.870000
## 14 19 Former 2.907857aggregate(FEV1 ~ Time + Smoker, data = dta, var, na.rm=T)## Time Smoker FEV1
## 1 0 Current 0.3402876
## 2 3 Current 0.3007072
## 3 6 Current 0.3191581
## 4 9 Current 0.3165019
## 5 12 Current 0.2969943
## 6 15 Current 0.2642588
## 7 19 Current 0.2546563
## 8 0 Former 0.2358265
## 9 3 Former 0.3941026
## 10 6 Former 0.3871962
## 11 9 Former 0.3380489
## 12 12 Former 0.3760530
## 13 15 Former 0.3709391
## 14 19 Former 0.41787676 covariance matrices
dtaW %>%
filter( Smoker == "Former") %>%
select( -(1:2)) %>%
var(na.rm = T)## Year0 Year3 Year6 Year9 Year12 Year15 Year19
## Year0 0.3617433 0.3442611 0.3550833 0.2358222 0.2398333 0.3091844 0.2144011
## Year3 0.3442611 0.3566944 0.3479167 0.2434444 0.2455556 0.3014333 0.2141278
## Year6 0.3550833 0.3479167 0.3990278 0.2583333 0.2502778 0.3320000 0.2674722
## Year9 0.2358222 0.2434444 0.2583333 0.1926667 0.1827778 0.2322889 0.1953444
## Year12 0.2398333 0.2455556 0.2502778 0.1827778 0.2022222 0.2201111 0.1406667
## Year15 0.3091844 0.3014333 0.3320000 0.2322889 0.2201111 0.3046711 0.2351044
## Year19 0.2144011 0.2141278 0.2674722 0.1953444 0.1406667 0.2351044 0.3509656dtaW %>%
filter( Smoker == "Current") %>%
select( -(1:2)) %>%
var(na.rm = T)## Year0 Year3 Year6 Year9 Year12 Year15 Year19
## Year0 0.2068623 0.1958723 0.1550242 0.1740866 0.1580108 0.1527199 0.1550177
## Year3 0.1958723 0.2418398 0.1676342 0.1770779 0.1647835 0.1818398 0.1675022
## Year6 0.1550242 0.1676342 0.1618623 0.1328485 0.1484870 0.1479152 0.1571320
## Year9 0.1740866 0.1770779 0.1328485 0.2496970 0.1501407 0.1508160 0.1533355
## Year12 0.1580108 0.1647835 0.1484870 0.1501407 0.1894372 0.1709978 0.1528485
## Year15 0.1527199 0.1818398 0.1479152 0.1508160 0.1709978 0.1873327 0.1654093
## Year19 0.1550177 0.1675022 0.1571320 0.1533355 0.1528485 0.1654093 0.1876729str(dtaW)## 'data.frame': 133 obs. of 9 variables:
## $ ID : chr "P1" "P10" "P100" "P101" ...
## $ Smoker: chr "Former" "Current" "Current" "Current" ...
## $ Year0 : num 3.4 2.65 3 2.7 2.8 NA NA 4.1 2.07 3.25 ...
## $ Year3 : num 3.4 NA 2.75 3 3.35 3.1 4.4 3.65 2.2 3.3 ...
## $ Year6 : num 3.45 3.3 2.45 NA NA NA 3.5 3.9 2.35 3.2 ...
## $ Year9 : num 3.2 2.4 NA NA 2.95 2.75 3.6 2.75 2.25 3.3 ...
## $ Year12: num NA 2.15 NA 2.25 2.5 2.5 NA 3.35 2.2 2.85 ...
## $ Year15: num 2.95 2.15 2.2 1.85 2.45 1.95 3.4 2.95 NA 3.2 ...
## $ Year19: num 2.4 2.25 2.18 1.6 NA 2.65 3.2 3.17 NA 3.1 ...7 relevel factor level to produce consistent legend in the plot below
dtaW <- dtaW %>%
mutate(Smoker = factor(Smoker))
dta <- dta %>%
mutate(Smoker = factor(Smoker))
dtaW$Smoker <- relevel(dtaW$Smoker, ref="Former")
dta$Smoker <- relevel(dta$Smoker, ref="Former")theme_set(theme_bw())8 means with error bars
ggplot(dta, aes(Time, FEV1,
group = Smoker, shape = Smoker, linetype = Smoker)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
scale_shape_manual(values = c(1, 2)) +
labs(x = "Time (years since baseline)", y = "Mean FEV1 (liters)",
linetype = "Smoker", shape = "Smoker") +
theme(legend.justification = c(1, 1), legend.position = c(1, 1),
legend.background = element_rect(fill = "white", color = "black"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
從圖中可以看出，無論是之前有吸菸者，或是最近有吸菸者，兩族群的呼氣量測量值(FEV1)皆是遞減趨勢；唯一不同的是最近有吸菸者所測量到的數據要比之前有吸菸者來的低。
#Fit Linear Model
dta <- mutate(dta, Time_f = as.factor(Time))8.1 m0
summary(m0 <- gls(FEV1 ~ Smoker*Time_f, data = dta))## Generalized least squares fit by REML
## Model: FEV1 ~ Smoker * Time_f
## Data: dta
## AIC BIC logLik
## 1359.281 1428.722 -664.6407
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.519130 0.1171462 30.040506 0.0000
## SmokerCurrent -0.289836 0.1320476 -2.194938 0.0285
## Time_f3 0.058647 0.1594157 0.367889 0.7131
## Time_f6 -0.254488 0.1581008 -1.609653 0.1079
## Time_f9 -0.347464 0.1557060 -2.231537 0.0259
## Time_f12 -0.382924 0.1568667 -2.441076 0.0149
## Time_f15 -0.649130 0.1639349 -3.959684 0.0001
## Time_f19 -0.611273 0.1581008 -3.866351 0.0001
## SmokerCurrent:Time_f3 -0.168468 0.1801366 -0.935222 0.3500
## SmokerCurrent:Time_f6 0.116654 0.1795986 0.649527 0.5162
## SmokerCurrent:Time_f9 -0.011124 0.1779636 -0.062510 0.9502
## SmokerCurrent:Time_f12 -0.048469 0.1794916 -0.270037 0.7872
## SmokerCurrent:Time_f15 0.100795 0.1868469 0.539454 0.5897
## SmokerCurrent:Time_f19 -0.120318 0.1815889 -0.662585 0.5078
##
## Correlation:
## (Intr) SmkrCr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerCurrent -0.887
## Time_f3 -0.735 0.652
## Time_f6 -0.741 0.657 0.544
## Time_f9 -0.752 0.667 0.553 0.557
## Time_f12 -0.747 0.663 0.549 0.553 0.562
## Time_f15 -0.715 0.634 0.525 0.529 0.538 0.534
## Time_f19 -0.741 0.657 0.544 0.549 0.557 0.553 0.529
## SmokerCurrent:Time_f3 0.650 -0.733 -0.885 -0.482 -0.489 -0.486 -0.465 -0.482
## SmokerCurrent:Time_f6 0.652 -0.735 -0.479 -0.880 -0.491 -0.487 -0.466 -0.483
## SmokerCurrent:Time_f9 0.658 -0.742 -0.484 -0.488 -0.875 -0.492 -0.470 -0.488
## SmokerCurrent:Time_f12 0.653 -0.736 -0.480 -0.484 -0.491 -0.874 -0.466 -0.484
## SmokerCurrent:Time_f15 0.627 -0.707 -0.461 -0.465 -0.472 -0.468 -0.877 -0.465
## SmokerCurrent:Time_f19 0.645 -0.727 -0.474 -0.478 -0.485 -0.482 -0.461 -0.871
## SC:T_3 SC:T_6 SC:T_9 SC:T_12 SC:T_15
## SmokerCurrent
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerCurrent:Time_f3
## SmokerCurrent:Time_f6 0.539
## SmokerCurrent:Time_f9 0.544 0.546
## SmokerCurrent:Time_f12 0.539 0.541 0.546
## SmokerCurrent:Time_f15 0.518 0.520 0.524 0.520
## SmokerCurrent:Time_f19 0.533 0.535 0.540 0.535 0.514
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.97158077 -0.67377531 0.05043194 0.66748149 3.13946469
##
## Residual standard error: 0.5618133
## Degrees of freedom: 771 total; 757 residual8.2 m1
summary(m1 <- gls(FEV1 ~ Smoker*Time_f,
weights = varIdent(form = ~ 1 | Time_f*Smoker),
data = dta))## Generalized least squares fit by REML
## Model: FEV1 ~ Smoker * Time_f
## Data: dta
## AIC BIC logLik
## 1377.981 1507.604 -660.9907
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time_f * Smoker
## Parameter estimates:
## 0*Former 3*Former 6*Former 9*Former 15*Former 19*Former 0*Current
## 1.000000 1.292730 1.281354 1.197259 1.254158 1.331145 1.201226
## 3*Current 6*Current 9*Current 12*Current 15*Current 19*Current 12*Former
## 1.129204 1.163337 1.158488 1.122216 1.058562 1.039166 1.262785
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.519130 0.1012590 34.75377 0.0000
## SmokerCurrent -0.289836 0.1194016 -2.42741 0.0154
## Time_f3 0.058647 0.1576382 0.37204 0.7100
## Time_f6 -0.254488 0.1551834 -1.63992 0.1014
## Time_f9 -0.347464 0.1467019 -2.36850 0.0181
## Time_f12 -0.382924 0.1523839 -2.51289 0.0122
## Time_f15 -0.649130 0.1603404 -4.04845 0.0001
## Time_f19 -0.611273 0.1586741 -3.85238 0.0001
## SmokerCurrent:Time_f3 -0.168468 0.1789371 -0.94149 0.3468
## SmokerCurrent:Time_f6 0.116654 0.1779643 0.65549 0.5124
## SmokerCurrent:Time_f9 -0.011124 0.1710216 -0.06505 0.9482
## SmokerCurrent:Time_f12 -0.048469 0.1757578 -0.27577 0.7828
## SmokerCurrent:Time_f15 0.100795 0.1825716 0.55209 0.5811
## SmokerCurrent:Time_f19 -0.120318 0.1806162 -0.66615 0.5055
##
## Correlation:
## (Intr) SmkrCr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerCurrent -0.848
## Time_f3 -0.642 0.545
## Time_f6 -0.653 0.553 0.419
## Time_f9 -0.690 0.585 0.443 0.450
## Time_f12 -0.664 0.564 0.427 0.434 0.459
## Time_f15 -0.632 0.536 0.406 0.412 0.436 0.420
## Time_f19 -0.638 0.541 0.410 0.416 0.440 0.424 0.403
## SmokerCurrent:Time_f3 0.566 -0.667 -0.881 -0.369 -0.391 -0.376 -0.357 -0.361
## SmokerCurrent:Time_f6 0.569 -0.671 -0.365 -0.872 -0.393 -0.378 -0.359 -0.363
## SmokerCurrent:Time_f9 0.592 -0.698 -0.380 -0.386 -0.858 -0.393 -0.374 -0.378
## SmokerCurrent:Time_f12 0.576 -0.679 -0.370 -0.376 -0.398 -0.867 -0.364 -0.368
## SmokerCurrent:Time_f15 0.555 -0.654 -0.356 -0.362 -0.383 -0.369 -0.878 -0.354
## SmokerCurrent:Time_f19 0.561 -0.661 -0.360 -0.366 -0.387 -0.373 -0.354 -0.879
## SC:T_3 SC:T_6 SC:T_9 SC:T_12 SC:T_15
## SmokerCurrent
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerCurrent:Time_f3
## SmokerCurrent:Time_f6 0.448
## SmokerCurrent:Time_f9 0.466 0.468
## SmokerCurrent:Time_f12 0.453 0.456 0.474
## SmokerCurrent:Time_f15 0.436 0.439 0.457 0.444
## SmokerCurrent:Time_f19 0.441 0.444 0.462 0.449 0.432
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.04445773 -0.68139172 0.04873177 0.67419760 2.87621234
##
## Residual standard error: 0.4856209
## Degrees of freedom: 771 total; 757 residual8.3 m2
summary(m2 <- gls(FEV1 ~ Smoker*Time_f,
weights = varExp(form = ~ Time),
data = dta))## Generalized least squares fit by REML
## Model: FEV1 ~ Smoker * Time_f
## Data: dta
## AIC BIC logLik
## 1360.927 1434.996 -664.4633
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~Time
## Parameter estimates:
## expon
## -0.002528783
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.519130 0.1197648 29.383681 0.0000
## SmokerCurrent -0.289836 0.1349993 -2.146947 0.0321
## Time_f3 0.058647 0.1624138 0.361098 0.7181
## Time_f6 -0.254488 0.1605419 -1.585179 0.1133
## Time_f9 -0.347464 0.1576421 -2.204130 0.0278
## Time_f12 -0.382924 0.1582709 -2.419418 0.0158
## Time_f15 -0.649130 0.1645762 -3.944255 0.0001
## Time_f19 -0.611273 0.1582604 -3.862453 0.0001
## SmokerCurrent:Time_f3 -0.168468 0.1835206 -0.917977 0.3589
## SmokerCurrent:Time_f6 0.116654 0.1823482 0.639733 0.5225
## SmokerCurrent:Time_f9 -0.011124 0.1801132 -0.061764 0.9508
## SmokerCurrent:Time_f12 -0.048469 0.1810081 -0.267774 0.7889
## SmokerCurrent:Time_f15 0.100795 0.1874982 0.537580 0.5910
## SmokerCurrent:Time_f19 -0.120318 0.1815964 -0.662558 0.5078
##
## Correlation:
## (Intr) SmkrCr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerCurrent -0.887
## Time_f3 -0.737 0.654
## Time_f6 -0.746 0.662 0.550
## Time_f9 -0.760 0.674 0.560 0.567
## Time_f12 -0.757 0.671 0.558 0.565 0.575
## Time_f15 -0.728 0.646 0.537 0.543 0.553 0.551
## Time_f19 -0.757 0.671 0.558 0.565 0.575 0.573 0.551
## SmokerCurrent:Time_f3 0.653 -0.736 -0.885 -0.487 -0.496 -0.494 -0.475 -0.494
## SmokerCurrent:Time_f6 0.657 -0.740 -0.484 -0.880 -0.499 -0.497 -0.478 -0.497
## SmokerCurrent:Time_f9 0.665 -0.750 -0.490 -0.496 -0.875 -0.503 -0.484 -0.503
## SmokerCurrent:Time_f12 0.662 -0.746 -0.488 -0.494 -0.503 -0.874 -0.481 -0.501
## SmokerCurrent:Time_f15 0.639 -0.720 -0.471 -0.477 -0.485 -0.483 -0.878 -0.483
## SmokerCurrent:Time_f19 0.660 -0.743 -0.486 -0.492 -0.501 -0.499 -0.480 -0.871
## SC:T_3 SC:T_6 SC:T_9 SC:T_12 SC:T_15
## SmokerCurrent
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerCurrent:Time_f3
## SmokerCurrent:Time_f6 0.545
## SmokerCurrent:Time_f9 0.551 0.555
## SmokerCurrent:Time_f12 0.549 0.552 0.559
## SmokerCurrent:Time_f15 0.530 0.533 0.540 0.537
## SmokerCurrent:Time_f19 0.547 0.550 0.557 0.554 0.535
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.92874262 -0.66677628 0.05044755 0.66745312 3.16543499
##
## Residual standard error: 0.5743718
## Degrees of freedom: 771 total; 757 residual8.4 m3
summary(m3 <- gls(FEV1 ~ Smoker*Time_f,
weights = varExp(form = ~ Time | Smoker),
data = dta))## Generalized least squares fit by REML
## Model: FEV1 ~ Smoker * Time_f
## Data: dta
## AIC BIC logLik
## 1358.216 1436.915 -662.1079
##
## Variance function:
## Structure: Exponential of variance covariate, different strata
## Formula: ~Time | Smoker
## Parameter estimates:
## Former Current
## 0.005425870 -0.006170002
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.519130 0.1201359 29.292914 0.0000
## SmokerCurrent -0.289836 0.1354176 -2.140315 0.0326
## Time_f3 0.058647 0.1647238 0.356034 0.7219
## Time_f6 -0.254488 0.1645770 -1.546313 0.1224
## Time_f9 -0.347464 0.1631958 -2.129122 0.0336
## Time_f12 -0.382924 0.1657444 -2.310326 0.0211
## Time_f15 -0.649130 0.1752396 -3.704246 0.0002
## Time_f19 -0.611273 0.1703016 -3.589357 0.0004
## SmokerCurrent:Time_f3 -0.168468 0.1854898 -0.908232 0.3640
## SmokerCurrent:Time_f6 0.116654 0.1856193 0.628459 0.5299
## SmokerCurrent:Time_f9 -0.011124 0.1844802 -0.060302 0.9519
## SmokerCurrent:Time_f12 -0.048469 0.1868438 -0.259411 0.7954
## SmokerCurrent:Time_f15 0.100795 0.1959415 0.514415 0.6071
## SmokerCurrent:Time_f19 -0.120318 0.1909351 -0.630152 0.5288
##
## Correlation:
## (Intr) SmkrCr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerCurrent -0.887
## Time_f3 -0.729 0.647
## Time_f6 -0.730 0.648 0.532
## Time_f9 -0.736 0.653 0.537 0.537
## Time_f12 -0.725 0.643 0.529 0.529 0.534
## Time_f15 -0.686 0.608 0.500 0.500 0.505 0.497
## Time_f19 -0.705 0.626 0.514 0.515 0.519 0.511 0.484
## SmokerCurrent:Time_f3 0.648 -0.730 -0.888 -0.473 -0.477 -0.469 -0.444 -0.457
## SmokerCurrent:Time_f6 0.647 -0.730 -0.472 -0.887 -0.476 -0.469 -0.444 -0.457
## SmokerCurrent:Time_f9 0.651 -0.734 -0.475 -0.475 -0.885 -0.472 -0.446 -0.459
## SmokerCurrent:Time_f12 0.643 -0.725 -0.469 -0.469 -0.473 -0.887 -0.441 -0.454
## SmokerCurrent:Time_f15 0.613 -0.691 -0.447 -0.448 -0.451 -0.444 -0.894 -0.433
## SmokerCurrent:Time_f19 0.629 -0.709 -0.459 -0.459 -0.463 -0.456 -0.431 -0.892
## SC:T_3 SC:T_6 SC:T_9 SC:T_12 SC:T_15
## SmokerCurrent
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerCurrent:Time_f3
## SmokerCurrent:Time_f6 0.533
## SmokerCurrent:Time_f9 0.536 0.536
## SmokerCurrent:Time_f12 0.529 0.529 0.532
## SmokerCurrent:Time_f15 0.505 0.504 0.507 0.501
## SmokerCurrent:Time_f19 0.518 0.517 0.521 0.514 0.490
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.95176426 -0.67527408 0.04683313 0.67280125 2.86836115
##
## Residual standard error: 0.5761515
## Degrees of freedom: 771 total; 757 residual8.5 m4
summary(m4 <- gls(FEV1 ~ Smoker*Time_f,
correlation = corSymm(form = ~ 1 | ID),
weights = varIdent(form = ~ 1 | Time_f*Smoker),
data = dta))## Generalized least squares fit by REML
## Model: FEV1 ~ Smoker * Time_f
## Data: dta
## AIC BIC logLik
## 361.1039 587.9427 -131.5519
##
## Correlation Structure: General
## Formula: ~1 | ID
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5 6
## 2 0.852
## 3 0.856 0.888
## 4 0.849 0.846 0.894
## 5 0.841 0.853 0.891 0.892
## 6 0.861 0.866 0.896 0.872 0.935
## 7 0.753 0.756 0.844 0.789 0.789 0.866
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time_f * Smoker
## Parameter estimates:
## 0*Former 3*Former 6*Former 9*Former 15*Former 19*Former 0*Current
## 1.0000000 1.1548253 1.0981217 1.0095093 1.1970028 1.0796092 1.1241763
## 6*Current 9*Current 12*Current 15*Current 19*Current 3*Current 12*Former
## 1.0518756 1.1013750 1.0380913 1.0158861 0.9676265 1.0363118 1.0527522
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.484680 0.09850685 35.37500 0.0000
## SmokerCurrent -0.257911 0.11599433 -2.22348 0.0265
## Time_f3 -0.024095 0.06609965 -0.36453 0.7156
## Time_f6 -0.221219 0.06278074 -3.52368 0.0005
## Time_f9 -0.328856 0.05936269 -5.53977 0.0000
## Time_f12 -0.394901 0.06223303 -6.34552 0.0000
## Time_f15 -0.479012 0.06777450 -7.06774 0.0000
## Time_f19 -0.598701 0.06580692 -9.09784 0.0000
## SmokerCurrent:Time_f3 -0.065321 0.07460797 -0.87552 0.3816
## SmokerCurrent:Time_f6 0.063890 0.07168101 0.89130 0.3730
## SmokerCurrent:Time_f9 -0.052926 0.06933374 -0.76336 0.4455
## SmokerCurrent:Time_f12 -0.014889 0.07199623 -0.20680 0.8362
## SmokerCurrent:Time_f15 -0.103704 0.07642577 -1.35692 0.1752
## SmokerCurrent:Time_f19 -0.127043 0.07531976 -1.68672 0.0921
##
## Correlation:
## (Intr) SmkrCr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerCurrent -0.849
## Time_f3 -0.129 0.110
## Time_f6 -0.197 0.168 0.608
## Time_f9 -0.334 0.283 0.521 0.662
## Time_f12 -0.289 0.245 0.524 0.647 0.700
## Time_f15 -0.062 0.053 0.515 0.588 0.572 0.719
## Time_f19 -0.258 0.219 0.493 0.627 0.613 0.653 0.690
## SmokerCurrent:Time_f3 0.115 -0.204 -0.886 -0.539 -0.462 -0.465 -0.456 -0.437
## SmokerCurrent:Time_f6 0.173 -0.249 -0.533 -0.876 -0.579 -0.566 -0.515 -0.549
## SmokerCurrent:Time_f9 0.286 -0.332 -0.446 -0.566 -0.856 -0.600 -0.490 -0.524
## SmokerCurrent:Time_f12 0.250 -0.326 -0.453 -0.559 -0.605 -0.864 -0.622 -0.565
## SmokerCurrent:Time_f15 0.055 -0.154 -0.457 -0.522 -0.507 -0.638 -0.887 -0.612
## SmokerCurrent:Time_f19 0.225 -0.322 -0.431 -0.548 -0.535 -0.571 -0.603 -0.874
## SC:T_3 SC:T_6 SC:T_9 SC:T_12 SC:T_15
## SmokerCurrent
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerCurrent:Time_f3
## SmokerCurrent:Time_f6 0.620
## SmokerCurrent:Time_f9 0.528 0.661
## SmokerCurrent:Time_f12 0.539 0.650 0.692
## SmokerCurrent:Time_f15 0.532 0.600 0.579 0.724
## SmokerCurrent:Time_f19 0.512 0.632 0.604 0.654 0.701
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.00857813 -0.65224622 0.05368483 0.64660517 3.17724882
##
## Residual standard error: 0.5411956
## Degrees of freedom: 771 total; 757 residual9 comparing models
anova(m0, m1, m2, m3, m4)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 15 1359.2814 1428.7218 -664.6407
## m1 2 28 1377.9814 1507.6035 -660.9907 1 vs 2 7.3000 0.8860
## m2 3 16 1360.9267 1434.9965 -664.4633 2 vs 3 6.9453 0.8612
## m3 4 17 1358.2158 1436.9149 -662.1079 3 vs 4 4.7109 0.0300
## m4 5 49 361.1039 587.9427 -131.5519 4 vs 5 1061.1119 <.00019.1 show estimated covariance
m4 <- update(m4, method="ML")
intervals(m4, which="var-cov")## Approximate 95% confidence intervals
##
## Correlation structure:
## lower est. upper
## cor(1,2) 0.7984104 0.8523042 0.8926498
## cor(1,3) 0.8033256 0.8559242 0.8952727
## cor(1,4) 0.7947990 0.8498786 0.8910711
## cor(1,5) 0.7853424 0.8419941 0.8846604
## cor(1,6) 0.8045082 0.8615122 0.9027897
## cor(1,7) 0.5939791 0.7501549 0.8518415
## cor(2,3) 0.8467963 0.8884937 0.9193380
## cor(2,4) 0.7906006 0.8468000 0.8888527
## cor(2,5) 0.8008437 0.8539526 0.8937347
## cor(2,6) 0.8117837 0.8671055 0.9070014
## cor(2,7) 0.6008394 0.7537021 0.8533860
## cor(3,4) 0.8536486 0.8939198 0.9235665
## cor(3,5) 0.8508037 0.8914253 0.9214563
## cor(3,6) 0.8481543 0.8962899 0.9297465
## cor(3,7) 0.7174308 0.8426634 0.9151364
## cor(4,5) 0.8516577 0.8924253 0.9224583
## cor(4,6) 0.8179564 0.8724600 0.9114402
## cor(4,7) 0.6580127 0.7872924 0.8714891
## cor(5,6) 0.9056436 0.9358770 0.9566434
## cor(5,7) 0.6478585 0.7869055 0.8752062
## cor(6,7) 0.7618984 0.8638875 0.9240729
## attr(,"label")
## [1] "Correlation structure:"
##
## Variance function:
## lower est. upper
## 3*Former 0.9687025 1.1565059 1.380719
## 6*Former 0.9312166 1.0999878 1.299347
## 9*Former 0.8537168 1.0114297 1.198278
## 15*Former 0.9942021 1.1984605 1.444684
## 19*Former 0.9035607 1.0802535 1.291499
## 0*Current 0.9718621 1.1411996 1.340042
## 6*Current 0.9271746 1.0677355 1.229605
## 9*Current 0.9585674 1.1179806 1.303905
## 12*Current 0.9069306 1.0541140 1.225183
## 15*Current 0.8851592 1.0311232 1.201157
## 19*Current 0.8413179 0.9817466 1.145615
## 3*Current 0.9129852 1.0521247 1.212469
## 12*Former 0.8828825 1.0546222 1.259769
## attr(,"label")
## [1] "Variance function:"
##
## Residual standard error:
## lower est. upper
## 0.4530627 0.5310278 0.6224094dta_m4 <- data.frame(dta, fit = fitted(m4), fit0 = fitted(m0))10 use fitted values to get error bars and see if they cover the data points
ggplot(dta_m4, aes(Time, fit, linetype = Smoker, shape = Smoker,
group = Smoker)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(aes(Time, fit0), fun.y = mean, geom = "point", col = "skyblue") +
stat_summary(aes(Time, FEV1), fun.data = mean_se,
geom = "pointrange", col = "firebrick", alpha = .3) +
scale_shape_manual(values = c(1, 2)) +
labs(x = "Time (years since baseline)", y = "FEV1", linetype = "Smoker",
shape = "Smoker") +
theme(legend.position = c(0.9, 0.9),
legend.background = element_rect(fill = "white", color = "black"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
11 residual plot
plot(m4, resid(., type = "pearson") ~ fitted(.) | Smoker,
layout = c(2, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
xlab = "Ftted values", ylab = "Pearson residuals")