LLM_Week7
WEN_SUNG
2016年11月7日
library
library(pacman)
pacman::p_load(tidyverse,nlme,lme4,car,Hmisc,magrittr,lattice,foreign)Question 1
dta1 <- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/weightlifting.txt",header = T,na.strings = ".")
str(dta1)## 'data.frame': 37 obs. of 9 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Program: Factor w/ 2 levels "R","W": 1 1 1 1 1 1 1 1 1 1 ...
## $ Day0 : int 79 83 81 81 80 76 81 77 84 74 ...
## $ Day2 : int NA 83 83 81 81 76 84 78 85 75 ...
## $ Day4 : int 79 85 82 81 82 76 83 79 87 78 ...
## $ Day6 : int 80 85 82 82 82 76 83 79 89 78 ...
## $ Day8 : int 80 86 83 82 82 76 85 81 NA 79 ...
## $ Day10 : int 78 87 83 83 NA 76 85 82 NA 78 ...
## $ Day12 : int 80 87 82 81 86 75 85 81 86 78 ...## wide format to long format
dta1L <- gather(dta1, key = "Day", value = "Liftweight", 3:9) %>%
mutate(Day = as.factor(Day))
str(dta1L)## 'data.frame': 259 obs. of 4 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Program : Factor w/ 2 levels "R","W": 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : Factor w/ 7 levels "Day0","Day10",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Liftweight: int 79 83 81 81 80 76 81 77 84 74 ...# set plot theme to black and white
# this removes the default gray background
ot <- theme_set(theme_bw())
#
ggplot(data = dta1L, aes(reorder(Day, Liftweight, mean,na.rm = T), Liftweight)) +
geom_point(pch = 20, col="gray") +
stat_summary(fun.data = "mean_cl_boot", size = rel(1.1)) +
geom_hline(yintercept = mean(dta1L$Liftweight,na.rm=T), lty = 2) +
coord_flip()+
labs(x = "Day", y = "Liftweight")
## mixed-effects analysis
# covariance patterns
# unequal variances and same correlation
summary(r1_1 <- gls(Liftweight ~ Program,
weights = varIdent(form = ~ 1 | Day),
correlation = corCompSymm(form = ~ 1 | ID),
data = dta1L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: Liftweight ~ Program
## Data: dta1L
## AIC BIC logLik
## 926.9555 961.6361 -453.4777
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | ID
## Parameter estimate(s):
## Rho
## 0.9441032
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Day
## Parameter estimates:
## Day0 Day4 Day6 Day8 Day10 Day12 Day2
## 1.000000 1.211367 1.260039 1.337123 1.354596 1.403072 1.142275
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.64349 0.7445963 104.276 0.0000
## ProgramW 1.02776 0.9939638 1.034 0.3022
##
## Correlation:
## (Intr)
## ProgramW -0.749
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.19837982 0.08789843 0.67773094 1.13146200 2.61098316
##
## Residual standard error: 3.897976
## Degrees of freedom: 239 total; 237 residual# unequal variances and unstructured correlations
summary(r1_2 <- gls(Liftweight ~ Program,
weights = varIdent(form = ~ 1 | Day),
correlation = corSymm(form = ~ 1 | ID),
data = dta1L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: Liftweight ~ Program
## Data: dta1L
## AIC BIC logLik
## 851.2957 955.3375 -395.6478
##
## Correlation Structure: General
## Formula: ~1 | ID
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5 6
## 2 0.944
## 3 0.881 0.932
## 4 0.785 0.843 0.948
## 5 0.786 0.852 0.945 0.944
## 6 0.706 0.783 0.892 0.916 0.946
## 7 0.709 0.783 0.893 0.942 0.949 0.961
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Day
## Parameter estimates:
## Day0 Day4 Day6 Day8 Day10 Day12 Day2
## 1.000000 1.098971 1.034351 1.144642 1.179752 1.183181 1.031949
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 80.35433 0.7145933 112.44764 0.0000
## ProgramW 1.48549 0.9475769 1.56767 0.1183
##
## Correlation:
## (Intr)
## ProgramW -0.754
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.4306621 -0.6187217 0.0480145 0.6952864 2.5914899
##
## Residual standard error: 3.225383
## Degrees of freedom: 239 total; 237 residual# unequal variances and unstructured correlations
summary(r1_3 <- lme(Liftweight ~ Program, random = ~ 1 | ID,
weights = varIdent(form = ~ 1 | Day),
correlation = corCompSymm(form = ~ 1 | ID),
data = dta1L,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta1L
## AIC BIC logLik
## 910.9188 949.0675 -444.4594
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 3.341608 1.730967
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | ID
## Parameter estimate(s):
## Rho
## -0.1666666
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Day
## Parameter estimates:
## Day0 Day4 Day6 Day8 Day10 Day12 Day2
## 1.0000000 0.3777633 0.4409408 0.4987618 0.6520739 0.7361915 0.7992358
## Fixed effects: Liftweight ~ Program
## Value Std.Error DF t-value p-value
## (Intercept) 80.90847 0.8358707 202 96.79544 0.0000
## ProgramW 1.35442 1.1096439 35 1.22059 0.2304
## Correlation:
## (Intr)
## ProgramW -0.753
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.09167743 -0.66831698 -0.05790854 0.64805069 2.98611281
##
## Number of Observations: 239
## Number of Groups: 37# comparing models and choose model 2
anova(r1_1, r1_2, r1_3)## Model df AIC BIC logLik Test L.Ratio p-value
## r1_1 1 10 926.9555 961.6361 -453.4777
## r1_2 2 30 851.2957 955.3375 -395.6478 1 vs 2 115.65979 <.0001
## r1_3 3 11 910.9188 949.0675 -444.4594 2 vs 3 97.62311 <.0001Question 2
dta2 <- read.csv("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/claudication.csv",header = T,sep = ",")
str(dta2)## 'data.frame': 38 obs. of 8 variables:
## $ Treatment: Factor w/ 2 levels "ACT","PBO": 1 1 1 1 1 1 1 1 1 1 ...
## $ PID : Factor w/ 38 levels "P101","P102",..: 1 5 9 12 17 22 23 28 29 32 ...
## $ Gender : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
## $ D0 : int 190 98 155 245 182 140 196 162 195 167 ...
## $ D1 : int 212 137 145 228 205 138 185 176 232 187 ...
## $ D2 : int 213 185 196 280 218 187 185 192 199 228 ...
## $ D3 : int 195 215 189 274 194 195 227 230 185 192 ...
## $ D4 : int 248 225 176 260 193 205 180 215 200 210 ...## wide format to long format
dta2L <- gather(dta2, key = "Month", value = "Distance", 4:8) %>%
mutate(Month = as.factor(Month))%>%mutate(Month = as.numeric(Month)-1)
#
theme_set(theme_bw())
# set dodge amount
pd <- position_dodge(width=.2)
#
ggplot(dta2L, aes(Month, Distance,
color = Treatment,group = Treatment)) +
geom_point(alpha = 0.5) +
stat_summary(fun.y = mean, geom = "line") +
facet_wrap(~ Gender) +
stat_summary(fun.y = mean, geom = "point", position = pd) +
stat_summary(fun.data = mean_se, geom = "errorbar", position = pd, linetype = "solid", width = 0.5) +
labs(x = "Month", y = "Walking Distance") +
theme(legend.position = c(.1, .8))
#mean
rst_m <- dta2 %>% group_by(Gender, Treatment) %>% summarise(D0 = mean(D0),
D1 = mean(D1),
D2 = mean(D2),
D3 = mean(D3),
D4 = mean(D4)
)
rst_m[, c(3:7)] = round(rst_m[,c(3:7)], 2)
rst_m## Source: local data frame [4 x 7]
## Groups: Gender [?]
##
## Gender Treatment D0 D1 D2 D3 D4
## <fctr> <fctr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 F ACT 180.38 193.38 207.88 196.75 208.50
## 2 F PBO 165.56 170.11 175.78 168.11 171.22
## 3 M ACT 163.17 179.50 195.58 204.08 207.67
## 4 M PBO 178.33 165.56 173.33 193.67 194.22#
rst_sd <- dta2 %>% group_by(Gender, Treatment) %>% summarise(D0 = sd(D0),
D1 = sd(D1),
D2 = sd(D2),
D3 = sd(D3),
D4 = sd(D4)
)
rst_sd[, c(3:7)] = round(rst_sd[,c(3:7)], 2)
rst_sd## Source: local data frame [4 x 7]
## Groups: Gender [?]
##
## Gender Treatment D0 D1 D2 D3 D4
## <fctr> <fctr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 F ACT 42.18 33.11 58.21 41.66 57.11
## 2 F PBO 41.53 39.13 47.25 33.43 45.34
## 3 M ACT 42.35 34.51 40.15 28.40 28.04
## 4 M PBO 45.01 45.53 43.68 55.38 47.35#cor
cor(subset(dta2L, Treatment == "ACT" & Gender == "F")[,-(1:3)])## Month Distance
## Month 1.0000000 0.1849544
## Distance 0.1849544 1.0000000cor(subset(dta2L, Treatment == "ACT" & Gender == "M")[,-(1:3)])## Month Distance
## Month 1.0000000 0.4279212
## Distance 0.4279212 1.0000000cor(subset(dta2L, Treatment == "PBO" & Gender == "F")[,-(1:3)])## Month Distance
## Month 1.00000000 0.03350857
## Distance 0.03350857 1.00000000cor(subset(dta2L, Treatment == "PBO" & Gender == "M")[,-(1:3)])## Month Distance
## Month 1.0000000 0.1831019
## Distance 0.1831019 1.0000000#gls
summary(r2 <- gls(Distance ~ Month+Gender*Treatment+Month*Treatment,
weights = varIdent(form = ~ 1 | Month),
correlation = corSymm(form = ~ 1 | PID),
data = dta2L))## Generalized least squares fit by REML
## Model: Distance ~ Month + Gender * Treatment + Month * Treatment
## Data: dta2L
## AIC BIC logLik
## 1781.019 1848.533 -869.5095
##
## Correlation Structure: General
## Formula: ~1 | PID
## Parameter estimate(s):
## Correlation:
## 1 2 3 4
## 2 0.876
## 3 0.774 0.721
## 4 0.777 0.640 0.672
## 5 0.750 0.730 0.704 0.846
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Month
## Parameter estimates:
## 0 1 2 3 4
## 1.0000000 0.8838348 1.0627783 0.9500921 1.0337160
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 177.00300 12.719156 13.916254 0.0000
## Month 10.28326 1.679140 6.124126 0.0000
## GenderM 0.01961 15.640605 0.001254 0.9990
## TreatmentPBO -11.14603 17.575891 -0.634166 0.5268
## GenderM:TreatmentPBO 0.20352 22.484786 0.009051 0.9928
## Month:TreatmentPBO -7.78950 2.439733 -3.192766 0.0017
##
## Correlation:
## (Intr) Month GendrM TrtPBO GM:TPB
## Month -0.304
## GenderM -0.738 0.000
## TreatmentPBO -0.724 0.220 0.534
## GenderM:TreatmentPBO 0.513 0.000 -0.696 -0.702
## Month:TreatmentPBO 0.210 -0.688 0.000 -0.320 0.000
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.82526759 -0.58373590 -0.06888869 0.53349861 3.70599149
##
## Residual standard error: 43.29371
## Degrees of freedom: 190 total; 184 residual# residual plot
plot(r2, pch = 20, cex = 0.5, col = "steelblue",
xlab = "Fitted values", ylab = "Standard residuals")
Question 3
dta3 <- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/adas.txt",header = T,na.strings = ".")
str(dta3)## 'data.frame': 80 obs. of 8 variables:
## $ Treatment: Factor w/ 3 levels "H","L","P": 2 2 2 2 2 2 2 2 2 2 ...
## $ PID : int 1 5 8 12 13 15 19 21 24 28 ...
## $ adas02 : int 22 34 40 24 29 31 22 43 18 25 ...
## $ adas04 : int 30 35 41 NA 26 36 27 49 28 24 ...
## $ adas06 : int NA 46 41 21 29 41 28 42 29 27 ...
## $ adas08 : int 33 37 46 28 26 46 24 48 NA 18 ...
## $ adas10 : int 28 31 52 30 NA 52 27 48 25 21 ...
## $ adas12 : int 30 35 48 27 36 57 28 46 28 22 ...## wide format to long format
dta3L <- gather(dta3, key = "Month", value = "ADAS_Score", 3:8) %>%
mutate(Month = as.factor(Month))
str(dta3L)## 'data.frame': 480 obs. of 4 variables:
## $ Treatment : Factor w/ 3 levels "H","L","P": 2 2 2 2 2 2 2 2 2 2 ...
## $ PID : int 1 5 8 12 13 15 19 21 24 28 ...
## $ Month : Factor w/ 6 levels "adas02","adas04",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ADAS_Score: int 22 34 40 24 29 31 22 43 18 25 ...# set plot theme to black and white
# this removes the default gray background
ot <- theme_set(theme_bw())
#
ggplot(data = dta3L, aes(reorder(Month, ADAS_Score, mean,na.rm = T), ADAS_Score)) +
geom_point(pch = 20, col="gray") +
stat_summary(fun.data = "mean_cl_boot", size = rel(1.1)) +
geom_hline(yintercept = mean(dta3L$ADAS_Score,na.rm=T), lty = 2) +
coord_flip()+
labs(x = "Month", y = "ADAS_Score")
#
summary(r0a <- gls(ADAS_Score ~ PID + Treatment, data = dta3L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: ADAS_Score ~ PID + Treatment
## Data: dta3L
## AIC BIC logLik
## 3290.392 3310.939 -1640.196
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 31.79872 1.0884562 29.214511 0.0000
## PID 0.05823 0.0182114 3.197687 0.0015
## TreatmentL 1.12909 1.0660286 1.059156 0.2901
## TreatmentP 2.78567 1.0321796 2.698826 0.0072
##
## Correlation:
## (Intr) PID TrtmnL
## PID -0.712
## TreatmentL -0.537 0.047
## TreatmentP -0.558 0.053 0.534
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.0480058 -0.7374923 -0.1365497 0.6915341 3.2289073
##
## Residual standard error: 8.983972
## Degrees of freedom: 454 total; 450 residual# unequal variances and same correlation
summary(r3_1 <- gls(ADAS_Score ~ Treatment,
weights = varIdent(form = ~ 1 | Month),
correlation = corCompSymm(form = ~ 1 | PID),
data = dta3L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: ADAS_Score ~ Treatment
## Data: dta3L
## AIC BIC logLik
## 2916.627 2957.742 -1448.314
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | PID
## Parameter estimate(s):
## Rho
## 0.7556982
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Month
## Parameter estimates:
## adas02 adas04 adas08 adas10 adas12 adas06
## 1.0000000 0.8015304 0.8184054 0.8225707 0.8885696 0.7852920
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 35.73185 1.601851 22.306600 0.0000
## TreatmentL 0.62621 2.222792 0.281721 0.7783
## TreatmentP 2.83402 2.149556 1.318422 0.1880
##
## Correlation:
## (Intr) TrtmnL
## TreatmentL -0.721
## TreatmentP -0.745 0.537
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.0483831 -0.8808792 -0.2953486 0.5155460 3.0078538
##
## Residual standard error: 10.95083
## Degrees of freedom: 454 total; 451 residual# unequal variances and unstructured correlations
summary(r3_2 <- gls(ADAS_Score ~ Treatment,
weights = varIdent(form = ~ 1 | Month),
correlation = corSymm(form = ~ 1 | PID),
data = dta3L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: ADAS_Score ~ Treatment
## Data: dta3L
## AIC BIC logLik
## 2812.098 2910.773 -1382.049
##
## Correlation Structure: General
## Formula: ~1 | PID
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5
## 2 0.886
## 3 0.717 0.806
## 4 0.704 0.759 0.869
## 5 0.624 0.689 0.728 0.857
## 6 0.530 0.580 0.691 0.781 0.911
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Month
## Parameter estimates:
## adas02 adas04 adas08 adas10 adas12 adas06
## 1.0000000 0.8422718 0.9016296 0.8819580 0.9387067 0.8328143
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 35.41130 1.617306 21.895238 0.0000
## TreatmentL 0.98772 2.245404 0.439884 0.6602
## TreatmentP 1.83945 2.171285 0.847171 0.3973
##
## Correlation:
## (Intr) TrtmnL
## TreatmentL -0.720
## TreatmentP -0.745 0.537
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.9780505 -0.8344523 -0.2181716 0.5682578 3.0553916
##
## Residual standard error: 10.31642
## Degrees of freedom: 454 total; 451 residual# unequal variances and unstructured correlations
summary(r3_3 <- lme(ADAS_Score ~ Treatment, random = ~ 1 | PID,
weights = varIdent(form = ~ 1 | Month),
correlation = corCompSymm(form = ~ 1 | PID),
data = dta3L,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta3L
## AIC BIC logLik
## 2907.195 2952.421 -1442.597
##
## Random effects:
## Formula: ~1 | PID
## (Intercept) Residual
## StdDev: 8.210979 5.806203
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | PID
## Parameter estimate(s):
## Rho
## -0.199997
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Month
## Parameter estimates:
## adas02 adas04 adas08 adas10 adas12 adas06
## 1.0000000 0.7042305 0.5212886 0.7059348 0.8118391 0.6161154
## Fixed effects: ADAS_Score ~ Treatment
## Value Std.Error DF t-value p-value
## (Intercept) 34.64769 1.678060 374 20.647464 0.0000
## TreatmentL 0.57193 2.327562 77 0.245722 0.8066
## TreatmentP 2.84263 2.252024 77 1.262256 0.2107
## Correlation:
## (Intr) TrtmnL
## TreatmentL -0.721
## TreatmentP -0.745 0.537
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.65800067 -0.62167375 -0.02375653 0.63506177 4.02101735
##
## Number of Observations: 454
## Number of Groups: 80# comparing models and choose model 2
anova(r3_1, r3_2, r3_3)## Model df AIC BIC logLik Test L.Ratio p-value
## r3_1 1 10 2916.628 2957.742 -1448.314
## r3_2 2 24 2812.098 2910.773 -1382.049 1 vs 2 132.5296 <.0001
## r3_3 3 11 2907.195 2952.421 -1442.597 2 vs 3 121.0969 <.0001# residual plots
plot(r0a, resid(., type = "pearson") ~ fitted(.) | Treatment,
layout = c(3, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
ylim = c(-2, 3.3),
xlab = "Ftted values", ylab = "Pearson residuals", main = "Linear model")
plot(r3_2, resid(., type = "pearson") ~ fitted(.) | Treatment,
layout = c(3, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
ylim = c(-2, 3.3),
xlab = "Ftted values", ylab = "Pearson residuals", main = "Covariance patterns")
plot(r3_3, resid(., type = "pearson") ~ fitted(.) | Treatment,
layout = c(3, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
ylim = c(-2, 3.3),
xlab = "Ftted values", ylab = "Pearson residuals",
main = "Covariance patterns + random effects")
Question 4
##Q4
dta4 <- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/grip_strength.txt",header = T)
dta4 <- dta4 %>% spread(key = Time, value = Strength)
str(dta4)## 'data.frame': 67 obs. of 7 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Baseline: int 120 80 60 64 40 50 80 40 70 80 ...
## $ Treat : int 2 1 2 1 1 2 2 1 1 2 ...
## $ Gender : Factor w/ 2 levels "F","M": 2 1 1 1 1 1 1 1 1 1 ...
## $ 1 : int 130 80 80 80 60 70 75 50 90 100 ...
## $ 2 : int 150 86 60 80 NA 70 90 30 110 80 ...
## $ 3 : int 120 80 60 70 NA 70 90 40 90 90 ...dta4L <- dta4%>%gather(key = "Visit",value = "Grip_strength",c(5,6,7))
str(dta4L)## 'data.frame': 201 obs. of 6 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Baseline : int 120 80 60 64 40 50 80 40 70 80 ...
## $ Treat : int 2 1 2 1 1 2 2 1 1 2 ...
## $ Gender : Factor w/ 2 levels "F","M": 2 1 1 1 1 1 1 1 1 1 ...
## $ Visit : chr "1" "1" "1" "1" ...
## $ Grip_strength: int 130 80 80 80 60 70 75 50 90 100 ...dta4L$Treat <- as.factor(dta4L$Treat)
dta4L$Visit <- as.factor(dta4L$Visit)
#xyplot
xyplot(Grip_strength ~ Visit | Treat+Gender, groups = ID, data = dta4L,
type = c("p", "l", "g"),layout = c(2, 2),
xlab = "Time", ylab = "Grip strength (in mmHg)")
#
summary(r4_1 <- lme(Grip_strength ~ Baseline + Treat*Gender*Visit+Baseline:Treat+Baseline:Gender+Baseline:Visit - 1,
random = ~ Visit | ID, data = dta4L,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta4L
## AIC BIC logLik
## 1739.298 1814.838 -845.649
##
## Random effects:
## Formula: ~Visit | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 23.10026 (Intr) Visit2
## Visit2 27.52581 -0.309
## Visit3 34.45645 -0.279 0.868
## Residual 10.25896
##
## Fixed effects: Grip_strength ~ Baseline + Treat * Gender * Visit + Baseline:Treat + Baseline:Gender + Baseline:Visit - 1
## Value Std.Error DF t-value p-value
## Baseline 0.776908 0.115521 60 6.725241 0.0000
## Treat1 29.078977 11.314780 60 2.569999 0.0127
## Treat2 25.464119 12.393932 60 2.054563 0.0443
## GenderM 20.889637 23.360301 60 0.894237 0.3748
## Visit2 10.444154 10.271547 113 1.016804 0.3114
## Visit3 13.500464 12.743263 113 1.059420 0.2917
## Treat2:GenderM 20.640664 18.874695 60 1.093563 0.2785
## Treat2:Visit2 6.840335 10.348386 113 0.661005 0.5100
## Treat2:Visit3 -1.032639 12.824416 113 -0.080521 0.9360
## GenderM:Visit2 11.235879 14.585693 113 0.770336 0.4427
## GenderM:Visit3 -7.292730 17.807437 113 -0.409533 0.6829
## Baseline:Treat2 -0.037886 0.129920 60 -0.291613 0.7716
## Baseline:GenderM 0.027066 0.134761 60 0.200846 0.8415
## Baseline:Visit2 -0.096869 0.083737 113 -1.156819 0.2498
## Baseline:Visit3 -0.056400 0.101600 113 -0.555116 0.5799
## Treat2:GenderM:Visit2 -7.174264 15.718339 113 -0.456426 0.6490
## Treat2:GenderM:Visit3 -3.659046 19.206881 113 -0.190507 0.8493
## Correlation:
## Baseln Treat1 Treat2 GendrM Visit2 Visit3 Tr2:GM
## Treat1 -0.859
## Treat2 -0.420 0.361
## GenderM 0.098 -0.212 -0.518
## Visit2 0.139 -0.256 -0.110 -0.023
## Visit3 0.119 -0.222 -0.095 -0.021 0.772
## Treat2:GenderM 0.278 -0.081 0.280 -0.612 -0.081 -0.072
## Treat2:Visit2 0.004 0.132 -0.127 -0.052 -0.540 -0.417 0.162
## Treat2:Visit3 0.004 0.115 -0.110 -0.046 -0.418 -0.557 0.142
## GenderM:Visit2 0.141 -0.025 -0.102 -0.193 0.073 0.057 0.150
## GenderM:Visit3 0.124 -0.021 -0.090 -0.172 0.058 0.042 0.133
## Baseline:Treat2 -0.392 0.337 -0.551 0.525 0.007 0.007 -0.748
## Baseline:GenderM -0.573 0.492 0.667 -0.807 0.013 0.013 0.245
## Baseline:Visit2 -0.215 0.185 0.164 0.125 -0.703 -0.543 -0.006
## Baseline:Visit3 -0.189 0.163 0.145 0.113 -0.555 -0.689 -0.004
## Treat2:GenderM:Visit2 -0.013 -0.078 0.089 0.106 0.326 0.252 -0.251
## Treat2:GenderM:Visit3 -0.011 -0.069 0.078 0.094 0.256 0.344 -0.223
## Tr2:V2 Tr2:V3 GnM:V2 GnM:V3 Bsl:T2 Bsl:GM Bsl:V2
## Treat1
## Treat2
## GenderM
## Visit2
## Visit3
## Treat2:GenderM
## Treat2:Visit2
## Treat2:Visit3 0.773
## GenderM:Visit2 0.321 0.248
## GenderM:Visit3 0.252 0.341 0.785
## Baseline:Treat2 0.000 -0.001 -0.005 -0.004
## Baseline:GenderM -0.014 -0.012 0.002 0.003 -0.280
## Baseline:Visit2 0.054 0.041 -0.612 -0.481 -0.002 -0.007
## Baseline:Visit3 0.042 0.052 -0.484 -0.606 -0.004 -0.009 0.790
## Treat2:GenderM:Visit2 -0.656 -0.507 -0.585 -0.459 0.003 0.007 0.006
## Treat2:GenderM:Visit3 -0.514 -0.666 -0.459 -0.590 0.003 0.005 0.005
## Bsl:V3 T2:GM:V2
## Treat1
## Treat2
## GenderM
## Visit2
## Visit3
## Treat2:GenderM
## Treat2:Visit2
## Treat2:Visit3
## GenderM:Visit2
## GenderM:Visit3
## Baseline:Treat2
## Baseline:GenderM
## Baseline:Visit2
## Baseline:Visit3
## Treat2:GenderM:Visit2 0.005
## Treat2:GenderM:Visit3 0.006 0.784
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.357879712 -0.208960189 0.001527691 0.219684500 1.888380279
##
## Number of Observations: 189
## Number of Groups: 67#
summary(r4_2 <- lme(Grip_strength ~ Baseline + Treat*Gender*Visit+Baseline:Treat+Baseline:Gender+Baseline:Visit - 1,
random = ~ Visit | ID,
correlation = corAR1(form = ~ Visit | ID),
control = lmeControl(opt="optim"),
data = dta4L,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta4L
## AIC BIC logLik
## 1741.298 1819.985 -845.649
##
## Random effects:
## Formula: ~Visit | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 23.10027 (Intr) Visit2
## Visit2 27.52595 -0.309
## Visit3 34.45625 -0.279 0.868
## Residual 10.25900
##
## Correlation Structure: AR(1)
## Formula: ~Visit | ID
## Parameter estimate(s):
## Phi
## 7.849881e-07
## Fixed effects: Grip_strength ~ Baseline + Treat * Gender * Visit + Baseline:Treat + Baseline:Gender + Baseline:Visit - 1
## Value Std.Error DF t-value p-value
## Baseline 0.776909 0.115521 60 6.725269 0.0000
## Treat1 29.078883 11.314759 60 2.569996 0.0127
## Treat2 25.463980 12.393902 60 2.054557 0.0443
## GenderM 20.889953 23.360224 60 0.894253 0.3748
## Visit2 10.444153 10.271591 113 1.016800 0.3114
## Visit3 13.500423 12.743208 113 1.059421 0.2917
## Treat2:GenderM 20.640587 18.874651 60 1.093561 0.2785
## Treat2:Visit2 6.840333 10.348431 113 0.661002 0.5100
## Treat2:Visit3 -1.032604 12.824360 113 -0.080519 0.9360
## GenderM:Visit2 11.235876 14.585759 113 0.770332 0.4427
## GenderM:Visit3 -7.292701 17.807360 113 -0.409533 0.6829
## Baseline:Treat2 -0.037886 0.129919 60 -0.291609 0.7716
## Baseline:GenderM 0.027064 0.134761 60 0.200831 0.8415
## Baseline:Visit2 -0.096869 0.083738 113 -1.156814 0.2498
## Baseline:Visit3 -0.056400 0.101599 113 -0.555118 0.5799
## Treat2:GenderM:Visit2 -7.174257 15.718407 113 -0.456424 0.6490
## Treat2:GenderM:Visit3 -3.659142 19.206797 113 -0.190513 0.8492
## Correlation:
## Baseln Treat1 Treat2 GendrM Visit2 Visit3 Tr2:GM
## Treat1 -0.859
## Treat2 -0.420 0.361
## GenderM 0.098 -0.212 -0.518
## Visit2 0.139 -0.256 -0.110 -0.023
## Visit3 0.119 -0.222 -0.095 -0.021 0.772
## Treat2:GenderM 0.278 -0.081 0.280 -0.612 -0.081 -0.072
## Treat2:Visit2 0.004 0.132 -0.127 -0.052 -0.540 -0.417 0.162
## Treat2:Visit3 0.004 0.115 -0.110 -0.046 -0.418 -0.557 0.142
## GenderM:Visit2 0.141 -0.025 -0.102 -0.193 0.073 0.057 0.150
## GenderM:Visit3 0.124 -0.021 -0.090 -0.172 0.058 0.042 0.133
## Baseline:Treat2 -0.392 0.337 -0.551 0.525 0.007 0.007 -0.748
## Baseline:GenderM -0.573 0.492 0.667 -0.807 0.013 0.013 0.245
## Baseline:Visit2 -0.215 0.185 0.164 0.125 -0.703 -0.543 -0.006
## Baseline:Visit3 -0.189 0.163 0.145 0.113 -0.555 -0.689 -0.004
## Treat2:GenderM:Visit2 -0.013 -0.078 0.089 0.106 0.326 0.252 -0.251
## Treat2:GenderM:Visit3 -0.011 -0.069 0.078 0.094 0.256 0.344 -0.223
## Tr2:V2 Tr2:V3 GnM:V2 GnM:V3 Bsl:T2 Bsl:GM Bsl:V2
## Treat1
## Treat2
## GenderM
## Visit2
## Visit3
## Treat2:GenderM
## Treat2:Visit2
## Treat2:Visit3 0.773
## GenderM:Visit2 0.321 0.248
## GenderM:Visit3 0.252 0.341 0.785
## Baseline:Treat2 0.000 -0.001 -0.005 -0.004
## Baseline:GenderM -0.014 -0.012 0.002 0.003 -0.280
## Baseline:Visit2 0.054 0.041 -0.612 -0.481 -0.002 -0.007
## Baseline:Visit3 0.042 0.052 -0.484 -0.606 -0.004 -0.009 0.790
## Treat2:GenderM:Visit2 -0.656 -0.507 -0.585 -0.459 0.003 0.007 0.006
## Treat2:GenderM:Visit3 -0.514 -0.666 -0.459 -0.590 0.003 0.005 0.005
## Bsl:V3 T2:GM:V2
## Treat1
## Treat2
## GenderM
## Visit2
## Visit3
## Treat2:GenderM
## Treat2:Visit2
## Treat2:Visit3
## GenderM:Visit2
## GenderM:Visit3
## Baseline:Treat2
## Baseline:GenderM
## Baseline:Visit2
## Baseline:Visit3
## Treat2:GenderM:Visit2 0.005
## Treat2:GenderM:Visit3 0.006 0.784
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.357874314 -0.208985790 0.001524979 0.219679711 1.888346656
##
## Number of Observations: 189
## Number of Groups: 67#
summary(r4_3 <- lme(Grip_strength ~ Baseline + Treat*Gender*Visit+Baseline:Treat+Baseline:Gender+Baseline:Visit - 1,
random = ~ Visit | ID,
weights = varIdent(form = ~ 1 | factor(Visit)),
correlation = corAR1(form = ~ Visit | ID),
control = lmeControl(opt="optim"),
data = dta4L,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta4L
## AIC BIC logLik
## 1745.298 1830.28 -845.649
##
## Random effects:
## Formula: ~Visit | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 23.10028 (Intr) Visit2
## Visit2 27.52596 -0.309
## Visit3 34.45627 -0.279 0.868
## Residual 10.25900
##
## Correlation Structure: AR(1)
## Formula: ~Visit | ID
## Parameter estimate(s):
## Phi
## 7.849881e-07
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | factor(Visit)
## Parameter estimates:
## 1 2 3
## 1.0000000 0.9999964 0.9999983
## Fixed effects: Grip_strength ~ Baseline + Treat * Gender * Visit + Baseline:Treat + Baseline:Gender + Baseline:Visit - 1
## Value Std.Error DF t-value p-value
## Baseline 0.776909 0.115521 60 6.725267 0.0000
## Treat1 29.078881 11.314764 60 2.569994 0.0127
## Treat2 25.463975 12.393908 60 2.054556 0.0443
## GenderM 20.889965 23.360234 60 0.894253 0.3748
## Visit2 10.444153 10.271592 113 1.016800 0.3114
## Visit3 13.500425 12.743212 113 1.059421 0.2917
## Treat2:GenderM 20.640582 18.874660 60 1.093560 0.2785
## Treat2:Visit2 6.840333 10.348432 113 0.661002 0.5100
## Treat2:Visit3 -1.032606 12.824364 113 -0.080519 0.9360
## GenderM:Visit2 11.235876 14.585760 113 0.770332 0.4427
## GenderM:Visit3 -7.292702 17.807366 113 -0.409533 0.6829
## Baseline:Treat2 -0.037886 0.129919 60 -0.291609 0.7716
## Baseline:GenderM 0.027064 0.134761 60 0.200830 0.8415
## Baseline:Visit2 -0.096869 0.083738 113 -1.156814 0.2498
## Baseline:Visit3 -0.056400 0.101599 113 -0.555118 0.5799
## Treat2:GenderM:Visit2 -7.174257 15.718409 113 -0.456424 0.6490
## Treat2:GenderM:Visit3 -3.659136 19.206803 113 -0.190513 0.8492
## Correlation:
## Baseln Treat1 Treat2 GendrM Visit2 Visit3 Tr2:GM
## Treat1 -0.859
## Treat2 -0.420 0.361
## GenderM 0.098 -0.212 -0.518
## Visit2 0.139 -0.256 -0.110 -0.023
## Visit3 0.119 -0.222 -0.095 -0.021 0.772
## Treat2:GenderM 0.278 -0.081 0.280 -0.612 -0.081 -0.072
## Treat2:Visit2 0.004 0.132 -0.127 -0.052 -0.540 -0.417 0.162
## Treat2:Visit3 0.004 0.115 -0.110 -0.046 -0.418 -0.557 0.142
## GenderM:Visit2 0.141 -0.025 -0.102 -0.193 0.073 0.057 0.150
## GenderM:Visit3 0.124 -0.021 -0.090 -0.172 0.058 0.042 0.133
## Baseline:Treat2 -0.392 0.337 -0.551 0.525 0.007 0.007 -0.748
## Baseline:GenderM -0.573 0.492 0.667 -0.807 0.013 0.013 0.245
## Baseline:Visit2 -0.215 0.185 0.164 0.125 -0.703 -0.543 -0.006
## Baseline:Visit3 -0.189 0.163 0.145 0.113 -0.555 -0.689 -0.004
## Treat2:GenderM:Visit2 -0.013 -0.078 0.089 0.106 0.326 0.252 -0.251
## Treat2:GenderM:Visit3 -0.011 -0.069 0.078 0.094 0.256 0.344 -0.223
## Tr2:V2 Tr2:V3 GnM:V2 GnM:V3 Bsl:T2 Bsl:GM Bsl:V2
## Treat1
## Treat2
## GenderM
## Visit2
## Visit3
## Treat2:GenderM
## Treat2:Visit2
## Treat2:Visit3 0.773
## GenderM:Visit2 0.321 0.248
## GenderM:Visit3 0.252 0.341 0.785
## Baseline:Treat2 0.000 -0.001 -0.005 -0.004
## Baseline:GenderM -0.014 -0.012 0.002 0.003 -0.280
## Baseline:Visit2 0.054 0.041 -0.612 -0.481 -0.002 -0.007
## Baseline:Visit3 0.042 0.052 -0.484 -0.606 -0.004 -0.009 0.790
## Treat2:GenderM:Visit2 -0.656 -0.507 -0.585 -0.459 0.003 0.007 0.006
## Treat2:GenderM:Visit3 -0.514 -0.666 -0.459 -0.590 0.003 0.005 0.005
## Bsl:V3 T2:GM:V2
## Treat1
## Treat2
## GenderM
## Visit2
## Visit3
## Treat2:GenderM
## Treat2:Visit2
## Treat2:Visit3
## GenderM:Visit2
## GenderM:Visit3
## Baseline:Treat2
## Baseline:GenderM
## Baseline:Visit2
## Baseline:Visit3
## Treat2:GenderM:Visit2 0.005
## Treat2:GenderM:Visit3 0.006 0.784
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.357874959 -0.208984186 0.001524996 0.219680069 1.888344109
##
## Number of Observations: 189
## Number of Groups: 67#
anova(r4_1, r4_2, r4_3)## Model df AIC BIC logLik Test L.Ratio p-value
## r4_1 1 24 1739.298 1814.838 -845.649
## r4_2 2 25 1741.298 1819.985 -845.649 1 vs 2 2.69174e-08 0.9999
## r4_3 3 27 1745.298 1830.280 -845.649 2 vs 3 2.25964e-09 1.0000Question 5
dta5 <- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/partners.txt",header = T)
str(dta5)## 'data.frame': 464 obs. of 5 variables:
## $ Sbj : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Gender : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
## $ Session: Factor w/ 4 levels "T1","T2","T3",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Time : Factor w/ 2 levels "B","E": 1 1 1 1 2 2 2 2 1 1 ...
## $ PANAS : int 15 16 16 14 35 35 34 41 14 14 ...#
theme_set(theme_bw())
# set dodge amount
pd <- position_dodge(width=.2)
# means with error bars
ggplot(dta5, aes(Session,PANAS ,
group = Gender:Time,col = Gender:Time)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point" , position = pd) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1,position = pd)+
labs(x = "Session", y = "Mean PANAS",
linetype = "Gender:Time", shape = "Gender:Time") +
theme(legend.justification = c(1, 1), legend.position = "top",
legend.background = element_rect(fill = "white", color = "black"))
# means with error bars
ggplot(dta5, aes(Session,PANAS ,
group = Gender,col = Gender)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point",position = pd) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1,position = pd)+
labs(x = "Session", y = "Mean PANAS",
linetype = "Gender:Time", shape = "Gender:Time") +
theme(legend.justification = c(1, 1), legend.position = "top",
legend.background = element_rect(fill = "white", color = "black"))
# means with error bars
ggplot(dta5, aes(as.numeric(Session),PANAS)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point",position = pd) +
stat_summary(fun.data = mean_se, geom = "errorbar",width = 0.1,position = pd) +
labs(x = "Session", y = "Mean PANAS",
linetype = "Gender:Time", shape = "Gender:Time") +
theme(legend.justification = c(1, 1), legend.position = "top",
legend.background = element_rect(fill = "white", color = "black"))
#不會做MANOVA所以改作其他的
summary(r5_1 <- lme(PANAS ~ Gender+Session+Time - 1,
random = ~ 1 | Sbj, data = dta5,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta5
## AIC BIC logLik
## 3253.617 3286.632 -1618.808
##
## Random effects:
## Formula: ~1 | Sbj
## (Intercept) Residual
## StdDev: 3.962854 7.657665
##
## Fixed effects: PANAS ~ Gender + Session + Time - 1
## Value Std.Error DF t-value p-value
## GenderF 22.578762 1.1772571 27 19.179126 0.0000
## GenderM 23.064963 1.8090177 27 12.749993 0.0000
## SessionT2 -0.603448 1.0055007 432 -0.600147 0.5487
## SessionT3 -1.250000 1.0055007 432 -1.243162 0.2145
## SessionT4 -0.827586 1.0055007 432 -0.823059 0.4109
## TimeE 8.228448 0.7109964 432 11.573123 0.0000
## Correlation:
## GendrF GendrM SssnT2 SssnT3 SssnT4
## GenderM 0.237
## SessionT2 -0.427 -0.278
## SessionT3 -0.427 -0.278 0.500
## SessionT4 -0.427 -0.278 0.500 0.500
## TimeE -0.302 -0.197 0.000 0.000 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.4566888 -0.6545219 -0.1100626 0.5937830 3.2426039
##
## Number of Observations: 464
## Number of Groups: 29Question 6
dta6 <- read.dta("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/bodyfat.dta")
str(dta6)## 'data.frame': 2273 obs. of 4 variables:
## $ id : num 1 1 1 2 2 2 3 3 3 4 ...
## $ sex : Factor w/ 2 levels "female","male": 2 2 2 1 1 1 1 1 1 1 ...
## $ age : num 11 13 15 11 13 15 11 13 15 11 ...
## $ bodyfat: num 4 6.2 10.5 8.1 10.4 ...
## - attr(*, "datalabel")= chr ""
## - attr(*, "time.stamp")= chr "21 May 2011 13:30"
## - attr(*, "formats")= chr "%9.0g" "%9.0g" "%9.0g" "%9.0g"
## - attr(*, "types")= int 254 254 254 254
## - attr(*, "val.labels")= chr "" "labsex" "" ""
## - attr(*, "var.labels")= chr "group(u)" "" "" ""
## - attr(*, "version")= int 8
## - attr(*, "label.table")=List of 1
## ..$ labsex: Named int 1 2
## .. ..- attr(*, "names")= chr "female" "male"dta6$id <- as.factor(dta6$id)
dta6$age <- as.factor(dta6$age)
#
xyplot(bodyfat ~ age | sex, groups = id, data = dta6,
type = c("p", "l", "g"),layout = c(2, 1),
xlab = "Time (year)", ylab = "Body fat (kg)")
#
summary(r6_1 <- lme(bodyfat ~ sex*age - 1,
random = ~ 1 | id, data = dta6,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta6
## AIC BIC logLik
## 10264.06 10309.87 -5124.029
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.141043 2.063225
##
## Fixed effects: bodyfat ~ sex * age - 1
## Value Std.Error DF t-value p-value
## sexfemale 8.771233 0.1234091 790 71.07444 0.0000
## sexmale 7.937471 0.1140984 790 69.56687 0.0000
## age13 1.374895 0.1552779 1478 8.85442 0.0000
## age15 2.542678 0.1569622 1478 16.19930 0.0000
## sexmale:age13 -0.463229 0.2114588 1478 -2.19064 0.0286
## sexmale:age15 -0.667826 0.2133657 1478 -3.12996 0.0018
## Correlation:
## sexfml sexmal age13 age15 sxm:13
## sexmale 0.000
## age13 -0.609 0.000
## age15 -0.602 0.000 0.484
## sexmale:age13 0.447 -0.413 -0.734 -0.355
## sexmale:age15 0.443 -0.410 -0.356 -0.736 0.485
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.9441717 -0.6248713 -0.1009279 0.5744466 3.8090258
##
## Number of Observations: 2273
## Number of Groups: 792#
summary(r6_2 <- lme(bodyfat ~ sex*age - 1,
random = ~ 1 | id,
correlation = corAR1(form = ~ 1 | id),
control = lmeControl(opt="optim"),
data = dta6,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta6
## AIC BIC logLik
## 10203.94 10255.48 -5092.971
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.04751357 2.362014
##
## Correlation Structure: AR(1)
## Formula: ~1 | id
## Parameter estimate(s):
## Phi
## 0.3394161
## Fixed effects: bodyfat ~ sex * age - 1
## Value Std.Error DF t-value p-value
## sexfemale 8.775579 0.1236299 790 70.98265 0.0000
## sexmale 7.937479 0.1143064 790 69.44038 0.0000
## age13 1.368614 0.1445918 1478 9.46537 0.0000
## age15 2.538273 0.1682051 1478 15.09034 0.0000
## sexmale:age13 -0.454637 0.1969132 1478 -2.30882 0.0211
## sexmale:age15 -0.660145 0.2287396 1478 -2.88601 0.0040
## Correlation:
## sexfml sexmal age13 age15 sxm:13
## sexmale 0.000
## age13 -0.564 0.000
## age15 -0.648 0.000 0.563
## sexmale:age13 0.414 -0.383 -0.734 -0.413
## sexmale:age15 0.477 -0.441 -0.414 -0.735 0.564
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.34888377 -0.65674924 -0.09117541 0.58628198 4.34306107
##
## Number of Observations: 2273
## Number of Groups: 792#choose model 2
anova(r6_1,r6_2)## Model df AIC BIC logLik Test L.Ratio p-value
## r6_1 1 8 10264.06 10309.87 -5124.029
## r6_2 2 9 10203.94 10255.48 -5092.971 1 vs 2 62.11657 <.0001(a)Body fat increases with age for boys and girls (b)Yes, boys increase less body fat with age (c)Yes
dta7 <- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/respiratory.txt",header = T)
str(dta7)## 'data.frame': 72 obs. of 11 variables:
## $ Patient : int 201 202 203 204 205 206 207 208 209 210 ...
## $ Basefev1: num 2.46 3.5 1.96 3.44 2.8 2.36 1.77 2.64 2.3 2.27 ...
## $ Fev11h : num 2.68 3.95 2.28 4.08 4.09 3.79 3.82 3.67 4.12 2.77 ...
## $ Fev12h : num 2.76 3.65 2.34 3.87 3.9 3.97 3.44 3.47 3.71 2.77 ...
## $ Fev13h : num 2.5 2.93 2.29 3.79 3.54 3.78 3.46 3.19 3.57 2.75 ...
## $ Fev14h : num 2.3 2.53 2.43 3.3 3.35 3.69 3.02 2.19 3.49 2.75 ...
## $ Fev15h : num 2.14 3.04 2.06 3.8 3.15 3.31 2.98 2.85 3.64 2.71 ...
## $ Fev16h : num 2.4 3.37 2.18 3.24 3.23 2.83 3.1 2.68 3.38 2.75 ...
## $ Fev17h : num 2.33 3.14 2.28 2.98 3.46 2.72 2.79 2.6 2.28 2.52 ...
## $ Fev18h : num 2.2 2.62 2.29 2.91 3.27 3 2.88 2.73 3.72 2.6 ...
## $ Drug : Factor w/ 3 levels "a","c","p": 1 1 1 1 1 1 1 1 1 1 ...dta7L <- gather(dta7,key = "Fev", value = "Respiratory",3:10)
str(dta7L)## 'data.frame': 576 obs. of 5 variables:
## $ Patient : int 201 202 203 204 205 206 207 208 209 210 ...
## $ Basefev1 : num 2.46 3.5 1.96 3.44 2.8 2.36 1.77 2.64 2.3 2.27 ...
## $ Drug : Factor w/ 3 levels "a","c","p": 1 1 1 1 1 1 1 1 1 1 ...
## $ Fev : chr "Fev11h" "Fev11h" "Fev11h" "Fev11h" ...
## $ Respiratory: num 2.68 3.95 2.28 4.08 4.09 3.79 3.82 3.67 4.12 2.77 ...dta7L$Fev <- dta7L$Fev%>%as.factor()%>%as.numeric()
dta7L$Fev <- dta7L$Fev+10
dta7L$Fev <- dta7L$Fev%>%as.factor()
#
theme_set(theme_bw())
# set dodge amount
pd <- position_dodge(width=.2)
#means with error bars
ggplot(dta7L, aes(Fev,Respiratory ,
group = Drug,col = Drug)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point",position = pd) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0,position = pd) +
labs(x = "Time (hour)", y = "Respiratory measure (FEV1)") +
theme(legend.justification = c(1, 1), legend.position = c(1,1),
legend.background = element_rect(fill = "white", color = "black"))
#simple
summary(r7_1 <- lm(Respiratory ~ Basefev1 + Drug*Fev , data = dta7L))##
## Call:
## lm(formula = Respiratory ~ Basefev1 + Drug * Fev, data = dta7L)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.31978 -0.36568 -0.04599 0.29073 1.57268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.07964 0.15146 7.128 3.21e-12 ***
## Basefev1 0.90285 0.04076 22.153 < 2e-16 ***
## Drugc 0.21845 0.14910 1.465 0.143471
## Drugp -0.64441 0.14910 -4.322 1.84e-05 ***
## Fev12 -0.07667 0.14910 -0.514 0.607317
## Fev13 -0.28958 0.14910 -1.942 0.052618 .
## Fev14 -0.42708 0.14910 -2.864 0.004337 **
## Fev15 -0.42000 0.14910 -2.817 0.005022 **
## Fev16 -0.49417 0.14910 -3.314 0.000979 ***
## Fev17 -0.60500 0.14910 -4.058 5.67e-05 ***
## Fev18 -0.61542 0.14910 -4.128 4.23e-05 ***
## Drugc:Fev12 0.01208 0.21086 0.057 0.954322
## Drugp:Fev12 0.14250 0.21086 0.676 0.499442
## Drugc:Fev13 0.17583 0.21086 0.834 0.404699
## Drugp:Fev13 0.36083 0.21086 1.711 0.087594 .
## Drugc:Fev14 0.17958 0.21086 0.852 0.394760
## Drugp:Fev14 0.47167 0.21086 2.237 0.025692 *
## Drugc:Fev15 -0.02167 0.21086 -0.103 0.918195
## Drugp:Fev15 0.36167 0.21086 1.715 0.086866 .
## Drugc:Fev16 -0.11167 0.21086 -0.530 0.596613
## Drugp:Fev16 0.48292 0.21086 2.290 0.022383 *
## Drugc:Fev17 -0.10917 0.21086 -0.518 0.604855
## Drugp:Fev17 0.56292 0.21086 2.670 0.007817 **
## Drugc:Fev18 -0.06542 0.21086 -0.310 0.756494
## Drugp:Fev18 0.52083 0.21086 2.470 0.013810 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5165 on 551 degrees of freedom
## Multiple R-squared: 0.5507, Adjusted R-squared: 0.5311
## F-statistic: 28.14 on 24 and 551 DF, p-value: < 2.2e-16# Compound symmetric
summary(r7_2 <- gls(Respiratory ~ Basefev1 + Drug*Fev,
weights = varIdent(form = ~ 1 | Fev),
correlation = corCompSymm(form = ~ 1 | Patient),
data = dta7L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: Respiratory ~ Basefev1 + Drug * Fev
## Data: dta7L
## AIC BIC logLik
## 525.4094 672.0084 -228.7047
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | Patient
## Parameter estimate(s):
## Rho
## 0.6989026
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Fev
## Parameter estimates:
## 11 12 13 14 15 16 17 18
## 1.000000 1.051613 1.033677 1.088561 1.131269 1.041099 1.130513 1.091243
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.2819259 0.23055804 9.897403 0.0000
## Basefev1 0.4522764 0.07614196 5.939910 0.0000
## Drugc 0.2073685 0.08459570 2.451288 0.0145
## Drugp -0.6590510 0.08461119 -7.789171 0.0000
## Fev12 -0.0766667 0.08691233 -0.882115 0.3781
## Fev13 -0.2895833 0.08606559 -3.364682 0.0008
## Fev14 -0.4270833 0.08876691 -4.811290 0.0000
## Fev15 -0.4200000 0.09108549 -4.611053 0.0000
## Fev16 -0.4941667 0.08641165 -5.718750 0.0000
## Fev17 -0.6050000 0.09104289 -6.645220 0.0000
## Fev18 -0.6154167 0.08890710 -6.922019 0.0000
## Drugc:Fev12 0.0120833 0.12273254 0.098453 0.9216
## Drugp:Fev12 0.1425000 0.12273254 1.161061 0.2461
## Drugc:Fev13 0.1758333 0.12163775 1.445549 0.1489
## Drugp:Fev13 0.3608333 0.12163775 2.966458 0.0031
## Drugc:Fev14 0.1795833 0.12501562 1.436487 0.1514
## Drugp:Fev14 0.4716667 0.12501562 3.772862 0.0002
## Drugc:Fev15 -0.0216667 0.12769899 -0.169670 0.8653
## Drugp:Fev15 0.3616667 0.12769899 2.832181 0.0048
## Drugc:Fev16 -0.1116667 0.12208973 -0.914628 0.3608
## Drugp:Fev16 0.4829167 0.12208973 3.955424 0.0001
## Drugc:Fev17 -0.1091667 0.12765109 -0.855196 0.3928
## Drugp:Fev17 0.5629167 0.12765109 4.409807 0.0000
## Drugc:Fev18 -0.0654167 0.12518274 -0.522569 0.6015
## Drugp:Fev18 0.5208333 0.12518274 4.160584 0.0000
##
## Correlation:
## (Intr) Basfv1 Drugc Drugp Fev12 Fev13 Fev14 Fev15 Fev16
## Basefev1 -0.881
## Drugc -0.203 0.022
## Drugp -0.209 0.029 0.500
## Fev12 -0.157 0.000 0.486 0.486
## Fev13 -0.166 0.000 0.491 0.491 0.480
## Fev14 -0.139 0.000 0.476 0.476 0.468 0.471
## Fev15 -0.118 0.000 0.464 0.464 0.459 0.461 0.454
## Fev16 -0.162 0.000 0.489 0.489 0.479 0.482 0.470 0.460
## Fev17 -0.119 0.000 0.464 0.464 0.459 0.461 0.454 0.448 0.460
## Fev18 -0.138 0.000 0.476 0.475 0.468 0.471 0.462 0.454 0.470
## Drugc:Fev12 0.126 0.000 -0.689 -0.344 -0.706 -0.339 -0.328 -0.320 -0.337
## Drugp:Fev12 0.126 0.000 -0.344 -0.689 -0.706 -0.339 -0.328 -0.320 -0.337
## Drugc:Fev13 0.128 0.000 -0.695 -0.348 -0.338 -0.707 -0.331 -0.323 -0.340
## Drugp:Fev13 0.128 0.000 -0.348 -0.695 -0.338 -0.707 -0.331 -0.323 -0.340
## Drugc:Fev14 0.124 0.000 -0.676 -0.338 -0.329 -0.332 -0.704 -0.314 -0.331
## Drugp:Fev14 0.124 0.000 -0.338 -0.676 -0.329 -0.332 -0.704 -0.314 -0.331
## Drugc:Fev15 0.121 0.000 -0.662 -0.331 -0.322 -0.325 -0.316 -0.701 -0.324
## Drugp:Fev15 0.121 0.000 -0.331 -0.662 -0.322 -0.325 -0.316 -0.701 -0.324
## Drugc:Fev16 0.127 0.000 -0.693 -0.346 -0.337 -0.340 -0.330 -0.322 -0.706
## Drugp:Fev16 0.127 0.000 -0.346 -0.692 -0.337 -0.340 -0.330 -0.322 -0.706
## Drugc:Fev17 0.122 0.000 -0.662 -0.331 -0.322 -0.326 -0.316 -0.308 -0.324
## Drugp:Fev17 0.122 0.000 -0.331 -0.662 -0.322 -0.326 -0.316 -0.308 -0.324
## Drugc:Fev18 0.124 0.000 -0.675 -0.338 -0.329 -0.332 -0.322 -0.314 -0.331
## Drugp:Fev18 0.124 0.000 -0.338 -0.675 -0.329 -0.332 -0.322 -0.314 -0.331
## Fev17 Fev18 Drgc:F12 Drgp:F12 Drgc:F13 Drgp:F13 Drgc:F14
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18 0.454
## Drugc:Fev12 -0.320 -0.328
## Drugp:Fev12 -0.320 -0.328 0.500
## Drugc:Fev13 -0.323 -0.331 0.479 0.240
## Drugp:Fev13 -0.323 -0.331 0.240 0.479 0.500
## Drugc:Fev14 -0.314 -0.322 0.466 0.233 0.470 0.235
## Drugp:Fev14 -0.314 -0.322 0.233 0.466 0.235 0.470 0.500
## Drugc:Fev15 -0.308 -0.315 0.456 0.228 0.460 0.230 0.448
## Drugp:Fev15 -0.308 -0.315 0.228 0.456 0.230 0.460 0.224
## Drugc:Fev16 -0.322 -0.329 0.477 0.239 0.482 0.241 0.469
## Drugp:Fev16 -0.322 -0.329 0.239 0.477 0.241 0.482 0.234
## Drugc:Fev17 -0.701 -0.315 0.457 0.228 0.461 0.230 0.448
## Drugp:Fev17 -0.701 -0.315 0.228 0.457 0.230 0.461 0.224
## Drugc:Fev18 -0.314 -0.704 0.466 0.233 0.470 0.235 0.457
## Drugp:Fev18 -0.314 -0.704 0.233 0.466 0.235 0.470 0.229
## Drgp:F14 Drgc:F15 Drgp:F15 Drgc:F16 Drgp:F16 Drgc:F17 Drgp:F17
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18
## Drugc:Fev12
## Drugp:Fev12
## Drugc:Fev13
## Drugp:Fev13
## Drugc:Fev14
## Drugp:Fev14
## Drugc:Fev15 0.224
## Drugp:Fev15 0.448 0.500
## Drugc:Fev16 0.234 0.459 0.229
## Drugp:Fev16 0.469 0.229 0.459 0.500
## Drugc:Fev17 0.224 0.439 0.219 0.459 0.229
## Drugp:Fev17 0.448 0.219 0.439 0.229 0.459 0.500
## Drugc:Fev18 0.229 0.447 0.224 0.468 0.234 0.448 0.224
## Drugp:Fev18 0.457 0.224 0.447 0.234 0.468 0.224 0.448
## Drgc:F18
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18
## Drugc:Fev12
## Drugp:Fev12
## Drugc:Fev13
## Drugp:Fev13
## Drugc:Fev14
## Drugp:Fev14
## Drugc:Fev15
## Drugp:Fev15
## Drugc:Fev16
## Drugp:Fev16
## Drugc:Fev17
## Drugp:Fev17
## Drugc:Fev18
## Drugp:Fev18 0.500
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.7066266 -0.6213056 0.0123327 0.6798032 2.9802010
##
## Residual standard error: 0.5339236
## Degrees of freedom: 576 total; 551 residual# Autoregressive (1)
summary(r7_3 <- gls(Respiratory ~ Basefev1 + Drug*Fev,
weights = varIdent(form = ~ 1 | Fev),
correlation = corAR1(form = ~ 1 | Patient),
data = dta7L,na.action = "na.exclude"))## Generalized least squares fit by REML
## Model: Respiratory ~ Basefev1 + Drug * Fev
## Data: dta7L
## AIC BIC logLik
## 768.2548 914.8538 -350.1274
##
## Correlation Structure: AR(1)
## Formula: ~1 | Patient
## Parameter estimate(s):
## Phi
## 0.6313349
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Fev
## Parameter estimates:
## 11 12 13 14 15 16 17 18
## 1.000000 1.089682 1.147065 1.330363 1.295462 1.152723 1.342478 1.214749
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.7711765 0.18096930 9.787165 0.0000
## Basefev1 0.6436878 0.05820899 11.058219 0.0000
## Drugc 0.2120740 0.07976043 2.658887 0.0081
## Drugp -0.6528301 0.10187416 -6.408202 0.0000
## Fev12 -0.0766667 0.11889811 -0.644810 0.5193
## Fev13 -0.2895833 0.13682494 -2.116451 0.0348
## Fev14 -0.4270833 0.15337920 -2.784493 0.0055
## Fev15 -0.4200000 0.15169308 -2.768748 0.0058
## Fev16 -0.4941667 0.14165532 -3.488515 0.0005
## Fev17 -0.6050000 0.15544862 -3.891961 0.0001
## Fev18 -0.6154167 0.14612160 -4.211675 0.0000
## Drugc:Fev12 0.0120833 0.12218644 0.098893 0.9213
## Drugp:Fev12 0.1425000 0.16427524 0.867447 0.3861
## Drugc:Fev13 0.1758333 0.12246315 1.435806 0.1516
## Drugp:Fev13 0.3608333 0.15862693 2.274729 0.0233
## Drugc:Fev14 0.1795833 0.13301948 1.350053 0.1776
## Drugp:Fev14 0.4716667 0.17049504 2.766454 0.0059
## Drugc:Fev15 -0.0216667 0.13058297 -0.165923 0.8683
## Drugp:Fev15 0.3616667 0.16693480 2.166515 0.0307
## Drugc:Fev16 -0.1116667 0.12171503 -0.917443 0.3593
## Drugp:Fev16 0.4829167 0.15549487 3.105676 0.0020
## Drugc:Fev17 -0.1091667 0.13350150 -0.817719 0.4139
## Drugp:Fev17 0.5629167 0.17052267 3.301125 0.0010
## Drugc:Fev18 -0.0654167 0.12547668 -0.521345 0.6023
## Drugp:Fev18 0.5208333 0.16026581 3.249809 0.0012
##
## Correlation:
## (Intr) Basfv1 Drugc Drugp Fev12 Fev13 Fev14 Fev15 Fev16
## Basefev1 -0.858
## Drugc -0.236 0.018
## Drugp -0.297 0.019 0.639
## Fev12 -0.291 0.000 0.481 0.723
## Fev13 -0.323 0.000 0.325 0.440 0.513
## Fev14 -0.304 0.000 0.269 0.350 0.377 0.530
## Fev15 -0.313 0.000 0.265 0.340 0.355 0.422 0.522
## Fev16 -0.336 0.000 0.282 0.361 0.374 0.422 0.427 0.550
## Fev17 -0.307 0.000 0.257 0.328 0.339 0.378 0.364 0.404 0.544
## Fev18 -0.326 0.000 0.273 0.349 0.360 0.401 0.379 0.397 0.453
## Drugc:Fev12 0.104 0.000 -0.705 -0.522 -0.514 -0.065 -0.104 -0.120 -0.132
## Drugp:Fev12 0.127 0.000 -0.459 -0.748 -0.691 -0.027 -0.113 -0.142 -0.160
## Drugc:Fev13 0.133 0.000 -0.665 -0.444 -0.366 -0.448 -0.060 -0.135 -0.164
## Drugp:Fev13 0.168 0.000 -0.427 -0.677 -0.531 -0.580 -0.004 -0.152 -0.203
## Drugc:Fev14 0.129 0.000 -0.603 -0.390 -0.303 -0.258 -0.434 -0.047 -0.140
## Drugp:Fev14 0.165 0.000 -0.386 -0.607 -0.450 -0.344 -0.556 0.019 -0.159
## Drugc:Fev15 0.134 0.000 -0.612 -0.392 -0.297 -0.214 -0.232 -0.430 -0.069
## Drugp:Fev15 0.171 0.000 -0.391 -0.612 -0.446 -0.289 -0.300 -0.550 -0.012
## Drugc:Fev16 0.144 0.000 -0.655 -0.419 -0.316 -0.217 -0.192 -0.238 -0.430
## Drugp:Fev16 0.184 0.000 -0.419 -0.656 -0.475 -0.293 -0.250 -0.305 -0.549
## Drugc:Fev17 0.132 0.000 -0.597 -0.382 -0.288 -0.195 -0.165 -0.176 -0.234
## Drugp:Fev17 0.168 0.000 -0.381 -0.597 -0.432 -0.264 -0.214 -0.225 -0.299
## Drugc:Fev18 0.140 0.000 -0.635 -0.406 -0.306 -0.207 -0.172 -0.173 -0.195
## Drugp:Fev18 0.179 0.000 -0.406 -0.635 -0.460 -0.280 -0.224 -0.221 -0.249
## Fev17 Fev18 Drgc:F12 Drgp:F12 Drgc:F13 Drgp:F13 Drgc:F14
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18 0.535
## Drugc:Fev12 -0.121 -0.129
## Drugp:Fev12 -0.147 -0.157 0.672
## Drugc:Fev13 -0.153 -0.164 0.429 0.241
## Drugp:Fev13 -0.192 -0.207 0.298 0.408 0.648
## Drugc:Fev14 -0.144 -0.158 0.415 0.260 0.357 0.174
## Drugp:Fev14 -0.178 -0.200 0.302 0.427 0.213 0.300 0.641
## Drugc:Fev15 -0.129 -0.159 0.428 0.277 0.396 0.240 0.323
## Drugp:Fev15 -0.144 -0.198 0.316 0.452 0.258 0.387 0.180
## Drugc:Fev16 -0.064 -0.153 0.461 0.300 0.433 0.275 0.384
## Drugp:Fev16 -0.005 -0.177 0.342 0.489 0.288 0.437 0.242
## Drugc:Fev17 -0.429 -0.057 0.421 0.274 0.397 0.254 0.357
## Drugp:Fev17 -0.548 0.006 0.312 0.447 0.264 0.403 0.229
## Drugc:Fev18 -0.230 -0.429 0.448 0.292 0.422 0.271 0.383
## Drugp:Fev18 -0.294 -0.548 0.332 0.476 0.282 0.430 0.247
## Drgp:F14 Drgc:F15 Drgp:F15 Drgc:F16 Drgp:F16 Drgc:F17 Drgp:F17
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18
## Drugc:Fev12
## Drugp:Fev12
## Drugc:Fev13
## Drugp:Fev13
## Drugc:Fev14
## Drugp:Fev14
## Drugc:Fev15 0.144
## Drugp:Fev15 0.253 0.639
## Drugc:Fev16 0.231 0.357 0.168
## Drugp:Fev16 0.369 0.201 0.288 0.639
## Drugc:Fev17 0.225 0.354 0.210 0.347 0.160
## Drugp:Fev17 0.355 0.219 0.336 0.193 0.277 0.639
## Drugc:Fev18 0.244 0.386 0.243 0.406 0.244 0.334 0.150
## Drugp:Fev18 0.384 0.245 0.382 0.252 0.389 0.184 0.262
## Drgc:F18
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18
## Drugc:Fev12
## Drugp:Fev12
## Drugc:Fev13
## Drugp:Fev13
## Drugc:Fev14
## Drugp:Fev14
## Drugc:Fev15
## Drugp:Fev15
## Drugc:Fev16
## Drugp:Fev16
## Drugc:Fev17
## Drugp:Fev17
## Drugc:Fev18
## Drugp:Fev18 0.639
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.89396645 -0.63998818 -0.01229434 0.61277367 3.09066202
##
## Residual standard error: 0.4549799
## Degrees of freedom: 576 total; 551 residual# Autoregressive (1) with random efect for patient
summary(r7_4 <- lme(Respiratory ~ Basefev1 + Drug*Fev, random = ~1|Patient,
weights = varIdent(form = ~ 1 | Fev),
correlation = corAR1(form = ~ 1 | Patient),
data = dta7L,na.action = "na.exclude"))## Linear mixed-effects model fit by REML
## Data: dta7L
## AIC BIC logLik
## 516.3043 667.215 -223.1521
##
## Random effects:
## Formula: ~1 | Patient
## (Intercept) Residual
## StdDev: 0.4810912 0.3000533
##
## Correlation Structure: AR(1)
## Formula: ~1 | Patient
## Parameter estimate(s):
## Phi
## -0.1054362
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Fev
## Parameter estimates:
## 11 12 13 14 15 16 17
## 1.0000000 0.9254995 0.9229697 1.0365893 1.0566369 1.0101240 1.2839431
## 18
## 1.0815457
## Fixed effects: Respiratory ~ Basefev1 + Drug * Fev
## Value Std.Error DF t-value p-value
## (Intercept) 2.2868401 0.23919525 528 9.560558 0.0000
## Basefev1 0.4504347 0.07844998 528 5.741681 0.0000
## Drugc 0.2073232 0.09109027 528 2.276019 0.0232
## Drugp -0.6591109 0.08617285 528 -7.648707 0.0000
## Fev12 -0.0766667 0.08350253 528 -0.918136 0.3590
## Fev13 -0.2895833 0.08334856 528 -3.474365 0.0006
## Fev14 -0.4270833 0.08821680 528 -4.841292 0.0000
## Fev15 -0.4200000 0.08910458 528 -4.713562 0.0000
## Fev16 -0.4941667 0.08705749 528 -5.676326 0.0000
## Fev17 -0.6050000 0.09967669 528 -6.069624 0.0000
## Fev18 -0.6154167 0.09021874 528 -6.821384 0.0000
## Drugc:Fev12 0.0120833 0.12446696 528 0.097081 0.9227
## Drugp:Fev12 0.1425000 0.11427287 528 1.247015 0.2129
## Drugc:Fev13 0.1758333 0.12393064 528 1.418804 0.1565
## Drugp:Fev13 0.3608333 0.11721930 528 3.078276 0.0022
## Drugc:Fev14 0.1795833 0.13116959 528 1.369093 0.1716
## Drugp:Fev14 0.4716667 0.12406201 528 3.801862 0.0002
## Drugc:Fev15 -0.0216667 0.13248962 528 -0.163535 0.8702
## Drugp:Fev15 0.3616667 0.12531052 528 2.886164 0.0041
## Drugc:Fev16 -0.1116667 0.12944581 528 -0.862652 0.3887
## Drugp:Fev16 0.4829167 0.12243163 528 3.944378 0.0001
## Drugc:Fev17 -0.1091667 0.14820929 528 -0.736571 0.4617
## Drugp:Fev17 0.5629167 0.14017840 528 4.015716 0.0001
## Drugc:Fev18 -0.0654167 0.13414627 528 -0.487652 0.6260
## Drugp:Fev18 0.5208333 0.12687739 528 4.105013 0.0000
## Correlation:
## (Intr) Basfv1 Drugc Drugp Fev12 Fev13 Fev14 Fev15 Fev16
## Basefev1 -0.875
## Drugc -0.209 0.021
## Drugp -0.206 0.030 0.473
## Fev12 -0.188 0.000 0.551 0.465
## Fev13 -0.188 0.000 0.546 0.517 0.539
## Fev14 -0.178 0.000 0.516 0.488 0.510 0.510
## Fev15 -0.176 0.000 0.511 0.483 0.505 0.505 0.477
## Fev16 -0.180 0.000 0.523 0.494 0.517 0.517 0.488 0.483
## Fev17 -0.157 0.000 0.457 0.432 0.451 0.452 0.427 0.422 0.432
## Fev18 -0.174 0.000 0.505 0.477 0.498 0.499 0.471 0.467 0.478
## Drugc:Fev12 0.139 0.000 -0.736 -0.309 -0.745 -0.396 -0.378 -0.374 -0.383
## Drugp:Fev12 0.136 0.000 -0.360 -0.717 -0.684 -0.425 -0.368 -0.365 -0.373
## Drugc:Fev13 0.140 0.000 -0.735 -0.347 -0.404 -0.743 -0.375 -0.376 -0.384
## Drugp:Fev13 0.132 0.000 -0.347 -0.735 -0.342 -0.703 -0.395 -0.355 -0.364
## Drugc:Fev14 0.132 0.000 -0.694 -0.328 -0.383 -0.379 -0.743 -0.350 -0.363
## Drugp:Fev14 0.125 0.000 -0.328 -0.694 -0.323 -0.358 -0.703 -0.374 -0.343
## Drugc:Fev15 0.131 0.000 -0.687 -0.325 -0.379 -0.376 -0.354 -0.743 -0.355
## Drugp:Fev15 0.124 0.000 -0.325 -0.687 -0.320 -0.355 -0.335 -0.703 -0.378
## Drugc:Fev16 0.134 0.000 -0.703 -0.333 -0.388 -0.384 -0.363 -0.359 -0.743
## Drugp:Fev16 0.127 0.000 -0.333 -0.703 -0.327 -0.364 -0.343 -0.340 -0.703
## Drugc:Fev17 0.117 0.000 -0.614 -0.290 -0.339 -0.336 -0.317 -0.314 -0.321
## Drugp:Fev17 0.111 0.000 -0.291 -0.614 -0.286 -0.318 -0.300 -0.297 -0.304
## Drugc:Fev18 0.129 0.000 -0.679 -0.321 -0.374 -0.371 -0.350 -0.347 -0.355
## Drugp:Fev18 0.122 0.000 -0.321 -0.679 -0.316 -0.351 -0.331 -0.328 -0.336
## Fev17 Fev18 Drgc:F12 Drgp:F12 Drgc:F13 Drgp:F13 Drgc:F14
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18 0.416
## Drugc:Fev12 -0.335 -0.370
## Drugp:Fev12 -0.326 -0.360 0.459
## Drugc:Fev13 -0.336 -0.371 0.538 0.291
## Drugp:Fev13 -0.317 -0.351 0.224 0.552 0.473
## Drugc:Fev14 -0.317 -0.350 0.511 0.250 0.507 0.268
## Drugp:Fev14 -0.300 -0.331 0.214 0.498 0.238 0.536 0.473
## Drugc:Fev15 -0.314 -0.347 0.506 0.248 0.505 0.239 0.474
## Drugp:Fev15 -0.297 -0.328 0.212 0.493 0.239 0.505 0.223
## Drugc:Fev16 -0.317 -0.355 0.518 0.254 0.517 0.244 0.488
## Drugp:Fev16 -0.346 -0.336 0.217 0.505 0.245 0.517 0.231
## Drugc:Fev17 -0.743 -0.305 0.452 0.222 0.452 0.214 0.427
## Drugp:Fev17 -0.703 -0.336 0.190 0.441 0.214 0.452 0.202
## Drugc:Fev18 -0.310 -0.743 0.500 0.245 0.499 0.236 0.471
## Drugp:Fev18 -0.293 -0.703 0.210 0.487 0.236 0.499 0.223
## Drgp:F14 Drgc:F15 Drgp:F15 Drgc:F16 Drgp:F16 Drgc:F17 Drgp:F17
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18
## Drugc:Fev12
## Drugp:Fev12
## Drugc:Fev13
## Drugp:Fev13
## Drugc:Fev14
## Drugp:Fev14
## Drugc:Fev15 0.255
## Drugp:Fev15 0.504 0.473
## Drugc:Fev16 0.231 0.480 0.257
## Drugp:Fev16 0.488 0.226 0.511 0.473
## Drugc:Fev17 0.202 0.422 0.200 0.429 0.235
## Drugp:Fev17 0.427 0.200 0.422 0.201 0.462 0.473
## Drugc:Fev18 0.223 0.467 0.221 0.478 0.226 0.414 0.229
## Drugp:Fev18 0.471 0.221 0.467 0.226 0.478 0.194 0.447
## Drgc:F18
## Basefev1
## Drugc
## Drugp
## Fev12
## Fev13
## Fev14
## Fev15
## Fev16
## Fev17
## Fev18
## Drugc:Fev12
## Drugp:Fev12
## Drugc:Fev13
## Drugp:Fev13
## Drugc:Fev14
## Drugp:Fev14
## Drugc:Fev15
## Drugp:Fev15
## Drugc:Fev16
## Drugp:Fev16
## Drugc:Fev17
## Drugp:Fev17
## Drugc:Fev18
## Drugp:Fev18 0.473
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.63057306 -0.57913799 0.04592047 0.62812617 3.03337920
##
## Number of Observations: 576
## Number of Groups: 24# unequal variances and unstructured correlations (unconverge)
#summary(r7_5 <- lme(Respiratory ~ Basefev1 + Drug*Fev,random = ~1|Patient,
# weights = varIdent(form = ~ 1 | Fev),
# correlation = corSymm(form = ~ 1 | Patient),
# control = lmeControl(opt="optim"),
# data = dta7L,na.action = "na.exclude"))