Week7 In-class exercises 3
Ching-Fang Wu
Mon Oct 26 10:59:55 2020
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.
Source: van der Lende, R., Kok, T.J., Peset, R., Quanjer, P.H., Schouten, J.P. & Orie, N.G.M. (1981). Decreases in VC and FEV1 with time: Indicators for effects of smoking and air pollution. Bulletin of European Physiopathology and Respiration, 17, 775-792.
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)
1 Data management
pacman::p_load(tidyverse, nlme, car)# input data
dta <- read.csv("C:/Users/Ching-Fang Wu/Documents/data/fev_smoker.csv", h = T)# show structure of data
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 ...# show first 8 line
head(dta,8)## 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.40
## 7 P2 Current 0 3.10
## 8 P2 Current 3 3.15每個參與者依據不同Time，重複測量FEV1，總計六次。
#save a copy of data in wide format
dtaW <- dta %>%
spread(key = Time, value = FEV1) #key is the column whose values will become variable names；value is the column where values will fill in under the new variables created from key.names(dtaW)[3:9] <- paste0("Year", names(dtaW)[3:9])paste將字串連接起來，如「Year 0」。
若不要有分隔符號，則用paste0，如「Year0」。
2 Data visualization
2.1 make ellipses in pairwise plot
scatterplotMatrix(~ Year0 + Year3 + Year6 + Year9 + Year12 + Year15 + Year19 | Smoker, #依照smoker分群
data = dtaW, smooth = F, by.group = TRUE,
reg.line = F, ellipse = T, diag = "none", pch = c(1, 20),
col = c("firebrick", "forestgreen")) ## Warning in applyDefaults(diagonal, defaults = list(method =
## "adaptiveDensity"), : unnamed diag arguments, will be ignored## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in plot.window(...): "reg.line" 不是一個繪圖參數## Warning in plot.xy(xy, type, ...): "reg.line" 不是一個繪圖參數## Warning in title(...): "reg.line" 不是一個繪圖參數## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數
## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" 不是一個
## 繪圖參數
3 Analysis
3.1 means & variances
aggregate(FEV1 ~ Time + Smoker, data = dta, mean, na.rm=T) #Time + Smoker為分群因子，計算FEV1的平均數## 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的變異數## 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.41787673.2 covariance matrices
#covariance matrices of "Former"
dtaW %>%
filter( Smoker == "Former") %>% #選要分析的觀察值，觀察值子集 (Row)
select( -(1:2)) %>% #選要分析的欄位，欄位子集 (Column)，-(1:2)表示剔除前兩欄：ID、Smoker
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.3509656#covariance matrices of "Current"
dtaW %>%
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 ...# relevel factor level to produce consistent legend in the plot below
#dta$Smoker <- relevel(dta$Smoker, ref="Former")4
theme_set(theme_bw())5 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.
任一時間點，Current的FEV1比former低，此外，這兩類吸菸者的FEV1隨時間遞減。
6
dta <- mutate(dta, Time_f = as.factor(Time)) #新增Time_f作為時間類別變數7 Modeling
7.1 m0
summary(m0 <- gls(FEV1 ~ Smoker*Time_f, data = dta)) #無權重的gls就等於lm## 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.229294 0.06093721 52.99380 0.0000
## SmokerFormer 0.289836 0.13204761 2.19494 0.0285
## Time_f3 -0.109820 0.08387973 -1.30926 0.1908
## Time_f6 -0.137833 0.08520443 -1.61768 0.1061
## Time_f9 -0.358588 0.08617823 -4.16101 0.0000
## Time_f12 -0.431393 0.08723567 -4.94514 0.0000
## Time_f15 -0.548335 0.08964988 -6.11641 0.0000
## Time_f19 -0.731591 0.08932341 -8.19037 0.0000
## SmokerFormer:Time_f3 0.168468 0.18013659 0.93522 0.3500
## SmokerFormer:Time_f6 -0.116654 0.17959864 -0.64953 0.5162
## SmokerFormer:Time_f9 0.011124 0.17796363 0.06251 0.9502
## SmokerFormer:Time_f12 0.048469 0.17949158 0.27004 0.7872
## SmokerFormer:Time_f15 -0.100795 0.18684687 -0.53945 0.5897
## SmokerFormer:Time_f19 0.120318 0.18158895 0.66259 0.5078
##
## Correlation:
## (Intr) SmkrFr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerFormer -0.461
## Time_f3 -0.726 0.335
## Time_f6 -0.715 0.330 0.520
## Time_f9 -0.707 0.326 0.514 0.506
## Time_f12 -0.699 0.322 0.507 0.500 0.494
## Time_f15 -0.680 0.314 0.494 0.486 0.481 0.475
## Time_f19 -0.682 0.315 0.496 0.488 0.482 0.477 0.464
## SmokerFormer:Time_f3 0.338 -0.733 -0.466 -0.242 -0.239 -0.236 -0.230 -0.231
## SmokerFormer:Time_f6 0.339 -0.735 -0.246 -0.474 -0.240 -0.237 -0.231 -0.231
## SmokerFormer:Time_f9 0.342 -0.742 -0.249 -0.245 -0.484 -0.239 -0.233 -0.234
## SmokerFormer:Time_f12 0.339 -0.736 -0.247 -0.243 -0.240 -0.486 -0.231 -0.232
## SmokerFormer:Time_f15 0.326 -0.707 -0.237 -0.233 -0.231 -0.228 -0.480 -0.222
## SmokerFormer:Time_f19 0.336 -0.727 -0.244 -0.240 -0.237 -0.234 -0.228 -0.492
## SF:T_3 SF:T_6 SF:T_9 SF:T_12 SF:T_15
## SmokerFormer
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerFormer:Time_f3
## SmokerFormer:Time_f6 0.539
## SmokerFormer:Time_f9 0.544 0.546
## SmokerFormer:Time_f12 0.539 0.541 0.546
## SmokerFormer:Time_f15 0.518 0.520 0.524 0.520
## SmokerFormer: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 residual7.2 m1
summary(m1 <- gls(FEV1 ~ Smoker*Time_f,
weights = varIdent(form = ~ 1 | Time_f*Smoker),
data = dta)) #varIdent代表Different variances per stratum，form = ~ 1 | Time_f*Smoker 的意思是Constant variance within stratum.(組內一致，組間有變異)## 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.292710 1.281344 1.197259 1.254151 1.331139 1.201218
## 3*Current 6*Current 9*Current 12*Current 15*Current 19*Current 12*Former
## 1.129204 1.163335 1.158479 1.122197 1.058560 1.039136 1.262762
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.229294 0.06327225 51.03808 0.0000
## SmokerFormer 0.289836 0.11940236 2.42739 0.0154
## Time_f3 -0.109820 0.08466840 -1.29707 0.1950
## Time_f6 -0.137833 0.08711755 -1.58215 0.1140
## Time_f9 -0.358588 0.08790306 -4.07936 0.0000
## Time_f12 -0.431393 0.08757802 -4.92581 0.0000
## Time_f15 -0.548335 0.08731204 -6.28018 0.0000
## Time_f19 -0.731591 0.08628221 -8.47905 0.0000
## SmokerFormer:Time_f3 0.168468 0.17893718 0.94149 0.3468
## SmokerFormer:Time_f6 -0.116654 0.17796497 -0.65549 0.5124
## SmokerFormer:Time_f9 0.011124 0.17102269 0.06505 0.9482
## SmokerFormer:Time_f12 0.048469 0.17575742 0.27577 0.7828
## SmokerFormer:Time_f15 -0.100795 0.18257248 -0.55208 0.5811
## SmokerFormer:Time_f19 0.120318 0.18061661 0.66615 0.5055
##
## Correlation:
## (Intr) SmkrFr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerFormer -0.530
## Time_f3 -0.747 0.396
## Time_f6 -0.726 0.385 0.543
## Time_f9 -0.720 0.381 0.538 0.523
## Time_f12 -0.722 0.383 0.540 0.525 0.520
## Time_f15 -0.725 0.384 0.542 0.526 0.522 0.524
## Time_f19 -0.733 0.389 0.548 0.533 0.528 0.530 0.531
## SmokerFormer:Time_f3 0.354 -0.667 -0.473 -0.257 -0.255 -0.255 -0.256 -0.259
## SmokerFormer:Time_f6 0.356 -0.671 -0.266 -0.490 -0.256 -0.257 -0.258 -0.261
## SmokerFormer:Time_f9 0.370 -0.698 -0.276 -0.269 -0.514 -0.267 -0.268 -0.271
## SmokerFormer:Time_f12 0.360 -0.679 -0.269 -0.261 -0.259 -0.498 -0.261 -0.264
## SmokerFormer:Time_f15 0.347 -0.654 -0.259 -0.252 -0.249 -0.250 -0.478 -0.254
## SmokerFormer:Time_f19 0.350 -0.661 -0.262 -0.254 -0.252 -0.253 -0.254 -0.478
## SF:T_3 SF:T_6 SF:T_9 SF:T_12 SF:T_15
## SmokerFormer
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerFormer:Time_f3
## SmokerFormer:Time_f6 0.448
## SmokerFormer:Time_f9 0.466 0.468
## SmokerFormer:Time_f12 0.453 0.456 0.474
## SmokerFormer:Time_f15 0.436 0.439 0.457 0.444
## SmokerFormer:Time_f19 0.441 0.444 0.462 0.449 0.432
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.04443149 -0.68139604 0.04873136 0.67419703 2.87624105
##
## Residual standard error: 0.4856249
## Degrees of freedom: 771 total; 757 residual7.3 m2
summary(m2 <- gls(FEV1 ~ Smoker*Time_f,
weights = varExp(form = ~ Time), #varExp是Exponential variance structure，varExp(form = ~ Hour)表示變異數隨時間而不同
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.229294 0.06229936 51.83511 0.0000
## SmokerFormer 0.289836 0.13499932 2.14695 0.0321
## Time_f3 -0.109820 0.08544929 -1.28521 0.1991
## Time_f6 -0.137833 0.08647075 -1.59399 0.1114
## Time_f9 -0.358588 0.08711898 -4.11607 0.0000
## Time_f12 -0.431393 0.08783081 -4.91163 0.0000
## Time_f15 -0.548335 0.08983467 -6.10383 0.0000
## Time_f19 -0.731591 0.08905559 -8.21500 0.0000
## SmokerFormer:Time_f3 0.168468 0.18352063 0.91798 0.3589
## SmokerFormer:Time_f6 -0.116654 0.18234823 -0.63973 0.5225
## SmokerFormer:Time_f9 0.011124 0.18011317 0.06176 0.9508
## SmokerFormer:Time_f12 0.048469 0.18100811 0.26777 0.7889
## SmokerFormer:Time_f15 -0.100795 0.18749823 -0.53758 0.5910
## SmokerFormer:Time_f19 0.120318 0.18159639 0.66256 0.5078
##
## Correlation:
## (Intr) SmkrFr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerFormer -0.461
## Time_f3 -0.729 0.336
## Time_f6 -0.720 0.332 0.525
## Time_f9 -0.715 0.330 0.521 0.515
## Time_f12 -0.709 0.327 0.517 0.511 0.507
## Time_f15 -0.693 0.320 0.506 0.500 0.496 0.492
## Time_f19 -0.700 0.323 0.510 0.504 0.500 0.496 0.485
## SmokerFormer:Time_f3 0.339 -0.736 -0.466 -0.245 -0.243 -0.241 -0.235 -0.237
## SmokerFormer:Time_f6 0.342 -0.740 -0.249 -0.474 -0.244 -0.242 -0.237 -0.239
## SmokerFormer:Time_f9 0.346 -0.750 -0.252 -0.249 -0.484 -0.245 -0.240 -0.242
## SmokerFormer:Time_f12 0.344 -0.746 -0.251 -0.248 -0.246 -0.485 -0.239 -0.241
## SmokerFormer:Time_f15 0.332 -0.720 -0.242 -0.239 -0.238 -0.236 -0.479 -0.232
## SmokerFormer:Time_f19 0.343 -0.743 -0.250 -0.247 -0.245 -0.243 -0.238 -0.490
## SF:T_3 SF:T_6 SF:T_9 SF:T_12 SF:T_15
## SmokerFormer
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerFormer:Time_f3
## SmokerFormer:Time_f6 0.545
## SmokerFormer:Time_f9 0.551 0.555
## SmokerFormer:Time_f12 0.549 0.552 0.559
## SmokerFormer:Time_f15 0.530 0.533 0.540 0.537
## SmokerFormer: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 residual1.SmokerFormer對FEV1有顯著影響 2.越遠期的時間類別變數Time_f3、Time_f6到Time_f19對FEV1的影響越小
7.4 m3
summary(m3 <- gls(FEV1 ~ Smoker*Time_f,
weights = varExp(form = ~ Time | Smoker),#Smoker為分群因子，varExp(form = ~ Time | Smoker)表示不同類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.229294 0.06249240 51.67499 0.0000
## SmokerFormer 0.289836 0.13541763 2.14031 0.0326
## Time_f3 -0.109820 0.08527911 -1.28778 0.1982
## Time_f6 -0.137833 0.08584232 -1.60566 0.1088
## Time_f9 -0.358588 0.08602381 -4.16848 0.0000
## Time_f12 -0.431393 0.08625192 -5.00155 0.0000
## Time_f15 -0.548335 0.08765943 -6.25529 0.0000
## Time_f19 -0.731591 0.08633408 -8.47396 0.0000
## SmokerFormer:Time_f3 0.168468 0.18548978 0.90823 0.3640
## SmokerFormer:Time_f6 -0.116654 0.18561925 -0.62846 0.5299
## SmokerFormer:Time_f9 0.011124 0.18448022 0.06030 0.9519
## SmokerFormer:Time_f12 0.048469 0.18684375 0.25941 0.7954
## SmokerFormer:Time_f15 -0.100795 0.19594152 -0.51441 0.6071
## SmokerFormer:Time_f19 0.120318 0.19093508 0.63015 0.5288
##
## Correlation:
## (Intr) SmkrFr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerFormer -0.461
## Time_f3 -0.733 0.338
## Time_f6 -0.728 0.336 0.533
## Time_f9 -0.726 0.335 0.532 0.529
## Time_f12 -0.725 0.334 0.531 0.527 0.526
## Time_f15 -0.713 0.329 0.522 0.519 0.518 0.517
## Time_f19 -0.724 0.334 0.530 0.527 0.526 0.524 0.516
## SmokerFormer:Time_f3 0.337 -0.730 -0.460 -0.245 -0.245 -0.244 -0.240 -0.244
## SmokerFormer:Time_f6 0.337 -0.730 -0.247 -0.462 -0.245 -0.244 -0.240 -0.244
## SmokerFormer:Time_f9 0.339 -0.734 -0.248 -0.247 -0.466 -0.245 -0.241 -0.245
## SmokerFormer:Time_f12 0.334 -0.725 -0.245 -0.243 -0.243 -0.462 -0.238 -0.242
## SmokerFormer:Time_f15 0.319 -0.691 -0.234 -0.232 -0.232 -0.231 -0.447 -0.231
## SmokerFormer:Time_f19 0.327 -0.709 -0.240 -0.238 -0.238 -0.237 -0.233 -0.452
## SF:T_3 SF:T_6 SF:T_9 SF:T_12 SF:T_15
## SmokerFormer
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerFormer:Time_f3
## SmokerFormer:Time_f6 0.533
## SmokerFormer:Time_f9 0.536 0.536
## SmokerFormer:Time_f12 0.529 0.529 0.532
## SmokerFormer:Time_f15 0.505 0.504 0.507 0.501
## SmokerFormer: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 residual7.5 m4
summary(m4 <- gls(FEV1 ~ Smoker*Time_f,
correlation = corSymm(form = ~ 1 | ID), #General Correlation Structure
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.1548288 1.0981265 1.0095142 1.1970072 1.0796136 1.1241875
## 6*Current 9*Current 12*Current 15*Current 19*Current 3*Current 12*Former
## 1.0518873 1.1013847 1.0381033 1.0158930 0.9676358 1.0363215 1.0527571
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 3.226769 0.06124627 52.68515 0.0000
## SmokerFormer 0.257911 0.11599381 2.22349 0.0265
## Time_f3 -0.089415 0.03460039 -2.58423 0.0099
## Time_f6 -0.157330 0.03459407 -4.54788 0.0000
## Time_f9 -0.381782 0.03582233 -10.65764 0.0000
## Time_f12 -0.409790 0.03620096 -11.31986 0.0000
## Time_f15 -0.582716 0.03532015 -16.49811 0.0000
## Time_f19 -0.725744 0.03664039 -19.80722 0.0000
## SmokerFormer:Time_f3 0.065320 0.07460758 0.87551 0.3816
## SmokerFormer:Time_f6 -0.063890 0.07168071 -0.89131 0.3730
## SmokerFormer:Time_f9 0.052926 0.06933340 0.76335 0.4455
## SmokerFormer:Time_f12 0.014888 0.07199584 0.20680 0.8362
## SmokerFormer:Time_f15 0.103703 0.07642538 1.35692 0.1752
## SmokerFormer:Time_f19 0.127043 0.07531936 1.68672 0.0921
##
## Correlation:
## (Intr) SmkrFr Tim_f3 Tim_f6 Tim_f9 Tm_f12 Tm_f15 Tm_f19
## SmokerFormer -0.528
## Time_f3 -0.434 0.229
## Time_f6 -0.399 0.211 0.664
## Time_f9 -0.329 0.174 0.556 0.661
## Time_f12 -0.429 0.227 0.588 0.663 0.667
## Time_f15 -0.439 0.232 0.592 0.639 0.606 0.744
## Time_f19 -0.509 0.269 0.578 0.649 0.581 0.656 0.739
## SmokerFormer:Time_f3 0.201 -0.204 -0.464 -0.308 -0.258 -0.273 -0.274 -0.268
## SmokerFormer:Time_f6 0.193 -0.249 -0.320 -0.483 -0.319 -0.320 -0.309 -0.313
## SmokerFormer:Time_f9 0.170 -0.332 -0.287 -0.341 -0.517 -0.345 -0.313 -0.300
## SmokerFormer:Time_f12 0.216 -0.326 -0.296 -0.333 -0.336 -0.503 -0.374 -0.330
## SmokerFormer:Time_f15 0.203 -0.154 -0.273 -0.295 -0.280 -0.344 -0.462 -0.342
## SmokerFormer:Time_f19 0.248 -0.322 -0.281 -0.316 -0.283 -0.319 -0.360 -0.486
## SF:T_3 SF:T_6 SF:T_9 SF:T_12 SF:T_15
## SmokerFormer
## Time_f3
## Time_f6
## Time_f9
## Time_f12
## Time_f15
## Time_f19
## SmokerFormer:Time_f3
## SmokerFormer:Time_f6 0.620
## SmokerFormer:Time_f9 0.528 0.661
## SmokerFormer:Time_f12 0.539 0.650 0.692
## SmokerFormer:Time_f15 0.532 0.600 0.579 0.724
## SmokerFormer:Time_f19 0.512 0.632 0.604 0.654 0.701
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.00857129 -0.65224535 0.05368457 0.64660355 3.17725689
##
## Residual standard error: 0.5411918
## Degrees of freedom: 771 total; 757 residual8 Model Comparison
#Likelihood ratio test
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 <.0001m4的AIC值最小，較佳。
#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.7983057 0.8523043 0.8927087
## cor(1,3) 0.8031536 0.8559241 0.8953689
## cor(1,4) 0.7945891 0.8498795 0.8911898
## cor(1,5) 0.7852223 0.8419956 0.8847307
## cor(1,6) 0.8043498 0.8615135 0.9028746
## cor(1,7) 0.5926745 0.7501572 0.8523958
## cor(2,3) 0.8467628 0.8884940 0.9193569
## cor(2,4) 0.7905851 0.8468004 0.8888620
## cor(2,5) 0.8009071 0.8539537 0.8937007
## cor(2,6) 0.8118988 0.8671070 0.9069437
## cor(2,7) 0.6010391 0.7537020 0.8533009
## cor(3,4) 0.8536316 0.8939201 0.9235763
## cor(3,5) 0.8508692 0.8914253 0.9214204
## cor(3,6) 0.8482345 0.8962907 0.9297087
## cor(3,7) 0.7172045 0.8426633 0.9152120
## cor(4,5) 0.8518477 0.8924256 0.9223554
## cor(4,6) 0.8182005 0.8724608 0.9113164
## cor(4,7) 0.6580460 0.7872937 0.8714766
## cor(5,6) 0.9056783 0.9358772 0.9566274
## cor(5,7) 0.6474578 0.7869054 0.8753675
## cor(6,7) 0.7611461 0.8638882 0.9243349
## attr(,"label")
## [1] "Correlation structure:"
##
## Variance function:
## lower est. upper
## 3*Former 0.9803992 1.1565109 1.364258
## 6*Former 0.9444468 1.0999910 1.281152
## 9*Former 0.8711339 1.0114322 1.174326
## 15*Former 1.0181613 1.1984629 1.410693
## 19*Former 0.9211522 1.0802560 1.266840
## 0*Current 1.0089939 1.1411998 1.290728
## 6*Current 0.9582028 1.0677374 1.189793
## 9*Current 0.9995390 1.1179829 1.250462
## 12*Current 0.9428850 1.0541173 1.178472
## 15*Current 0.9236644 1.0311238 1.151085
## 19*Current 0.8761923 0.9817489 1.100022
## 3*Current 0.9370533 1.0521248 1.181327
## 12*Former 0.9073768 1.0546249 1.225768
## attr(,"label")
## [1] "Variance function:"
##
## Residual standard error:
## lower est. upper
## 0.4639700 0.5310277 0.60777749 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")
9.1 m4
dta_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.