Week_6_Homework_Answers
Hunter B. Sturgill
2023-05-22
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: Matrix
##
##
## Attaching package: 'Matrix'
##
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
##
## Attaching package: 'psych'
##
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
##
##
##
## Attaching package: 'MASS'
##
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
##
## arm (Version 1.13-1, built: 2022-8-25)
##
##
## Working directory is /Users/user/Desktop/R - Spring
##
##
##
## Attaching package: 'arm'
##
##
## The following objects are masked from 'package:psych':
##
## logit, rescale, simHaptic Example
# Read in the data and create factor variables
hs_data <- read.csv("MLM_haptic_example_data.csv", header = TRUE)
hs_data <- as_tibble(hs_data) %>% mutate(Sub_Num = factor(Subject_Num))
hs_data <- hs_data %>% mutate(Ratio = factor(Relative_Size))
hs_data <- hs_data %>% mutate(Num_Dis = factor(Distractor_Number))why might we want to treat some of our variables as factors?
Answer: It makes them into categories instead of integers. In this case they are identifiers for different groups so it makes sense to make them into categories.
describe(hs_data)## vars n mean sd median trimmed mad min max range
## Experiment_Number 1 3370 0.49 0.50 0.00 0.49 0.00 0.00 1.00 1.00
## Subject_Num 2 3370 39.69 23.19 38.00 39.12 28.17 1.00 85.00 84.00
## Blk 3 3370 1.51 0.50 2.00 1.52 0.00 1.00 2.00 1.00
## Distractor_Number 4 3370 2.00 0.82 2.00 2.00 1.48 1.00 3.00 2.00
## Distractor_Lengths 5 3370 1.00 0.54 0.75 1.00 0.74 0.25 1.75 1.50
## Relative_Size 6 3370 0.98 0.94 0.75 0.82 0.37 0.25 3.00 2.75
## Times 7 3370 14.70 10.26 11.65 13.17 7.89 0.55 60.05 59.50
## Sub_Num* 8 3370 37.52 21.91 36.00 37.05 26.69 1.00 79.00 78.00
## Ratio* 9 3370 3.52 1.70 4.00 3.53 1.48 1.00 6.00 5.00
## Num_Dis* 10 3370 2.00 0.82 2.00 2.00 1.48 1.00 3.00 2.00
## skew kurtosis se
## Experiment_Number 0.03 -2.00 0.01
## Subject_Num 0.17 -1.07 0.40
## Blk -0.06 -2.00 0.01
## Distractor_Number -0.01 -1.51 0.01
## Distractor_Lengths 0.01 -1.51 0.01
## Relative_Size 1.50 0.64 0.02
## Times 1.39 1.83 0.18
## Sub_Num* 0.15 -1.08 0.38
## Ratio* -0.02 -1.27 0.03
## Num_Dis* -0.01 -1.51 0.01#summary(hs_data)Multi Level Visualization
I want to walk through how we visualize multilevel models as a segway from the previous work we have been doing in the course. I am going to fit previous models that we have discussed in class to recap and extend them to transition into the multilevel framework.
Start with the group level model
hs_mdl0 <- lm("Times ~ Distractor_Number", data = hs_data)
display(hs_mdl0)## lm(formula = "Times ~ Distractor_Number", data = hs_data)
## coef.est coef.se
## (Intercept) 12.52 0.47
## Distractor_Number 1.09 0.22
## ---
## n = 3370, k = 2
## residual sd = 10.23, R-Squared = 0.01Now plot that
ggplot(data = hs_data, mapping = aes(x = Distractor_Number, y = Times)) + geom_jitter() + geom_smooth(method = "lm")## `geom_smooth()` using formula = 'y ~ x'
ggplot(data = hs_data, mapping = aes(x = Distractor_Number, y = Times)) + geom_point() + geom_smooth(method = "lm")## `geom_smooth()` using formula = 'y ~ x'
#varying intercepts
hs_mdl0.5 <- lm("Times ~ Num_Dis-1", data = hs_data)
display(hs_mdl0.5)## lm(formula = "Times ~ Num_Dis-1", data = hs_data)
## coef.est coef.se
## Num_Dis1 13.62 0.31
## Num_Dis2 14.69 0.31
## Num_Dis3 15.80 0.30
## ---
## n = 3370, k = 3
## residual sd = 10.23, R-Squared = 0.67ggplot(data = hs_data, mapping = aes(x = 1:3370, y = Times)) + geom_jitter() + geom_smooth(method = "lm", aes(color = Num_Dis))## `geom_smooth()` using formula = 'y ~ x'
Now look at the lower level model
ggplot(data = hs_data, mapping = aes(x = Relative_Size, y = Times)) + geom_jitter(mapping = aes(color = Num_Dis)) + stat_smooth(method = "lm", se=FALSE, color = 'red')## `geom_smooth()` using formula = 'y ~ x'
How do we write this model using the lm() function?
hs_mdl1 <- lm("Times ~ Relative_Size", data = hs_data)
display(hs_mdl1)## lm(formula = "Times ~ Relative_Size", data = hs_data)
## coef.est coef.se
## (Intercept) 17.57 0.25
## Relative_Size -2.92 0.18
## ---
## n = 3370, k = 2
## residual sd = 9.89, R-Squared = 0.07hs_mdl1.5 <- lm("Times ~ Ratio-1", data = hs_data)
display(hs_mdl1.5)## lm(formula = "Times ~ Ratio-1", data = hs_data)
## coef.est coef.se
## Ratio0.25 17.45 0.42
## Ratio0.333333333333333 19.65 0.41
## Ratio0.5 14.88 0.41
## Ratio0.75 14.43 0.41
## Ratio1 12.43 0.41
## Ratio3 9.61 0.41
## ---
## n = 3370, k = 6
## residual sd = 9.74, R-Squared = 0.71ggplot(data = hs_data, mapping = aes(x = 1:3370, y = Times)) + geom_jitter() + stat_smooth(method = "lm", aes(color = Ratio))## `geom_smooth()` using formula = 'y ~ x'
Now how can we represent both
hs_mdl2 <- lm('Times ~ Ratio - 1 + Num_Dis - 1',data= hs_data)
display(hs_mdl2)## lm(formula = "Times ~ Ratio - 1 + Num_Dis - 1", data = hs_data)
## coef.est coef.se
## Ratio0.25 16.35 0.48
## Ratio0.333333333333333 18.53 0.47
## Ratio0.5 13.76 0.47
## Ratio0.75 13.31 0.47
## Ratio1 11.33 0.47
## Ratio3 8.48 0.47
## Num_Dis2 1.13 0.41
## Num_Dis3 2.20 0.41
## ---
## n = 3370, k = 8
## residual sd = 9.71, R-Squared = 0.71all ratios is for num_dis1 and if want 2 then add 1.13 to each, for 3 add 2.20
What is being ignored by all of the above models? Answer: group level contribution
ggplot(data = hs_data, aes(x = Relative_Size, y = Times, color = Num_Dis)) + geom_jitter() + geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
Part 2 Chapter 12.4
Now lets do multilevel modeling using lmer
Varying Intercept Model
hs_lm <- lmer(Times ~ 1 + (1 | Experiment_Number),data=hs_data)
display(hs_lm)## lmer(formula = Times ~ 1 + (1 | Experiment_Number), data = hs_data)
## coef.est coef.se
## 14.70 0.56
##
## Error terms:
## Groups Name Std.Dev.
## Experiment_Number (Intercept) 0.76
## Residual 10.25
## ---
## number of obs: 3370, groups: Experiment_Number, 2
## AIC = 25258.6, DIC = 25253.3
## deviance = 25253.0What is being included in this model?
Answer: This model simply includes a constant term and allows that term to vary by the grouping variable of Experiment.
Varying Intercept Model with 1 predictor
hs_lm_RI <- lmer(Times ~ 1 + Relative_Size +(1 | Experiment_Number),data=hs_data)
display(hs_lm_RI)## lmer(formula = Times ~ 1 + Relative_Size + (1 | Experiment_Number),
## data = hs_data)
## coef.est coef.se
## (Intercept) 17.58 0.68
## Relative_Size -2.94 0.18
##
## Error terms:
## Groups Name Std.Dev.
## Experiment_Number (Intercept) 0.90
## Residual 9.87
## ---
## number of obs: 3370, groups: Experiment_Number, 2
## AIC = 25005.8, DIC = 24995.8
## deviance = 24996.8This model starts with the no-pooling model $ y~ x$ and then adds the varying intercept from the previous model.
What does adding the varying intercept to the no-pooling model allow for?
Answer: This allows the coefficient for the predictor to vary by group.
# From the book
hs_c <- coef(hs_lm_RI)
hs_c## $Experiment_Number
## (Intercept) Relative_Size
## 0 18.19459 -2.942456
## 1 16.96378 -2.942456
##
## attr(,"class")
## [1] "coef.mer"Write the estimated regression line for each of the groups.
Answer: for Experiment 1 = y = 18.19 - 2.94X Experiment 2 = y = 16.96 - 2.94X
Alternatively we can look at the “fixed effect”
hs_cF <- fixef(hs_lm_RI)
hs_cF## (Intercept) Relative_Size
## 17.579188 -2.942456Write the regression line for the fixed effect and explain in words what this line means.
Answer: This is the estimated regression line averaging over the two experiments y = -2.94X + 17.5
Now lets look at the experiment level errors
hs_cE <-ranef(hs_lm_RI)
hs_cE ## $Experiment_Number
## (Intercept)
## 0 0.6154041
## 1 -0.6154041
##
## with conditional variances for "Experiment_Number"What do these two numbers tell us?
Answer: These two numbers tell us how much the intercept is shifted above or below the intercept of the previous regression line (Averaged over experiments line) depending on if you were in experiment 1 or experiment 2
Now lets look at the uncertainty in the estimates of our coefficients
se.fixef(hs_lm_RI)## (Intercept) Relative_Size
## 0.6825168 0.1802499se.ranef(hs_lm_RI)## $Experiment_Number
## (Intercept)
## 0 0.2306844
## 1 0.2339019What might be the cause of the similarity or difference between the standard errors of our coefficients.
Answer: because the number of observations within each group level (experiment) are similar the standard errors should be similar. If they were very different then it would be due to the weighting or ‘shrinking’ due to having different n’s.
summary(hs_lm_RI)## Linear mixed model fit by REML ['lmerMod']
## Formula: Times ~ 1 + Relative_Size + (1 | Experiment_Number)
## Data: hs_data
##
## REML criterion at convergence: 24997.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6295 -0.7019 -0.2522 0.4293 4.4198
##
## Random effects:
## Groups Name Variance Std.Dev.
## Experiment_Number (Intercept) 0.8115 0.9008
## Residual 97.3846 9.8684
## Number of obs: 3370, groups: Experiment_Number, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 17.5792 0.6825 25.76
## Relative_Size -2.9425 0.1802 -16.32
##
## Correlation of Fixed Effects:
## (Intr)
## Relative_Sz -0.259Interpret the fit of the two models and explain which one is better and how you know? Answer: The aic and bic are both reduced in the model which includes a predictor. However, more important is that the chi squared test has an appreciable p value.
Below is an effective way in which to test the difference between the two.
anova(hs_lm, hs_lm_RI)## refitting model(s) with ML (instead of REML)## Data: hs_data
## Models:
## hs_lm: Times ~ 1 + (1 | Experiment_Number)
## hs_lm_RI: Times ~ 1 + Relative_Size + (1 | Experiment_Number)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## hs_lm 3 25259 25277 -12626 25253
## hs_lm_RI 4 25005 25029 -12498 24997 256.16 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(hs_mdl2)
hs_lm_RI2 <- lmer(Times ~ 1 + Relative_Size +(1|Distractor_Number),data=hs_data)
summary(hs_lm_RI2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Times ~ 1 + Relative_Size + (1 | Distractor_Number)
## Data: hs_data
##
## REML criterion at convergence: 24987.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7583 -0.7043 -0.2350 0.4268 4.2826
##
## Random effects:
## Groups Name Variance Std.Dev.
## Distractor_Number (Intercept) 1.147 1.071
## Residual 97.022 9.850
## Number of obs: 3370, groups: Distractor_Number, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 17.5685 0.6651 26.42
## Relative_Size -2.9253 0.1798 -16.27
##
## Correlation of Fixed Effects:
## (Intr)
## Relative_Sz -0.265display(hs_lm_RI2)## lmer(formula = Times ~ 1 + Relative_Size + (1 | Distractor_Number),
## data = hs_data)
## coef.est coef.se
## (Intercept) 17.57 0.67
## Relative_Size -2.93 0.18
##
## Error terms:
## Groups Name Std.Dev.
## Distractor_Number (Intercept) 1.07
## Residual 9.85
## ---
## number of obs: 3370, groups: Distractor_Number, 3
## AIC = 24995.8, DIC = 24986.1
## deviance = 24986.9anova(hs_lm_RI, hs_lm_RI2)## refitting model(s) with ML (instead of REML)## Data: hs_data
## Models:
## hs_lm_RI: Times ~ 1 + Relative_Size + (1 | Experiment_Number)
## hs_lm_RI2: Times ~ 1 + Relative_Size + (1 | Distractor_Number)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## hs_lm_RI 4 25005 25029 -12498 24997
## hs_lm_RI2 4 24995 25019 -12494 24987 9.8545 0Homework
Repeate Part 2 above using your own data. If you do not have data feel free to use the data from our textbook.
Read in your data
my_data <- read.csv("wpa_lena_pos.csv")
describe(my_data)## vars n mean sd median trimmed mad min
## id_uni 1 1488 10383.11 3181.90 11103.00 11228.93 1332.86 801.00
## age 2 1488 11.76 0.94 11.83 11.64 1.27 10.74
## sit_time 3 1488 0.45 0.28 0.41 0.43 0.29 0.00
## held_time 4 1488 0.07 0.13 0.01 0.04 0.02 0.00
## prone_time 5 1488 0.15 0.18 0.10 0.11 0.14 0.00
## supine_time 6 1488 0.14 0.24 0.02 0.08 0.02 0.00
## upright_time 7 1488 0.20 0.22 0.13 0.16 0.19 0.00
## adult_word_cnt 8 1488 231.77 214.07 174.95 200.79 171.45 0.00
## child_utt_cnt 9 1488 26.38 22.61 21.00 23.41 19.27 0.00
## conv_turn_count 10 1488 7.22 7.08 5.00 6.16 5.93 0.00
## child_cry_scnds 11 1488 14.90 21.62 8.53 10.61 10.12 0.00
## max range skew kurtosis se
## id_uni 13001.00 12200.00 -2.40 4.42 82.49
## age 14.23 3.48 1.01 0.48 0.02
## sit_time 1.00 1.00 0.37 -0.70 0.01
## held_time 0.99 0.99 3.23 12.79 0.00
## prone_time 1.00 1.00 2.27 6.77 0.00
## supine_time 1.00 1.00 2.24 4.38 0.01
## upright_time 1.00 1.00 1.27 1.17 0.01
## adult_word_cnt 1567.23 1567.23 1.83 5.16 5.55
## child_utt_cnt 121.00 121.00 1.22 1.40 0.59
## conv_turn_count 39.00 39.00 1.28 1.39 0.18
## child_cry_scnds 240.28 240.28 4.09 25.92 0.56age_round <- my_data$age
round(age_round)## [1] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [25] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [49] 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [73] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [97] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [121] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [145] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [169] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [193] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [217] 11 11 11 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [241] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [265] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [289] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [313] 13 13 13 13 13 13 13 13 13 13 13 13 13 14 14 14 14 14 14 14 14 14 14 14
## [337] 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14
## [361] 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14
## [385] 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 11 11 11
## [409] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [433] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [457] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [481] 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [505] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [529] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 11 11
## [553] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [577] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [601] 11 11 11 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [625] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [649] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [673] 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [697] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [721] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [745] 13 13 13 14 14 14 14 14 14 14 14 14 14 14 14 14 12 12 12 12 12 12 12 12
## [769] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [793] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [817] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [841] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 11
## [865] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [889] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [913] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [937] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [961] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [985] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1009] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1033] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1057] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12
## [1081] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1105] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1129] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1153] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1177] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 11 11 11 11 11 11
## [1201] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1225] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1249] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1273] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1297] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1321] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 12 12 12 12 12 12
## [1345] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1369] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1393] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [1417] 12 12 12 12 12 12 12 12 12 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1441] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [1465] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11Varying Intercept Model
imu_lm <- lmer(child_utt_cnt ~ 1 + (1 | age_round), data = my_data)
display(imu_lm)## lmer(formula = child_utt_cnt ~ 1 + (1 | age_round), data = my_data)
## coef.est coef.se
## 27.74 2.46
##
## Error terms:
## Groups Name Std.Dev.
## age_round (Intercept) 10.12
## Residual 20.51
## ---
## number of obs: 1488, groups: age_round, 18
## AIC = 13268.2, DIC = 13269.4
## deviance = 13265.8Varying Intercept Model with 1 predictor
imu_lm_2 <- lmer(child_utt_cnt ~ 1 + adult_word_cnt + (1 | age_round), data = my_data)
display(imu_lm_2)## lmer(formula = child_utt_cnt ~ 1 + adult_word_cnt + (1 | age_round),
## data = my_data)
## coef.est coef.se
## (Intercept) 23.62 2.48
## adult_word_cnt 0.02 0.00
##
## Error terms:
## Groups Name Std.Dev.
## age_round (Intercept) 9.81
## Residual 20.27
## ---
## number of obs: 1488, groups: age_round, 18
## AIC = 13243.1, DIC = 13222.2
## deviance = 13228.7imu_c <- coef(imu_lm_2)
imu_c## $age_round
## (Intercept) adult_word_cnt
## 10.7433264887064 15.914892 0.01645784
## 10.8090349075975 18.018780 0.01645784
## 10.8418891170431 37.038038 0.01645784
## 10.9075975359343 21.801115 0.01645784
## 10.9404517453799 26.140723 0.01645784
## 10.9733059548255 33.709314 0.01645784
## 11.0390143737166 25.599561 0.01645784
## 11.1047227926078 8.288907 0.01645784
## 11.8275154004107 37.417691 0.01645784
## 11.8932238193018 27.176697 0.01645784
## 11.9260780287474 9.903938 0.01645784
## 12.0246406570842 18.527623 0.01645784
## 12.1560574948665 14.274815 0.01645784
## 12.2546201232033 28.249594 0.01645784
## 12.8459958932238 12.024857 0.01645784
## 12.9774127310062 34.028956 0.01645784
## 14.0616016427105 35.548874 0.01645784
## 14.2258726899384 21.491169 0.01645784
##
## attr(,"class")
## [1] "coef.mer"the estimated regression line at 10.74 months is 15.91 + 0.016X, at 10.81 months its 18.02 + 0.016X, and so on…
imu_cF <- fixef(imu_lm_2)
imu_cF## (Intercept) adult_word_cnt
## 23.61975245 0.01645784the estimated regression line in average age group? is y = 23.62 - 0.016x
ranef(imu_lm_2)## $age_round
## (Intercept)
## 10.7433264887064 -7.704860
## 10.8090349075975 -5.600972
## 10.8418891170431 13.418286
## 10.9075975359343 -1.818637
## 10.9404517453799 2.520970
## 10.9733059548255 10.089562
## 11.0390143737166 1.979809
## 11.1047227926078 -15.330845
## 11.8275154004107 13.797939
## 11.8932238193018 3.556944
## 11.9260780287474 -13.715814
## 12.0246406570842 -5.092130
## 12.1560574948665 -9.344937
## 12.2546201232033 4.629842
## 12.8459958932238 -11.594896
## 12.9774127310062 10.409203
## 14.0616016427105 11.929121
## 14.2258726899384 -2.128583
##
## with conditional variances for "age_round"se.ranef(imu_lm_2)## $age_round
## (Intercept)
## 10.7433264887064 2.471381
## 10.8090349075975 2.678487
## 10.8418891170431 2.276636
## 10.9075975359343 2.221279
## 10.9404517453799 1.571933
## 10.9733059548255 2.055209
## 11.0390143737166 2.276636
## 11.1047227926078 2.221279
## 11.8275154004107 2.320981
## 11.8932238193018 2.418047
## 11.9260780287474 1.789700
## 12.0246406570842 1.957065
## 12.1560574948665 2.121691
## 12.2546201232033 2.076671
## 12.8459958932238 1.930259
## 12.9774127310062 2.401022
## 14.0616016427105 4.878095
## 14.2258726899384 2.208059anova(imu_lm, imu_lm_2)## refitting model(s) with ML (instead of REML)## Data: my_data
## Models:
## imu_lm: child_utt_cnt ~ 1 + (1 | age_round)
## imu_lm_2: child_utt_cnt ~ 1 + adult_word_cnt + (1 | age_round)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## imu_lm 3 13272 13288 -6632.9 13266
## imu_lm_2 4 13237 13258 -6614.3 13229 37.167 1 1.084e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(imu_lm_2)