Final
Hailey Rousey
2023-06-12
library(ggplot2)
library(interactions)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.1 ✔ tidyr 1.3.0
## ── 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 errorslibrary(tidyr)
library(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(psych)##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(MASS)##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(arm)##
## 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, simlibrary(readxl)beauty_raw <- read.csv("http://www.stat.columbia.edu/~gelman/arm/examples/beauty/ProfEvaltnsBeautyPublic.csv")beauty <- beauty_raw
dim(beauty)## [1] 463 64table(beauty$profnumber)##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 4 3 2 8 6 7 5 7 7 10 3 5 7 4 4 6 5 8 9 10 6 1 5 7 2 5
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## 7 4 4 1 7 3 5 13 3 4 8 3 8 1 5 4 4 3 4 3 1 3 7 13 6 4
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## 7 6 3 5 2 10 2 3 1 1 2 3 7 6 2 3 1 9 10 6 5 4 2 2 6 4
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
## 3 4 4 11 5 5 8 3 2 7 3 2 3 7 6 4ggplot(beauty, aes(y=courseevaluation,
x=factor(female, labels=c("Male",
"Female")))) + geom_boxplot(outlier.shape = NULL) +
geom_jitter(width=.1, height=.05, alpha=.5) +
xlab("Gender") + ylab("Course Evaluation")
Descriptive Statistics
describe(beauty)## vars n mean sd median trimmed mad min max range
## tenured 1 463 0.55 0.50 1.00 0.56 0.00 0.00 1.00 1.00
## profnumber 2 463 45.43 27.51 44.00 44.96 37.06 1.00 94.00 93.00
## minority 3 463 0.14 0.35 0.00 0.05 0.00 0.00 1.00 1.00
## age 4 463 48.37 9.80 48.00 48.35 11.86 29.00 73.00 44.00
## beautyf2upper 5 463 5.21 2.02 5.00 5.16 1.48 1.00 10.00 9.00
## beautyflowerdiv 6 463 3.96 1.87 4.00 3.86 1.48 1.00 8.00 7.00
## beautyfupperdiv 7 463 5.02 1.93 5.00 5.01 1.48 1.00 9.00 8.00
## beautym2upper 8 463 4.75 1.58 5.00 4.73 1.48 1.00 9.00 8.00
## beautymlowerdiv 9 463 3.41 1.64 3.00 3.33 1.48 1.00 7.00 6.00
## beautymupperdiv 10 463 4.15 2.11 4.00 4.00 1.48 1.00 9.00 8.00
## btystdave 11 463 -0.09 0.79 -0.16 -0.13 0.87 -1.54 1.88 3.42
## btystdf2u 12 463 -0.09 0.96 -0.19 -0.11 0.71 -2.10 2.20 4.29
## btystdfl 13 463 -0.09 1.00 -0.07 -0.15 0.79 -1.67 2.05 3.72
## btystdfu 14 463 -0.08 0.94 -0.09 -0.08 0.72 -2.03 1.84 3.87
## btystdm2u 15 463 -0.07 0.95 0.08 -0.08 0.90 -2.34 2.50 4.84
## btystdml 16 463 -0.07 0.96 -0.31 -0.12 0.87 -1.49 2.04 3.53
## btystdmu 17 463 -0.13 0.97 -0.20 -0.19 0.68 -1.57 2.10 3.67
## class1 18 463 0.01 0.10 0.00 0.00 0.00 0.00 1.00 1.00
## class2 19 463 0.00 0.07 0.00 0.00 0.00 0.00 1.00 1.00
## class3 20 463 0.02 0.13 0.00 0.00 0.00 0.00 1.00 1.00
## class4 21 463 0.04 0.20 0.00 0.00 0.00 0.00 1.00 1.00
## class5 22 463 0.01 0.09 0.00 0.00 0.00 0.00 1.00 1.00
## class6 23 463 0.01 0.11 0.00 0.00 0.00 0.00 1.00 1.00
## class7 24 463 0.01 0.09 0.00 0.00 0.00 0.00 1.00 1.00
## class8 25 463 0.00 0.07 0.00 0.00 0.00 0.00 1.00 1.00
## class9 26 463 0.02 0.13 0.00 0.00 0.00 0.00 1.00 1.00
## class10 27 463 0.01 0.10 0.00 0.00 0.00 0.00 1.00 1.00
## class11 28 463 0.00 0.07 0.00 0.00 0.00 0.00 1.00 1.00
## class12 29 463 0.01 0.08 0.00 0.00 0.00 0.00 1.00 1.00
## class13 30 463 0.01 0.08 0.00 0.00 0.00 0.00 1.00 1.00
## class14 31 463 0.01 0.08 0.00 0.00 0.00 0.00 1.00 1.00
## class15 32 463 0.00 0.07 0.00 0.00 0.00 0.00 1.00 1.00
## class16 33 463 0.01 0.09 0.00 0.00 0.00 0.00 1.00 1.00
## class17 34 463 0.02 0.12 0.00 0.00 0.00 0.00 1.00 1.00
## class18 35 463 0.01 0.09 0.00 0.00 0.00 0.00 1.00 1.00
## class19 36 463 0.01 0.11 0.00 0.00 0.00 0.00 1.00 1.00
## class20 37 463 0.01 0.10 0.00 0.00 0.00 0.00 1.00 1.00
## class21 38 463 0.03 0.17 0.00 0.00 0.00 0.00 1.00 1.00
## class22 39 463 0.02 0.15 0.00 0.00 0.00 0.00 1.00 1.00
## class23 40 463 0.01 0.10 0.00 0.00 0.00 0.00 1.00 1.00
## class24 41 463 0.01 0.08 0.00 0.00 0.00 0.00 1.00 1.00
## class25 42 463 0.01 0.08 0.00 0.00 0.00 0.00 1.00 1.00
## class26 43 463 0.01 0.08 0.00 0.00 0.00 0.00 1.00 1.00
## class27 44 463 0.00 0.07 0.00 0.00 0.00 0.00 1.00 1.00
## class28 45 463 0.01 0.09 0.00 0.00 0.00 0.00 1.00 1.00
## class29 46 463 0.00 0.07 0.00 0.00 0.00 0.00 1.00 1.00
## class30 47 463 0.02 0.13 0.00 0.00 0.00 0.00 1.00 1.00
## courseevaluation 48 463 4.00 0.55 4.00 4.03 0.59 2.10 5.00 2.90
## didevaluation 49 463 36.62 45.02 23.00 27.54 14.83 5.00 380.00 375.00
## female 50 463 0.42 0.49 0.00 0.40 0.00 0.00 1.00 1.00
## formal 51 463 0.17 0.37 0.00 0.08 0.00 0.00 1.00 1.00
## fulldept 52 463 0.89 0.31 1.00 0.99 0.00 0.00 1.00 1.00
## lower 53 463 0.34 0.47 0.00 0.30 0.00 0.00 1.00 1.00
## multipleclass 54 463 0.34 0.47 0.00 0.30 0.00 0.00 1.00 1.00
## nonenglish 55 463 0.06 0.24 0.00 0.00 0.00 0.00 1.00 1.00
## onecredit 56 463 0.06 0.23 0.00 0.00 0.00 0.00 1.00 1.00
## percentevaluating 57 463 74.43 16.76 76.92 75.64 17.08 10.42 100.00 89.58
## profevaluation 58 463 4.17 0.54 4.30 4.22 0.59 2.30 5.00 2.70
## students 59 463 55.18 75.07 29.00 38.66 19.27 8.00 581.00 573.00
## tenuretrack 60 463 0.78 0.41 1.00 0.85 0.00 0.00 1.00 1.00
## blkandwhite 61 463 0.17 0.37 0.00 0.09 0.00 0.00 1.00 1.00
## btystdvariance 62 463 1.84 1.26 1.57 1.71 1.12 0.09 5.79 5.71
## btystdavepos 63 463 0.28 0.48 0.00 0.17 0.00 0.00 1.88 1.88
## btystdaveneg 64 463 -0.37 0.43 -0.16 -0.31 0.23 -1.54 0.00 1.54
## skew kurtosis se
## tenured -0.19 -1.97 0.02
## profnumber 0.12 -1.25 1.28
## minority 2.09 2.37 0.02
## age 0.05 -0.80 0.46
## beautyf2upper 0.23 -0.27 0.09
## beautyflowerdiv 0.37 -0.57 0.09
## beautyfupperdiv 0.02 -0.78 0.09
## beautym2upper 0.16 -0.23 0.07
## beautymlowerdiv 0.43 -0.87 0.08
## beautymupperdiv 0.56 -0.40 0.10
## btystdave 0.51 -0.40 0.04
## btystdf2u 0.23 -0.27 0.04
## btystdfl 0.37 -0.57 0.05
## btystdfu 0.03 -0.79 0.04
## btystdm2u 0.16 -0.23 0.04
## btystdml 0.43 -0.87 0.04
## btystdmu 0.56 -0.40 0.05
## class1 9.44 87.22 0.00
## class2 15.07 225.51 0.00
## class3 7.38 52.65 0.01
## class4 4.61 19.31 0.01
## class5 10.58 110.27 0.00
## class6 8.58 71.86 0.01
## class7 10.58 110.27 0.00
## class8 15.07 225.51 0.00
## class9 7.38 52.65 0.01
## class10 9.44 87.22 0.00
## class11 15.07 225.51 0.00
## class12 12.26 148.68 0.00
## class13 12.26 148.68 0.00
## class14 12.26 148.68 0.00
## class15 15.07 225.51 0.00
## class16 10.58 110.27 0.00
## class17 7.92 60.88 0.01
## class18 10.58 110.27 0.00
## class19 8.58 71.86 0.01
## class20 9.44 87.22 0.00
## class21 5.47 27.97 0.01
## class22 6.23 36.94 0.01
## class23 9.44 87.22 0.00
## class24 12.26 148.68 0.00
## class25 12.26 148.68 0.00
## class26 12.26 148.68 0.00
## class27 15.07 225.51 0.00
## class28 10.58 110.27 0.00
## class29 15.07 225.51 0.00
## class30 7.38 52.65 0.01
## courseevaluation -0.46 -0.13 0.03
## didevaluation 4.47 25.98 2.09
## female 0.32 -1.90 0.02
## formal 1.79 1.19 0.02
## fulldept -2.55 4.53 0.01
## lower 0.68 -1.54 0.02
## multipleclass 0.68 -1.54 0.02
## nonenglish 3.68 11.54 0.01
## onecredit 3.76 12.14 0.01
## percentevaluating -0.66 0.04 0.78
## profevaluation -0.70 0.04 0.03
## students 4.13 21.74 3.49
## tenuretrack -1.35 -0.19 0.02
## blkandwhite 1.77 1.12 0.02
## btystdvariance 0.92 0.13 0.06
## btystdavepos 1.70 1.84 0.02
## btystdaveneg -0.89 -0.37 0.02- Here are the actual items to be answered.
- [20 pts.] Write a varying-intercept model for these data with no group-level predictors. Fit this model using lmer() and interpret the results.
model_1 <- lmer(courseevaluation ~ 1 + (1|profnumber), data = beauty)
display(model_1)## lmer(formula = courseevaluation ~ 1 + (1 | profnumber), data = beauty)
## coef.est coef.se
## 3.94 0.05
##
## Error terms:
## Groups Name Std.Dev.
## profnumber (Intercept) 0.38
## Residual 0.41
## ---
## number of obs: 463, groups: profnumber, 94
## AIC = 650, DIC = 635.3
## deviance = 639.7#fixef(model_1)
#ranef(model_1)
#coef(model_1)Interpretation: The coef. est of 3.94 is the fixed intercept representing the average course evaluation for all the professors. The coef. se of .05 is the standard error of the coef. est. For the intercepts that vary across professors there is a standard deviation of .38. We can see the intercepts for each professor by using the coef function.
- [20 pts.] Write a varying-intercept model that you would like to fit including three group-level predictors. Fit this model using lmer() and interpret the results.
model_2 <- lmer(courseevaluation ~ female + age + btystdave + (1 | profnumber), data = beauty)
display(model_2)## lmer(formula = courseevaluation ~ female + age + btystdave +
## (1 | profnumber), data = beauty)
## coef.est coef.se
## (Intercept) 4.06 0.24
## female -0.21 0.09
## age 0.00 0.00
## btystdave 0.13 0.06
##
## Error terms:
## Groups Name Std.Dev.
## profnumber (Intercept) 0.37
## Residual 0.41
## ---
## number of obs: 463, groups: profnumber, 94
## AIC = 662.1, DIC = 609.2
## deviance = 629.7#fixef(model_2)
#ranef(model_2)
#coef(model_2)Interpretation: When all group level predictors are 0 the fixed intercept or average course evaluation is 4.06. Assuming 1 is female and 0 is male, being female is associated with a .21 decrease in course evaluation, age apparently had no impact, and for every one unit increase in average estimates of attractiveness there is a .13 increase in course evaluation. Again, we can see the intercepts for each professor by using the coef function. For these intercepts that vary across professors there is a standard deviation of .37.
- [10 pts.] How does the variation in average ratings across instructors compare to the variation in ratings across the different classes for the same instructor?
Model 1 - Across instructors the variation is .38 and for the same instructor across different classes the variation is .41. The variance ratio is:
.38^2 / .41^2 ## [1] 0.8590125Model 2 - Across instructors the variation is .37 and for the same instructor across different classes the variation is .41. The variance ratio is:
.37^2 / .41^2 ## [1] 0.8143962ICC
icc_by_hand_model_1 <-.38^2 / (.38^2 + .41^2)
icc_by_hand_model_1## [1] 0.46208icc_by_hand_model_2 <-.37^2 / (.37^2 + .41^2)
icc_by_hand_model_2## [1] 0.4488525This gives us an idea of how much information the grouping is telling us. If the ICC is 0 then grouping is not providing any additional information. However, if it is 1 then the individual level is not providing information and we might as well just look at group level data. In our case we are getting reasonable information from the group level that we don’t get from individual level.
- Fit a multilevel model to predict course evaluations from beauty and other predictors in the beauty dataset allowing the intercept and coefficient for beauty to vary by professor:
- [15 pts.] Write the model in statistical notation.
Individual level:
y_i ~ N (alpha_{j[i]} + beta_{j[i]}x_i, sigma_y^2) for i = 1,…, n
OR
courseevaluation_ij = β0j + β1j * btystdave_ij + β2 * age_ij + β3 * female_ij + u0j + u1j * btystdave_ij + error_ij
Group level: alpha_j ~ N (( mu_alpha),(sigma_alpha^2 rho sigma_alpha sigma_beta)) beta_j ~ N ((mu_beta),(rho sigma_alpha sigma_beta, sigma_beta^2))
I tried my best to format it like the book but I’m not sure how to properly do it using R. Th above should be all within the same parentheses rather than in two separate.
- [20 pts.] Fit the model using lmer() and discuss the results: the coefficient estimates and the estimated standard deviation and correlation parameters. Identify each of the estimated parameters with the notation in your model from (a).
model_3 <- lmer(courseevaluation ~ btystdave + age + female + (1 + btystdave | profnumber), data = beauty)
display(model_3)## lmer(formula = courseevaluation ~ btystdave + age + female +
## (1 + btystdave | profnumber), data = beauty)
## coef.est coef.se
## (Intercept) 4.03 0.24
## btystdave 0.13 0.06
## age 0.00 0.00
## female -0.22 0.09
##
## Error terms:
## Groups Name Std.Dev. Corr
## profnumber (Intercept) 0.36
## btystdave 0.10 -0.30
## Residual 0.41
## ---
## number of obs: 463, groups: profnumber, 94
## AIC = 665.6, DIC = 608.6
## deviance = 629.1Interpretation: courseevaluation_ij = 4.03 + 0.13 * btystdave_ij + 0 * age_ij + -0.22 * female_ij + u0j + u1j * btystdave_ij + error_ij
When all predictors are 0 the average course evaluation is 4.03. For every one unit increase in average attractiveness score there is a .13 increase in course evaluation. For every one unit increase in age there is no change in course evaluation. Being female is associated with a .22 decrease in course evaluation. (I am assuming female = 1 and male = 0)
The unexplained within-group variation has an estimated standard deviation of .41. The estimated standard deviation of the intercepts is .36. The estimated standard deviation of the slopes is .10. The estimated correlation between the intercepts and slopes is -.30.
#fixef(model_3)
#ranef(model_3)
#coef(model_3)- [15 pts.] Display the estimated model graphically in plots that also include the data.
ggplot(beauty, aes(x = btystdave, y = courseevaluation, color = profnumber)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, linetype = "dashed", formula = y ~ x, aes(group = profnumber), size = 1) +
labs(x = "Beauty", y = "Evaluation") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot(beauty, aes(x = age, y = courseevaluation, color = profnumber)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, linetype = "dashed", formula = y ~ x, aes(group = profnumber), size = 1) +
labs(x = "Age", y = "Evaluation") +
theme_minimal()
ggplot(beauty, aes(x = female, y = courseevaluation, color = profnumber)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, linetype = "dashed", formula = y ~ x, aes(group = profnumber), size = 1) +
labs(x = "Gender", y = "Evaluation") +
theme_minimal()