1012_In Class Exercise 2
Pei-Yun,Li
Thu Oct 15 20:54:00 2020
1 load package
library(lme4)## Loading required package: Matrix18 subjects (3hrs sleep) 10 observations on each, p64 of bates
the “new” random effects (vs nlme) see Bates book 2010
sleep deprivation 3hrs/night truckers
2 load Data
data(sleepstudy, package="lme4")
dim(sleepstudy) ## [1] 180 3透過dim函數顯示資料框,列與行的長度。
head(sleepstudy)## Reaction Days Subject
## 1 249.5600 0 308
## 2 258.7047 1 308
## 3 250.8006 2 308
## 4 321.4398 3 308
## 5 356.8519 4 308
## 6 414.6901 5 308with(sleepstudy, table(Subject))## Subject
## 308 309 310 330 331 332 333 334 335 337 349 350 351 352 369 370 371 372
## 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 103 Plot
lattice::xyplot(Reaction ~ Days | Subject, data=sleepstudy,
xlab="Days of sleep deprivation",
ylab="Average reaction time (ms)",
type=c("g","p","r"))
繪製“剝奪睡眠天數”對“反應時間”的XY散佈圖
data來源是“sleepstudy”
“g”為在圖形背景添加網格
“r”是為散佈圖添加迴歸線
“p”是標示出點
pacman::p_load(tidyverse, afex, segmented, foreign, haven, lme4, lmerTest)ggplot(sleepstudy, aes(Days, Reaction, group = Subject, color=Subject))+
geom_point()+
stat_smooth(method = "lm", formula = y ~ poly(x, 2)) + facet_wrap( ~ Subject)+
labs(x = "Days", y = "Reaction") +
theme(legend.position = c(1.3, .2))
這部分試了幾次仍無法畫出和老師一樣的“長條狀”的圖
(或許是我用錯函數了?!!)
library(sqldf)## Loading required package: gsubfn## Loading required package: proto## Loading required package: RSQLite sqldf('select * from sleepstudy
order by Subject')## Reaction Days Subject
## 1 249.5600 0 308
## 2 258.7047 1 308
## 3 250.8006 2 308
## 4 321.4398 3 308
## 5 356.8519 4 308
## 6 414.6901 5 308
## 7 382.2038 6 308
## 8 290.1486 7 308
## 9 430.5853 8 308
## 10 466.3535 9 308
## 11 222.7339 0 309
## 12 205.2658 1 309
## 13 202.9778 2 309
## 14 204.7070 3 309
## 15 207.7161 4 309
## 16 215.9618 5 309
## 17 213.6303 6 309
## 18 217.7272 7 309
## 19 224.2957 8 309
## 20 237.3142 9 309
## 21 199.0539 0 310
## 22 194.3322 1 310
## 23 234.3200 2 310
## 24 232.8416 3 310
## 25 229.3074 4 310
## 26 220.4579 5 310
## 27 235.4208 6 310
## 28 255.7511 7 310
## 29 261.0125 8 310
## 30 247.5153 9 310
## 31 321.5426 0 330
## 32 300.4002 1 330
## 33 283.8565 2 330
## 34 285.1330 3 330
## 35 285.7973 4 330
## 36 297.5855 5 330
## 37 280.2396 6 330
## 38 318.2613 7 330
## 39 305.3495 8 330
## 40 354.0487 9 330
## 41 287.6079 0 331
## 42 285.0000 1 331
## 43 301.8206 2 331
## 44 320.1153 3 331
## 45 316.2773 4 331
## 46 293.3187 5 331
## 47 290.0750 6 331
## 48 334.8177 7 331
## 49 293.7469 8 331
## 50 371.5811 9 331
## 51 234.8606 0 332
## 52 242.8118 1 332
## 53 272.9613 2 332
## 54 309.7688 3 332
## 55 317.4629 4 332
## 56 309.9976 5 332
## 57 454.1619 6 332
## 58 346.8311 7 332
## 59 330.3003 8 332
## 60 253.8644 9 332
## 61 283.8424 0 333
## 62 289.5550 1 333
## 63 276.7693 2 333
## 64 299.8097 3 333
## 65 297.1710 4 333
## 66 338.1665 5 333
## 67 332.0265 6 333
## 68 348.8399 7 333
## 69 333.3600 8 333
## 70 362.0428 9 333
## 71 265.4731 0 334
## 72 276.2012 1 334
## 73 243.3647 2 334
## 74 254.6723 3 334
## 75 279.0244 4 334
## 76 284.1912 5 334
## 77 305.5248 6 334
## 78 331.5229 7 334
## 79 335.7469 8 334
## 80 377.2990 9 334
## 81 241.6083 0 335
## 82 273.9472 1 335
## 83 254.4907 2 335
## 84 270.8021 3 335
## 85 251.4519 4 335
## 86 254.6362 5 335
## 87 245.4523 6 335
## 88 235.3110 7 335
## 89 235.7541 8 335
## 90 237.2466 9 335
## 91 312.3666 0 337
## 92 313.8058 1 337
## 93 291.6112 2 337
## 94 346.1222 3 337
## 95 365.7324 4 337
## 96 391.8385 5 337
## 97 404.2601 6 337
## 98 416.6923 7 337
## 99 455.8643 8 337
## 100 458.9167 9 337
## 101 236.1032 0 349
## 102 230.3167 1 349
## 103 238.9256 2 349
## 104 254.9220 3 349
## 105 250.7103 4 349
## 106 269.7744 5 349
## 107 281.5648 6 349
## 108 308.1020 7 349
## 109 336.2806 8 349
## 110 351.6451 9 349
## 111 256.2968 0 350
## 112 243.4543 1 350
## 113 256.2046 2 350
## 114 255.5271 3 350
## 115 268.9165 4 350
## 116 329.7247 5 350
## 117 379.4445 6 350
## 118 362.9184 7 350
## 119 394.4872 8 350
## 120 389.0527 9 350
## 121 250.5265 0 351
## 122 300.0576 1 351
## 123 269.8939 2 351
## 124 280.5891 3 351
## 125 271.8274 4 351
## 126 304.6336 5 351
## 127 287.7466 6 351
## 128 266.5955 7 351
## 129 321.5418 8 351
## 130 347.5655 9 351
## 131 221.6771 0 352
## 132 298.1939 1 352
## 133 326.8785 2 352
## 134 346.8555 3 352
## 135 348.7402 4 352
## 136 352.8287 5 352
## 137 354.4266 6 352
## 138 360.4326 7 352
## 139 375.6406 8 352
## 140 388.5417 9 352
## 141 271.9235 0 369
## 142 268.4369 1 369
## 143 257.2424 2 369
## 144 277.6566 3 369
## 145 314.8222 4 369
## 146 317.2135 5 369
## 147 298.1353 6 369
## 148 348.1229 7 369
## 149 340.2800 8 369
## 150 366.5131 9 369
## 151 225.2640 0 370
## 152 234.5235 1 370
## 153 238.9008 2 370
## 154 240.4730 3 370
## 155 267.5373 4 370
## 156 344.1937 5 370
## 157 281.1481 6 370
## 158 347.5855 7 370
## 159 365.1630 8 370
## 160 372.2288 9 370
## 161 269.8804 0 371
## 162 272.4428 1 371
## 163 277.8989 2 371
## 164 281.7895 3 371
## 165 279.1705 4 371
## 166 284.5120 5 371
## 167 259.2658 6 371
## 168 304.6306 7 371
## 169 350.7807 8 371
## 170 369.4692 9 371
## 171 269.4117 0 372
## 172 273.4740 1 372
## 173 297.5968 2 372
## 174 310.6316 3 372
## 175 287.1726 4 372
## 176 329.6076 5 372
## 177 334.4818 6 372
## 178 343.2199 7 372
## 179 369.1417 8 372
## 180 364.1236 9 372ggplot(sleepstudy, aes(Days, Reaction))+
geom_point()+
stat_smooth(method = "lm", formula = y ~ x) +
facet_wrap( ~ Subject) +
labs(x = "Days", y = "Reaction") +
theme(legend.position = c(1.3, .2))
qplot(Days, Reaction,
data = sleepstudy, geom= c("point","smooth"),
method= "glm", family= binomial, facets = . ~ Subject) +theme_bw()## Warning: Ignoring unknown parameters: method, family## Warning: Ignoring unknown parameters: family## `geom_smooth()` using formula 'y ~ x'
sleeplmList <- lmList(Reaction ~ Days |Subject, data = sleepstudy)sleeplmList ## Call: lmList(formula = Reaction ~ Days | Subject, data = sleepstudy)
## Coefficients:
## (Intercept) Days
## 308 244.1927 21.764702
## 309 205.0549 2.261785
## 310 203.4842 6.114899
## 330 289.6851 3.008073
## 331 285.7390 5.266019
## 332 264.2516 9.566768
## 333 275.0191 9.142045
## 334 240.1629 12.253141
## 335 263.0347 -2.881034
## 337 290.1041 19.025974
## 349 215.1118 13.493933
## 350 225.8346 19.504017
## 351 261.1470 6.433498
## 352 276.3721 13.566549
## 369 254.9681 11.348109
## 370 210.4491 18.056151
## 371 253.6360 9.188445
## 372 267.0448 11.298073
##
## Degrees of freedom: 180 total; 144 residual
## Residual standard error: 25.59182#mean int and slope match lmer Fixed effects results p.67
mean(coef(sleeplmList)[,1]) ## [1] 251.4051mean(coef(sleeplmList)[,2]) ## [1] 10.46729[]“中括號”用於將資料進行篩選,使用[列,行]方式取出直行欄位資料。
(這個區塊指的是取出上步驟–>sleeplmList的Coefficients的平均數)
[,1]表示是“第一行”的平均數(Intercept)
[,2]表示是“第二行”的平均數(Days)
#these are too big as they should be, #compare with variance random effects p.67 # mle subtracts off wobble in estimated indiv regressions Y on t
4 - SSR/SST
var(coef(sleeplmList)[,1]) ## [1] 838.3423var(coef(sleeplmList)[,2]) ## [1] 43.01034決定係數, r2 = SSR/SST,，用以解釋X軸與Y軸變項之直線關係的強弱。 SSR為y軸數值迴歸平方和, SSE為y軸數值殘差平方和, SST為y軸數值之總平方和, SSR + SSE = SST
coef(sleeplmList)[,1]## [1] 244.1927 205.0549 203.4842 289.6851 285.7390 264.2516 275.0191 240.1629
## [9] 263.0347 290.1041 215.1118 225.8346 261.1470 276.3721 254.9681 210.4491
## [17] 253.6360 267.0448coef(sleeplmList)[,2]## [1] 21.764702 2.261785 6.114899 3.008073 5.266019 9.566768 9.142045
## [8] 12.253141 -2.881034 19.025974 13.493933 19.504017 6.433498 13.566549
## [15] 11.348109 18.056151 9.188445 11.2980735 quantile
quantile(coef(sleeplmList)[,1])## 0% 25% 50% 75% 100%
## 203.4842 229.4167 258.0576 273.0255 290.1041quantile(coef(sleeplmList)[,2])## 0% 25% 50% 75% 100%
## -2.881034 6.194548 10.432421 13.548395 21.764702分位數，將資料分為相等的組
對應到第一、二、三、四分位數的數分別如上表所示
6 stem
stem(coef(sleeplmList)[,2])##
## The decimal point is 1 digit(s) to the right of the |
##
## -0 | 3
## 0 | 23
## 0 | 56699
## 1 | 011234
## 1 | 89
## 2 | 02透過stem函數產生莖葉圖。
本區塊以進行了分位數分組的資料作為莖葉圖之依據
圖中可以得知分組如下:
[-3]
[2、3]
[5、6、6、9、9]
[18、19]
[10、11、11、12、13、14]
[20、22]
cor(coef(sleeplmList)[,1], coef(sleeplmList)[,2]) ## [1] -0.1375534計算相關係數，得出相關係數為-0.1375534
因此兩者之間屬於弱負相關
7 first 2 do REML, second pair match Bates p.67 in doing MLE (REML FALSE)
sleeplmer <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
summary(sleeplmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1743.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9536 -0.4634 0.0231 0.4634 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 612.10 24.741
## Days 35.07 5.922 0.07
## Residual 654.94 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 251.405 6.825 17.000 36.838 < 2e-16 ***
## Days 10.467 1.546 17.000 6.771 3.26e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138運用lmer，討論“剝奪睡眠天數”對“反應時間”的影響
解釋:增加每單位的“Days”，對於“Reaction”有-0.138的減少
sleeplmer2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=FALSE)
summary(sleeplmer2)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## AIC BIC logLik deviance df.resid
## 1763.9 1783.1 -876.0 1751.9 174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9416 -0.4656 0.0289 0.4636 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 565.48 23.780
## Days 32.68 5.717 0.08
## Residual 654.95 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 251.405 6.632 18.001 37.907 < 2e-16 ***
## Days 10.467 1.502 18.000 6.968 1.65e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138