Q4
#random effect
dat = read.table("C:/Users/USER/Desktop/Multilevel Statistical/w4 exercise/skin_resistance.txt", header = TRUE)
str(dat)
## 'data.frame': 80 obs. of 3 variables:
## $ Electrode: Factor w/ 5 levels "E1","E2","E3",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ Subject : Factor w/ 16 levels "S1","S10","S11",..: 1 1 1 1 1 9 9 9 9 9 ...
## $ Kohm : int 500 400 98 200 250 660 600 600 75 310 ...
#random effect model
require(lme4)
## Loading required package: lme4
## Loading required package: Matrix
m0 = lmer(Kohm ~ Electrode + (1|Subject), data = dat)
summary(m0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Kohm ~ Electrode + (1 | Subject)
## Data: dat
##
## REML criterion at convergence: 999.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4000 -0.4998 -0.0682 0.3076 3.6676
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 14180 119.1
## Residual 22379 149.6
## Number of obs: 80, groups: Subject, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 181.69 47.80 3.801
## ElectrodeE2 105.63 52.89 1.997
## ElectrodeE3 76.31 52.89 1.443
## ElectrodeE4 -31.31 52.89 -0.592
## ElectrodeE5 -43.81 52.89 -0.828
##
## Correlation of Fixed Effects:
## (Intr) ElctE2 ElctE3 ElctE4
## ElectrodeE2 -0.553
## ElectrodeE3 -0.553 0.500
## ElectrodeE4 -0.553 0.500 0.500
## ElectrodeE5 -0.553 0.500 0.500 0.500
#intervals .95
require(nlme)
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
##
## lmList
confint(m0,level = 0.95)
## Computing profile confidence intervals ...
## Warning in if (parm == "theta_") {: 條件的長度 > 1，因此只能用其第一元素
## Warning in if (parm == "beta_") {: 條件的長度 > 1，因此只能用其第一元素
## 2.5 % 97.5 %
## .sig01 70.548336 184.6781
## .sigma 122.962425 174.0792
## (Intercept) 89.126933 274.2481
## ElectrodeE2 3.729397 207.5206
## ElectrodeE3 -25.583103 178.2081
## ElectrodeE4 -133.208103 70.5831
## ElectrodeE5 -145.708103 58.0831
#The interval of Electrode2 did not include zero point
#E2 may be having an undue influence on the results
Q5
#two model to compare
require(mets)
## Loading required package: mets
## Loading required package: timereg
## Loading required package: survival
## Loading required package: lava
## lava version 1.5.1
## mets version 1.2.2
data(dermalridgesMZ)
head(dermalridgesMZ)
## sex left right group id
## 1 female 95 83 psychotic&neurotic 1
## 2 female 90 85 psychotic&neurotic 1
## 3 female 93 90 psychotic&neurotic 2
## 4 female 80 80 psychotic&neurotic 2
## 5 male 99 94 psychotic&neurotic 3
## 6 male 94 99 psychotic&neurotic 3
str(dermalridgesMZ)
## 'data.frame': 36 obs. of 5 variables:
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 2 2 1 1 1 1 ...
## $ left : int 95 90 93 80 99 94 55 48 41 26 ...
## $ right: num 83 85 90 80 94 99 61 54 24 42 ...
## $ group: Factor w/ 2 levels "normal","psychotic&neurotic": 2 2 2 2 2 2 2 2 2 2 ...
## $ id : int 1 1 2 2 3 3 4 4 5 5 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:12] 37 38 39 40 41 42 43 44 45 46 ...
## .. ..- attr(*, "names")= chr [1:12] "37" "38" "39" "40" ...
as.factor(dermalridgesMZ$id)
## [1] 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12
## [24] 12 13 13 14 14 15 15 16 16 17 17 18 18
## Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
m0 = lmer(right~ left+(1|id), data = dermalridgesMZ)
summary(m0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: right ~ left + (1 | id)
## Data: dermalridgesMZ
##
## REML criterion at convergence: 250.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.66773 -0.46281 0.02278 0.59754 1.60380
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.00 0.000
## Residual 63.38 7.961
## Number of obs: 36, groups: id, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.07986 4.13771 2.436
## left 0.85750 0.06022 14.240
##
## Correlation of Fixed Effects:
## (Intr)
## left -0.947
m1 = lm ( right ~ left, data = dermalridgesMZ)
summary(m1)
##
## Call:
## lm(formula = right ~ left, data = dermalridgesMZ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.2374 -3.6843 0.1813 4.7570 12.7676
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.07986 4.13771 2.436 0.0202 *
## left 0.85750 0.06022 14.240 6.86e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.961 on 34 degrees of freedom
## Multiple R-squared: 0.8564, Adjusted R-squared: 0.8522
## F-statistic: 202.8 on 1 and 34 DF, p-value: 6.857e-16
#compare two model
anova(m0,m1)
## refitting model(s) with ML (instead of REML)
## Data: dermalridgesMZ
## Models:
## m1: right ~ left
## m0: right ~ left + (1 | id)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 3 255.47 260.22 -124.74 249.47
## m0 4 257.47 263.81 -124.74 249.47 0 1 1
#p>.05
Q6
#
# variance components
#
# set path to working directory
setwd("C:/Users/USER/Desktop/Multilevel Statistical/w4 exercise")
# package management
require(pacman)
## Loading required package: pacman
#
pacman::p_load(ggplot2, nlme, lme4, dplyr, tidyr)
# input data
dta <- read.table("dog_scans.txt", header = T)
# first 6 lines
head(dta)
## Dog Observer Scan Thickness
## 1 Midge C dog_67 2.60
## 2 Tara C dog_77 4.93
## 3 Tara B dog_77 3.47
## 4 Tara A dog_77 5.29
## 5 Tara C dog_78 2.54
## 6 Tara B dog_78 2.80
# data structure
str(dta)
## 'data.frame': 55 obs. of 4 variables:
## $ Dog : Factor w/ 11 levels "Aussie","Corrie",..: 9 11 11 11 11 11 11 7 7 7 ...
## $ Observer : Factor w/ 3 levels "A","B","C": 3 3 2 1 3 2 1 3 2 1 ...
## $ Scan : Factor w/ 29 levels "dog_100","dog_101",..: 14 15 15 15 16 16 16 17 17 17 ...
## $ Thickness: num 2.6 4.93 3.47 5.29 2.54 2.8 3 3.57 3.3 2.66 ...
# set theme to black and white and save parameters
ot <- theme_set(theme_bw())
# plot
ggplot(dta, aes(reorder(Dog, Thickness, mean), Thickness, color = Scan)) +
geom_point() +
geom_hline(yintercept =mean(dta$Thickness), linetype = "dashed") +
#虛線為mean of Thickness
facet_wrap( ~ Observer) +
labs(y = "Thickness", x = "Dog ID") +
theme(legend.position = "NONE")
# mixed-effect approach
summary(m2 <- lmer(Thickness ~ (1 | Observer) + (1 | Dog) +
(1 | Dog:Observer) + (1 | Scan), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Thickness ~ (1 | Observer) + (1 | Dog) + (1 | Dog:Observer) +
## (1 | Scan)
## Data: dta
##
## REML criterion at convergence: 177.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8144 -0.5141 -0.0410 0.3252 1.6897
##
## Random effects:
## Groups Name Variance Std.Dev.
## Scan (Intercept) 0.70477 0.8395
## Dog:Observer (Intercept) 0.30706 0.5541
## Dog (Intercept) 0.29388 0.5421
## Observer (Intercept) 0.08284 0.2878
## Residual 0.56554 0.7520
## Number of obs: 55, groups:
## Scan, 29; Dog:Observer, 28; Dog, 11; Observer, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.4157 0.3298 13.39
#take into account the interaction between Dog and Observer
# confidence intervals for model parameters
confint(m2)
## Computing profile confidence intervals ...
## Warning in if (parm == "theta_") {: 條件的長度 > 1，因此只能用其第一元素
## Warning in if (parm == "beta_") {: 條件的長度 > 1，因此只能用其第一元素
## 2.5 % 97.5 %
## .sig01 0.2863717 1.357591
## .sig02 0.0000000 1.071654
## .sig03 0.0000000 1.281156
## .sig04 0.0000000 1.065162
## .sigma 0.4839820 1.180072
## (Intercept) 3.7390502 5.170843
# random effects
ranef(m2)
## $Scan
## (Intercept)
## dog_100 0.80990693
## dog_101 -0.49295053
## dog_102 0.31435413
## dog_105 -0.03842648
## dog_106 -0.54508083
## dog_107 -0.74506368
## dog_108 0.67618155
## dog_110 0.92897796
## dog_111 -0.18539685
## dog_112 0.27473124
## dog_113 -0.81398077
## dog_114 0.19778951
## dog_115 0.09132913
## dog_67 -0.42626504
## dog_77 0.39241108
## dog_78 -1.01457482
## dog_79 -1.01069271
## dog_80 1.08532180
## dog_83 0.32642052
## dog_84 -0.14169816
## dog_85 -0.12955776
## dog_86 -0.36689329
## dog_90 0.39267038
## dog_92 -0.86257705
## dog_95 -0.16183142
## dog_96 -0.36696207
## dog_97 1.06929816
## dog_98 0.60834053
## dog_99 0.13421851
##
## $`Dog:Observer`
## (Intercept)
## Aussie:A 0.50616152
## Aussie:C -0.60257703
## Corrie:A -0.18028303
## Corrie:B 0.36867094
## Corrie:C -0.41877924
## Dance:A 0.29433508
## Dance:B 0.07070268
## Dance:C -0.29718507
## Gem:A 0.45073628
## Gem:B 0.04060181
## Gem:C -0.07998978
## Gus:C -0.07781310
## Isla:A -0.01727277
## Isla:B 0.38779768
## Isla:C 0.36040957
## Jenny:A -0.48395797
## Jenny:B 0.34611485
## Jenny:C 0.17035846
## Jos:A -0.34649793
## Jos:B 0.04627890
## Jos:C 0.38070116
## Midge:C -0.41439890
## Mist:A -0.01804194
## Mist:B -0.25223106
## Mist:C 0.03723091
## Tara:A 0.11156981
## Tara:B -0.44833391
## Tara:C 0.06569207
##
## $Dog
## (Intercept)
## Aussie -0.09227512
## Corrie -0.22049759
## Dance 0.06493888
## Gem 0.39368372
## Gus -0.07447156
## Isla 0.69954585
## Jenny 0.03111902
## Jos 0.07702598
## Midge -0.39660331
## Mist -0.22303453
## Tara -0.25943134
##
## $Observer
## (Intercept)
## A 0.08545233
## B 0.15096899
## C -0.23642131
# covariates summarized over groups
###