Q3
#Is there any distinction in exercise tolerance profiles between patients who receive Novafylline and those on placebo?
dat = read.csv("http://140.116.183.121/~sheu/lmm/Data/claudication.csv", header = TRUE)
str(dat)
## '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 ...
require(magrittr)
## Loading required package: magrittr
require(tidyr)
## Loading required package: tidyr
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
dat = dat %>% gather(key = month, value = dis, D0,D1,D2,D3,D4)
dat = dat %>%separate(month, c("n", "month"), sep = 1)
dat<-dat[-4]
#Plot
require(ggplot2)
## Loading required package: ggplot2
ggplot(dat, aes(x = month, y = dis, group = Gender, color = Treatment))+
stat_summary(aes(color = Treatment, group = Treatment),fun.y = mean, geom = "point")+
stat_summary(aes(color = Treatment, group = Treatment),fun.y = mean, geom = "line")+
stat_summary(aes(color = Treatment, group = Treatment),fun.data = mean_se, geom = "errorbar", width = .2)+
facet_wrap(~Gender)+
labs(x = "months", y = "Distance")+
theme_bw()
#有交互作用
require(nlme)
## Loading required package: nlme
summary(m0 <- gls(dis~month+Gender*Treatment+month:Treatment,
correlation = corSymm(form = ~1|PID),
weights = varIdent(form = ~1|month),
data = dat))
## Generalized least squares fit by REML
## Model: dis ~ month + Gender * Treatment + month:Treatment
## Data: dat
## AIC BIC logLik
## 1749.506 1835.414 -847.7529
##
## Correlation Structure: General
## Formula: ~1 | PID
## Parameter estimate(s):
## Correlation:
## 1 2 3 4
## 2 0.878
## 3 0.781 0.718
## 4 0.770 0.635 0.673
## 5 0.746 0.736 0.715 0.837
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | month
## Parameter estimates:
## 0 1 2 3 4
## 1.0000000 0.8915599 1.0795057 0.9486381 1.0287427
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 170.03789 13.440748 12.650925 0.0000
## month1 15.00000 4.607348 3.255669 0.0014
## month2 30.45000 6.661777 4.570852 0.0000
## month3 31.10000 6.377944 4.876180 0.0000
## month4 37.95000 6.955426 5.456172 0.0000
## GenderM 0.02018 15.640507 0.001290 0.9990
## TreatmentPBO 1.79485 18.675206 0.096109 0.9235
## GenderM:TreatmentPBO 0.20322 22.484646 0.009038 0.9928
## month1:TreatmentPBO -19.11111 6.694322 -2.854824 0.0048
## month2:TreatmentPBO -27.83889 9.679338 -2.876115 0.0045
## month3:TreatmentPBO -22.15556 9.266937 -2.390817 0.0179
## month4:TreatmentPBO -27.17222 10.106000 -2.688722 0.0079
##
## Correlation:
## (Intr) month1 month2 month3 month4 GendrM TrtPBO
## month1 -0.325
## month2 -0.162 0.195
## month3 -0.291 0.076 0.252
## month4 -0.230 0.360 0.366 0.665
## GenderM -0.698 0.000 0.000 0.000 0.000
## TreatmentPBO -0.720 0.234 0.117 0.210 0.165 0.503
## GenderM:TreatmentPBO 0.486 0.000 0.000 0.000 0.000 -0.696 -0.660
## month1:TreatmentPBO 0.223 -0.688 -0.134 -0.052 -0.248 0.000 -0.340
## month2:TreatmentPBO 0.112 -0.134 -0.688 -0.174 -0.252 0.000 -0.170
## month3:TreatmentPBO 0.200 -0.052 -0.174 -0.688 -0.457 0.000 -0.305
## month4:TreatmentPBO 0.158 -0.248 -0.252 -0.457 -0.688 0.000 -0.240
## GM:TPB m1:TPB m2:TPB m3:TPB
## month1
## month2
## month3
## month4
## GenderM
## TreatmentPBO
## GenderM:TreatmentPBO
## month1:TreatmentPBO 0.000
## month2:TreatmentPBO 0.000 0.195
## month3:TreatmentPBO 0.000 0.076 0.252
## month4:TreatmentPBO 0.000 0.360 0.366 0.665
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.8562836 -0.6190432 -0.1082266 0.5224223 3.5519928
##
## Residual standard error: 43.03221
## Degrees of freedom: 190 total; 178 residual
#fit
dat_f = data.frame(dat, fit = fitted(m0))
ggplot(dat_f, aes(x = month, y = fit, group = Gender, color = Treatment))+
stat_summary(aes(color = Treatment, group = Treatment), fun.y = mean, geom = "point")+
stat_summary(aes(color = Treatment, group = Treatment), fun.y = mean, geom = "line")+
facet_wrap(~Gender)+
theme_bw()+
labs(x = "month", y = "distance")+
stat_summary(aes(y = dis, color = Treatment, group = Treatment), fun.y = mean, geom = "point", alpha = .2)+
stat_summary(aes(y = dis, color = Treatment, group = Treatment), fun.y = mean, geom = "line", alpha = .2)+
stat_summary(aes(y = dis, color = Treatment, group = Treatment), fun.data = mean_se, geom = "errorbar", width = .2, alpha = .2)
Q4
dat = read.table("http://140.116.183.121/~sheu/lmm/Data/adas.txt", header = TRUE)
# "." represents missing values
dat$PID %<>% as.factor
#寬變長
require("tidyverse")
## Loading required package: tidyverse
## -- Attaching packages ----------------------------------------------- tidyverse 1.2.1 --
## √ tibble 1.3.4 √ dplyr 0.7.4
## √ readr 1.1.1 √ stringr 1.2.0
## √ purrr 0.2.4 √ forcats 0.2.0
## -- Conflicts -------------------------------------------------- tidyverse_conflicts() --
## x dplyr::collapse() masks nlme::collapse()
## x tidyr::extract() masks magrittr::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x purrr::set_names() masks magrittr::set_names()
dat_l <- dat %>% gather(key = month , value = adas, contains("adas")) %>%
separate(month, c("n", "month"), sep = 4)
## Warning: attributes are not identical across measure variables;
## they will be dropped
dat_l$month %<>% as.factor
dat_l<-dat_l[-3]
#
dat_l$PID <- as.factor(dat_l$PID)
ggplot(dat_l, aes( month, adas, group = PID)) +
geom_point(size= .1 , alpha = .3) +
geom_line(size= .1 , alpha = .2) +
stat_smooth(aes(group = Treatment), method = "lm") +
facet_wrap(~ Treatment) +
labs(x = "month", y = "adas score")+
theme_bw()+
theme(legend.position = "none")
#
dat_l = na.omit(dat_l)
dat_l$adas %<>% as.factor
dat_l$PID %<>% as.factor
dat_l$month %<>% as.integer
summary(m <- gls(adas~Treatment*month, data = dat_l, weights = varIdent(form = ~1|month)))
## Warning in Ops.factor(y, Fitted): '-' not meaningful for factors
## Generalized least squares fit by REML
## Model: adas ~ Treatment * month
## Data: dat_l
## AIC BIC logLik
## 3541.924 3591.858 -1758.962
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | month
## Parameter estimates:
## 1 2 3 4 5 6
## 1.000000 1.048331 1.041631 1.111783 1.034157 1.257259
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 16.102240 1.7571027 9.164086 0.0000
## TreatmentL -2.288507 2.4366630 -0.939197 0.3481
## TreatmentP -1.686697 2.3574007 -0.715490 0.4747
## month 0.400285 0.4724856 0.847190 0.3973
## TreatmentL:month 0.864952 0.6552196 1.320096 0.1874
## TreatmentP:month 1.119439 0.6339059 1.765939 0.0780
##
## Correlation:
## (Intr) TrtmnL TrtmnP month TrtmL:
## TreatmentL -0.721
## TreatmentP -0.745 0.537
## month -0.893 0.644 0.665
## TreatmentL:month 0.644 -0.893 -0.480 -0.721
## TreatmentP:month 0.665 -0.480 -0.893 -0.745 0.537
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## NA NA NA NA NA
##
## Residual standard error: 8.850246
## Degrees of freedom: 480 total; 474 residual
Q5
##The interest is to compare treatment effects and
#to examine if the treatment difference is independent of sex.
dat = read.table("http://140.116.183.121/~sheu/lmm/Data/grip_strength.txt", header = TRUE)
str(dat)
## 'data.frame': 201 obs. of 6 variables:
## $ ID : int 26 26 26 27 27 27 29 29 29 34 ...
## $ Baseline: int 175 175 175 165 165 165 175 175 175 178 ...
## $ Treat : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Gender : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
## $ Time : int 1 2 3 1 2 3 1 2 3 1 ...
## $ Strength: int 161 210 230 215 245 265 134 215 139 165 ...
#題目說有遺漏值
require(ggplot2)
dat$ID = as.factor(dat$ID)
dat$Treat = as.factor(dat$Treat)
ggplot(dat, aes(x = Time, y = Strength, group = ID, color = ID))+
geom_point(size =.5 , alpha = .5 ,na.rm = T)+
geom_line(size = .1, alpha = .5 ,na.rm = T)+
facet_wrap(c("Gender","Treat"))+
labs(x = "Time", y = "Grip Strength")+
theme_bw()+theme(legend.position = "none")
ggplot(dat, aes(x = Time, y = Strength, group = ID, color = Gender))+
geom_point(size =.5 , alpha = .5 ,na.rm = T)+
geom_line(size = .1, alpha = .5 ,na.rm = T)+
facet_wrap("Treat")+
labs(x = "Time", y = "Grip Strength")+
theme_bw()+theme()
#model
require(nlme)
#Strengthijk = ai + bitimek + β1baselinei + β2trtj + β3sexi + β4,5timek
# +β6trtj×sexi + β7,8trtj×timek + β9,10sexi×timek + β11,12trtj×sexj×timek
# +β13baselinei×trtj + β14baselinei×sexi + β15baselinei×timek,
summary(m <- lme(Strength~Baseline+Treat+Gender+Time+
Treat:Gender+Treat:Time+Gender:Time+Treat:Gender:Time+
Baseline:Treat+Baseline:Gender+Baseline:Time,
random = ~Time|Gender,
data = dat, na.action = na.omit))
## Warning in pt(-abs(tVal), fDF): NaNs produced
## Linear mixed-effects model fit by REML
## Data: dat
## AIC BIC logLik
## 1826.897 1877.715 -897.4483
##
## Random effects:
## Formula: ~Time | Gender
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 8.492844 (Intr)
## Time 4.016573 0
## Residual 30.959955
##
## Fixed effects: Strength ~ Baseline + Treat + Gender + Time + Treat:Gender + Treat:Time + Gender:Time + Treat:Gender:Time + Baseline:Treat + Baseline:Gender + Baseline:Time
## Value Std.Error DF t-value p-value
## (Intercept) 21.95848 18.661264 177 1.176688 0.2409
## Baseline 0.82468 0.145669 177 5.661289 0.0000
## Treat2 -5.32378 17.697469 177 -0.300822 0.7639
## GenderM 42.28818 29.771448 0 1.420427 NaN
## Time 5.06855 8.462133 177 0.598968 0.5500
## Treat2:GenderM 21.38202 26.167087 177 0.817134 0.4150
## Treat2:Time 0.95930 7.501764 177 0.127876 0.8984
## GenderM:Time -2.37254 11.878526 177 -0.199733 0.8419
## Baseline:Treat2 -0.00225 0.100730 177 -0.022358 0.9822
## Baseline:GenderM -0.07359 0.105431 177 -0.698010 0.4861
## Baseline:Time -0.02634 0.059517 177 -0.442547 0.6586
## Treat2:GenderM:Time -5.51655 11.204731 177 -0.492341 0.6231
## Correlation:
## (Intr) Baseln Treat2 GendrM Time Tr2:GM Trt2:T
## Baseline -0.658
## Treat2 -0.471 0.141
## GenderM -0.158 -0.296 0.084
## Time -0.663 0.463 0.384 -0.043
## Treat2:GenderM 0.185 0.111 -0.336 -0.535 -0.249
## Treat2:Time 0.414 -0.020 -0.814 -0.231 -0.489 0.555
## GenderM:Time -0.033 0.431 -0.251 -0.610 -0.130 0.433 0.309
## Baseline:Treat2 0.155 -0.241 -0.462 0.310 0.017 -0.421 -0.007
## Baseline:GenderM 0.225 -0.348 0.112 -0.498 0.030 0.129 -0.014
## Baseline:Time 0.522 -0.786 -0.022 0.439 -0.602 -0.011 0.040
## Treat2:GenderM:Time -0.254 -0.021 0.543 0.403 0.302 -0.840 -0.668
## GndM:T Bsl:T2 Bsl:GM Bsln:T
## Baseline
## Treat2
## GenderM
## Time
## Treat2:GenderM
## Treat2:Time
## GenderM:Time
## Baseline:Treat2 0.001
## Baseline:GenderM 0.013 -0.260
## Baseline:Time -0.539 -0.016 -0.033
## Treat2:GenderM:Time -0.519 0.007 0.005 0.016
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.954982711 -0.510757205 -0.002518863 0.533586765 3.715857760
##
## Number of Observations: 189
## Number of Groups: 2
anova(m)
## Warning in pf(Fval[i], nDF[i], dDF[i]): NaNs produced
## numDF denDF F-value p-value
## (Intercept) 1 177 390.2302 <.0001
## Baseline 1 177 279.8499 <.0001
## Treat 1 177 0.0934 0.7602
## Gender 1 0 7.2261 NaN
## Time 1 177 0.0181 0.8932
## Treat:Gender 1 177 1.3570 0.2456
## Treat:Time 1 177 0.0735 0.7866
## Gender:Time 1 177 1.0536 0.3061
## Baseline:Treat 1 177 0.0481 0.8266
## Baseline:Gender 1 177 0.5048 0.4783
## Baseline:Time 1 177 0.1890 0.6642
## Treat:Gender:Time 1 177 0.2424 0.6231
Q7
require(foreign)
## Loading required package: foreign
dat = read.dta("http://140.116.183.121/~sheu/lmm/Data/bodyfat.dta")
str(dat)
## '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"
summary(dat)
## id sex age bodyfat
## Min. : 1.0 female:1045 Min. :11.00 Min. : 2.900
## 1st Qu.:199.0 male :1228 1st Qu.:11.00 1st Qu.: 7.500
## Median :396.0 Median :13.00 Median : 9.100
## Mean :396.5 Mean :12.95 Mean : 9.394
## 3rd Qu.:595.0 3rd Qu.:15.00 3rd Qu.:11.100
## Max. :792.0 Max. :15.00 Max. :19.600
require(ggplot2)
ggplot(dat, aes(x = age, y = bodyfat, group = id, color = id))+
geom_point(size =.1 , alpha = .5)+
geom_line(size = .1, alpha = .5)+
facet_wrap(~sex)+
labs(x = "Time(year)", y = "Body fat (kg)")+
theme_bw()+theme(legend.position = "none")
ggplot(dat, aes(x = age, y = bodyfat, group = sex, color = sex))+
stat_summary(fun.y = mean, geom = "point")+
stat_summary(fun.y = mean, geom = "line")+
stat_summary(fun.data = mean_se, geom = "errorbar", width = .2)+
labs(x = "Time(year)", y = "Body fat (kg)")+
theme_bw()
#男女體重成長的截距與斜率有差異
#將資料進行分組，然後求出各組平均數與標準差
require(magrittr)
require(dplyr)
m <- dat %>% group_by(sex, age) %>% summarise(mean(bodyfat)) #男女平均有差
m
## # A tibble: 6 x 3
## # Groups: sex [?]
## sex age `mean(bodyfat)`
## <fctr> <dbl> <dbl>
## 1 female 11 8.771233
## 2 female 13 10.150289
## 3 female 15 11.323054
## 4 male 11 7.937471
## 5 male 13 8.844691
## 6 male 15 9.815152
s <- dat %>% group_by(sex, age) %>% summarise(sd(bodyfat)) #標準差還好
s
## # A tibble: 6 x 3
## # Groups: sex [?]
## sex age `sd(bodyfat)`
## <fctr> <dbl> <dbl>
## 1 female 11 1.925571
## 2 female 13 2.265062
## 3 female 15 2.797142
## 4 male 11 1.983441
## 5 male 13 2.217271
## 6 male 15 2.853605
#
dat$age_f = as.factor(dat$age)
dat$id = as.factor(dat$id)
## 各種嘗試
require(nlme)
summary(m0 <- gls(bodyfat ~ sex*age, data = dat))
## Generalized least squares fit by REML
## Model: bodyfat ~ sex * age
## Data: dat
## AIC BIC logLik
## 10370.27 10398.9 -5180.134
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.7780638 0.5818214 3.056030 0.0023
## sexmale 0.9879022 0.7914043 1.248290 0.2121
## age 0.6387135 0.0446062 14.318936 0.0000
## sexmale:age -0.1694899 0.0606550 -2.794329 0.0052
##
## Correlation:
## (Intr) sexmal age
## sexmale -0.735
## age -0.992 0.729
## sexmale:age 0.730 -0.992 -0.735
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.37650737 -0.66432158 -0.08668265 0.59228998 4.35814062
##
## Residual standard error: 2.3571
## Degrees of freedom: 2273 total; 2269 residual
summary(m1 <- gls(bodyfat ~ sex*age,
weights = varIdent(form = ~ 1 | age*sex),
data = dat))
## Generalized least squares fit by REML
## Model: bodyfat ~ sex * age
## Data: dat
## AIC BIC logLik
## 10272.54 10329.82 -5126.272
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | age * sex
## Parameter estimates:
## 11*male 13*male 15*male 11*female 13*female 15*female
## 1.0000000 1.1172027 1.4384400 0.9708505 1.1415509 1.4101107
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.6952671 0.5589721 3.032829 0.0025
## sexmale 1.0963401 0.7689107 1.425835 0.1541
## age 0.6453216 0.0444957 14.502996 0.0000
## sexmale:age -0.1781926 0.0611332 -2.914827 0.0036
##
## Correlation:
## (Intr) sexmal age
## sexmale -0.727
## age -0.992 0.721
## sexmale:age 0.722 -0.992 -0.728
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.8517689 -0.7062621 -0.1006567 0.5964325 5.1784674
##
## Residual standard error: 1.983207
## Degrees of freedom: 2273 total; 2269 residual
summary(m2 <- gls(bodyfat ~ sex*age,
correlation = corSymm(form = ~ 1 | id),
weights = varIdent(form = ~ 1 | age*sex),
data = dat))
## Generalized least squares fit by REML
## Model: bodyfat ~ sex * age
## Data: dat
## AIC BIC logLik
## 10082.83 10157.28 -5028.414
##
## Correlation Structure: General
## Formula: ~1 | id
## Parameter estimate(s):
## Correlation:
## 1 2
## 2 0.166
## 3 0.054 0.461
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | age * sex
## Parameter estimates:
## 11*male 13*male 15*male 11*female 13*female 15*female
## 1.0000000 1.1217878 1.4461103 0.9782322 1.1453847 1.4092153
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.8045021 0.5480472 3.292604 0.0010
## sexmale 0.9760194 0.7538982 1.294630 0.1956
## age 0.6358776 0.0444559 14.303576 0.0000
## sexmale:age -0.1677516 0.0611374 -2.743847 0.0061
##
## Correlation:
## (Intr) sexmal age
## sexmale -0.727
## age -0.989 0.719
## sexmale:age 0.719 -0.989 -0.727
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.8516936 -0.7063816 -0.1030066 0.5950403 5.1962283
##
## Residual standard error: 1.976451
## Degrees of freedom: 2273 total; 2269 residual
summary(m3 <- gls(bodyfat ~ sex*age,
correlation = corSymm(form = ~1| id),
weights = varExp(form = ~ age),
data = dat))
## Generalized least squares fit by REML
## Model: bodyfat ~ sex * age
## Data: dat
## AIC BIC logLik
## 10078.37 10129.92 -5030.187
##
## Correlation Structure: General
## Formula: ~1 | id
## Parameter estimate(s):
## Correlation:
## 1 2
## 2 0.169
## 3 0.057 0.467
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~age
## Parameter estimates:
## expon
## 0.09320165
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.8063673 0.5466440 3.304467 0.0010
## sexmale 0.9706584 0.7433525 1.305785 0.1918
## age 0.6355736 0.0444097 14.311598 0.0000
## sexmale:age -0.1670398 0.0603772 -2.766604 0.0057
##
## Correlation:
## (Intr) sexmal age
## sexmale -0.735
## age -0.989 0.727
## sexmale:age 0.727 -0.989 -0.736
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.8424922 -0.7177972 -0.1027411 0.6076323 5.3372886
##
## Residual standard error: 0.6901846
## Degrees of freedom: 2273 total; 2269 residual
summary(m4 <- gls(bodyfat ~ sex*age,
correlation = corSymm(form = ~1| id),
weights = varExp(form = ~ age | sex),
data = dat))
## Generalized least squares fit by REML
## Model: bodyfat ~ sex * age
## Data: dat
## AIC BIC logLik
## 10080.25 10137.52 -5030.124
##
## Correlation Structure: General
## Formula: ~1 | id
## Parameter estimate(s):
## Correlation:
## 1 2
## 2 0.169
## 3 0.057 0.468
## Variance function:
## Structure: Exponential of variance covariate, different strata
## Formula: ~age | sex
## Parameter estimates:
## male female
## 0.09356957 0.09276003
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.8069690 0.5434204 3.325177 0.0009
## sexmale 0.9701560 0.7426651 1.306317 0.1916
## age 0.6355265 0.0441394 14.398179 0.0000
## sexmale:age -0.1670016 0.0603208 -2.768556 0.0057
##
## Correlation:
## (Intr) sexmal age
## sexmale -0.732
## age -0.989 0.724
## sexmale:age 0.724 -0.989 -0.732
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.8612851 -0.7138193 -0.1032827 0.6051643 5.3156166
##
## Residual standard error: 0.6901994
## Degrees of freedom: 2273 total; 2269 residual
anova(m0, m1, m2, m3, m4)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 5 10370.27 10398.90 -5180.134
## m1 2 10 10272.55 10329.82 -5126.272 1 vs 2 107.72342 <.0001
## m2 3 13 10082.83 10157.28 -5028.414 2 vs 3 195.71664 <.0001
## m3 4 9 10078.37 10129.92 -5030.187 3 vs 4 3.54575 0.471
## m4 5 10 10080.25 10137.52 -5030.124 4 vs 5 0.12564 0.723
#m3 is the best
#擬合後的圖
dat_f = data.frame(dat, fit = fitted(m3))
ggplot(dat_f, aes(x = age, y = fit, group = sex, color = sex))+
stat_summary(fun.y = mean, geom = "point")+
stat_summary(fun.y = mean, geom = "line")+
stat_summary(aes(y = bodyfat), fun.y = mean, geom = "point", alpha = .5)+
stat_summary(aes(y = bodyfat), fun.y = mean, geom = "line", alpha = .5, linetype = 2)+
stat_summary(aes(y = bodyfat), fun.data = mean_se, geom = "errorbar", width = .3, alpha = .5, linetype = 2)+
labs(x = "Age", y = "Body fat (kg)")+
theme_bw()
