if (Sys.info()["sysname"] == "Windows") {
setwd("~/Masters/DATA606/Week8/Lab/Lab8")
} else {
setwd("~/Documents/Masters/DATA606/Week8/Lab/Lab8")
}
load("more/evals.RData")
require(ggplot2)
## Loading required package: ggplot2
Answer:
This is an observational study. No, since this is an observational study it is only generally sufficient to show associations between variables and should not be used to draw causal conclusions. The appropriate way to rephrase the question would be to ask: Is there a correlation between the instructor’s beauty and evaluation scores?
Answer:
ggplot(evals, aes(x = score)) + geom_histogram(binwidth = 0.125,
position = "identity", aes(y = ..density..)) + stat_function(fun = dnorm,
color = "black", args = list(mean = mean(evals$score), sd(evals$score))) +
labs(x = "Score")
qqnorm(evals$score)
qqline(evals$score)
The distribution has a left skew. Additionally, we can see that the residuals have a steped nature which show that there are relatively few possible unique values that were recorded. This shows that students will typically give a high score to an average class since the median is much closer to the maximum value than the minimum. It leads me to assume that if a student felt neutral towards a teacher’s performance, they generally give them a good score, possibly giving the professor the benefit of the doubt about their skills. I did expect to see something similar to this. In my experience, I find that people generally will lean towards being aggreeable in these type of surveys since it doesn’t require an argument and at least, in their minds, they would not have to strongly defend a higher score rather than a lower.
Answer:
for this exercise, I will compare the percent of students who completed the evaluations against the average beauty of the professor.
ggplot(evals, aes(x = cls_perc_eval, y = bty_avg)) + geom_point() +
geom_smooth(method = lm, fullrange = TRUE)
e3.lm <- lm(bty_avg ~ cls_perc_eval, data = evals)
summary(e3.lm)
##
## Call:
## lm(formula = bty_avg ~ cls_perc_eval, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7153 -1.2138 -0.1994 1.0002 3.9876
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.381014 0.320063 10.56 < 2e-16 ***
## cls_perc_eval 0.013931 0.004196 3.32 0.00097 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.511 on 461 degrees of freedom
## Multiple R-squared: 0.02336, Adjusted R-squared: 0.02124
## F-statistic: 11.02 on 1 and 461 DF, p-value: 0.0009704
From the linear regression line, there appears to be a slightly positive slope between the percent of students who completed the evaluation and the average beauty rating of the professor. However, it does not appear that there is a good r-squared correlation with a linear model and there are not any apparent patterns for a better fit in other models (e.g. quadratic, power, etc.). It is noted that there is a significantly low p-value for the slope of the fit which suggests that the class percentage may be a statistically significant predictor of average beauty even though the correlation is low.
Answer:
plot(evals$score ~ evals$bty_avg)
jitter_bty_avg <- jitter(evals$bty_avg, factor = 1.5)
plot(evals$score ~ jitter_bty_avg)
Due to the low number of unique variable values possible, there appear to be many overlapping data points. The jitter plot slightly shifts the vector values so that this can be more apparent. The initial scatterplot is misleading because it does not accurately depict the number of points of data.
Answer:
m_bty <- lm(score ~ bty_avg, data = evals)
summary(m_bty)
##
## Call:
## lm(formula = score ~ bty_avg, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9246 -0.3690 0.1420 0.3977 0.9309
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.88034 0.07614 50.96 < 2e-16 ***
## bty_avg 0.06664 0.01629 4.09 5.08e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5348 on 461 degrees of freedom
## Multiple R-squared: 0.03502, Adjusted R-squared: 0.03293
## F-statistic: 16.73 on 1 and 461 DF, p-value: 5.083e-05
plot(evals$score ~ jitter_bty_avg)
abline(m_bty)
\[evaluation\quad score\quad =\quad 3.88\quad +\quad 0.067\quad *\quad average\quad beauty\]
For every one point of increase in average beauty, a professor is expected to recieve an additional 0.067 points in their evaluation. Since the p-value is less than 0.05, it appears that average beauty is a statistically significant predictor of the evaluation score. It does not appear to be a practically significant predictor because the r-squared value is 0.03 which means that only 3% of the variability is explained by the linear model. Additionally, the slope of the equation is very small and may not have a significant impact on the results.
Answer:
hist(m_bty$residuals)
qqnorm(m_bty$residuals)
qqline(m_bty$residuals)
ggplot(m_bty, aes(y = abs(resid(m_bty)), x = m_bty$fitted.values)) +
geom_point() + geom_hline(yintercept = 0)
ggplot(m_bty, aes(y = resid(m_bty), x = bty_avg)) + geom_point() +
geom_hline(yintercept = 0)
The assumptions to check are the following:
The histogram and normal probability plot show a left skew which means that this linear model may not be a good fit for the data.
The plot of the absolute value of the residuals against the fitted value show that the variability is nearly constant.
We do not have the order of the collection data; however, since this was a published study we can assume that the residuals are independent.
there does not appear to be a trend in the residuals in the graph that shows the actual values (not absolute) that would suggest any other type of model would be a better fit.
Answer:
plot(evals$bty_avg ~ evals$bty_f1lower)
cor(evals$bty_avg, evals$bty_f1lower)
## [1] 0.8439112
plot(evals[, 13:19])
m_bty_gen <- lm(score ~ bty_avg + gender, data = evals)
summary(m_bty_gen)
##
## Call:
## lm(formula = score ~ bty_avg + gender, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8305 -0.3625 0.1055 0.4213 0.9314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.74734 0.08466 44.266 < 2e-16 ***
## bty_avg 0.07416 0.01625 4.563 6.48e-06 ***
## gendermale 0.17239 0.05022 3.433 0.000652 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5287 on 460 degrees of freedom
## Multiple R-squared: 0.05912, Adjusted R-squared: 0.05503
## F-statistic: 14.45 on 2 and 460 DF, p-value: 8.177e-07
The assumptions to check are the following:
hist(m_bty_gen$residuals)
qqnorm(m_bty_gen$residuals)
qqline(m_bty_gen$residuals)
The residuals have a left skew; therefore they are not nearly normal and a linear model may not be an appropriate fit for the data.
ggplot(m_bty_gen, aes(y = abs(resid(m_bty_gen)), x = m_bty_gen$fitted.values)) +
geom_point() + geom_hline(yintercept = 0)
The comparison of the the fitted values against the absolute value of the residuals shows a constant variability and there does not appear to be any underlying patterns.
As with the previous response for this assumption, we can assume independence.
ggplot(m_bty_gen, aes(y = resid(m_bty_gen), x = bty_avg)) + geom_point() +
geom_hline(yintercept = 0)
ggplot(m_bty_gen, aes(y = resid(m_bty_gen), x = evals$gender)) +
geom_boxplot() + geom_point()
To check this assumption, we plot the residuals against each variable. There does not appear to be any pattern that would suggest a linear fit is not appropriate, especially for the gender data which only has two categorical variables.
Answer:
summary(m_bty_gen)
##
## Call:
## lm(formula = score ~ bty_avg + gender, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8305 -0.3625 0.1055 0.4213 0.9314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.74734 0.08466 44.266 < 2e-16 ***
## bty_avg 0.07416 0.01625 4.563 6.48e-06 ***
## gendermale 0.17239 0.05022 3.433 0.000652 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5287 on 460 degrees of freedom
## Multiple R-squared: 0.05912, Adjusted R-squared: 0.05503
## F-statistic: 14.45 on 2 and 460 DF, p-value: 8.177e-07
summary(m_bty)
##
## Call:
## lm(formula = score ~ bty_avg, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9246 -0.3690 0.1420 0.3977 0.9309
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.88034 0.07614 50.96 < 2e-16 ***
## bty_avg 0.06664 0.01629 4.09 5.08e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5348 on 461 degrees of freedom
## Multiple R-squared: 0.03502, Adjusted R-squared: 0.03293
## F-statistic: 16.73 on 1 and 461 DF, p-value: 5.083e-05
Yes, since the p-value for the slope of the bty_avg variable is less than 0.05, it is still a significant predictor of the score. However, when compared to the linear model that just uses bty_avg as a predictor, we can see that the addition of gender has influenced the parameter estimate and the significance of bty_avg.
Answer:
The equation of the line for males is the following:
\[\widehat {evaluation\quad score\quad} =\quad 3.747\quad +\quad 0.074\quad *\quad average\quad beauty\quad +\quad 0.172\quad *\quad gender\\ \left( evaluation\quad score|\quad male \right) \quad =\quad 3.747\quad +\quad 0.074\quad *\quad average\quad beauty\quad +\quad 0.172\quad *\quad 1\\ \left( evaluation\quad score|\quad male \right) \quad =\quad 3.919\quad +\quad 0.074\quad *\quad average\quad beauty\]
For two professors with the same beauty rating, the male professor tends to have the higher course evaluation score.
Answer:
m_bty_rank <- lm(score ~ bty_avg + rank, data = evals)
summary(m_bty_rank)
##
## Call:
## lm(formula = score ~ bty_avg + rank, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8713 -0.3642 0.1489 0.4103 0.9525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.98155 0.09078 43.860 < 2e-16 ***
## bty_avg 0.06783 0.01655 4.098 4.92e-05 ***
## ranktenure track -0.16070 0.07395 -2.173 0.0303 *
## ranktenured -0.12623 0.06266 -2.014 0.0445 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5328 on 459 degrees of freedom
## Multiple R-squared: 0.04652, Adjusted R-squared: 0.04029
## F-statistic: 7.465 on 3 and 459 DF, p-value: 6.88e-05
It appears to select one level as a reference which is then compared to the remaining levels in the categorical set.
Answer:
I would think that the language of the school where the professor recieved their education would have the highest p-value since they obviously have shown enough proficiency with the necessary language (assuming english) to give their own lectures.
Answer:
m_full <- lm(score ~ rank + ethnicity + gender + language + age +
cls_perc_eval + cls_students + cls_level + cls_profs + cls_credits +
bty_avg + pic_outfit + pic_color, data = evals)
summary(m_full)
##
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age +
## cls_perc_eval + cls_students + cls_level + cls_profs + cls_credits +
## bty_avg + pic_outfit + pic_color, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.77397 -0.32432 0.09067 0.35183 0.95036
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0952141 0.2905277 14.096 < 2e-16 ***
## ranktenure track -0.1475932 0.0820671 -1.798 0.07278 .
## ranktenured -0.0973378 0.0663296 -1.467 0.14295
## ethnicitynot minority 0.1234929 0.0786273 1.571 0.11698
## gendermale 0.2109481 0.0518230 4.071 5.54e-05 ***
## languagenon-english -0.2298112 0.1113754 -2.063 0.03965 *
## age -0.0090072 0.0031359 -2.872 0.00427 **
## cls_perc_eval 0.0053272 0.0015393 3.461 0.00059 ***
## cls_students 0.0004546 0.0003774 1.205 0.22896
## cls_levelupper 0.0605140 0.0575617 1.051 0.29369
## cls_profssingle -0.0146619 0.0519885 -0.282 0.77806
## cls_creditsone credit 0.5020432 0.1159388 4.330 1.84e-05 ***
## bty_avg 0.0400333 0.0175064 2.287 0.02267 *
## pic_outfitnot formal -0.1126817 0.0738800 -1.525 0.12792
## pic_colorcolor -0.2172630 0.0715021 -3.039 0.00252 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.498 on 448 degrees of freedom
## Multiple R-squared: 0.1871, Adjusted R-squared: 0.1617
## F-statistic: 7.366 on 14 and 448 DF, p-value: 6.552e-14
My suspicions were not confirmed, the p-value for the language variable is 0.0397, which is less than 0.05, showing that it is a significant predictor. It appears that the cls_profs variable is the least significant predictor.
Answer:
The p-value for the ethnicity variable is 0.117 which shows that is not a significant predictor of the evaluation score. Additionally, the 95% confidence interval for the predictor variable is the following, which shows that the value of 0 is within the range:
slope <- 0.1235
se <- 0.0786
t_value <- 1.571
t_95 <- qt(p = 0.975, df = (nrow(evals) - 1))
slope - t_95 * se
## [1] -0.0309578
slope + t_95 * se
## [1] 0.2779578
Answer:
m_full_1 <- lm(score ~ rank + ethnicity + gender + language +
age + cls_perc_eval + cls_students + cls_level + cls_credits +
bty_avg + pic_outfit + pic_color, data = evals)
summary(m_full_1)
##
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age +
## cls_perc_eval + cls_students + cls_level + cls_credits +
## bty_avg + pic_outfit + pic_color, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7836 -0.3257 0.0859 0.3513 0.9551
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0872523 0.2888562 14.150 < 2e-16 ***
## ranktenure track -0.1476746 0.0819824 -1.801 0.072327 .
## ranktenured -0.0973829 0.0662614 -1.470 0.142349
## ethnicitynot minority 0.1274458 0.0772887 1.649 0.099856 .
## gendermale 0.2101231 0.0516873 4.065 5.66e-05 ***
## languagenon-english -0.2282894 0.1111305 -2.054 0.040530 *
## age -0.0089992 0.0031326 -2.873 0.004262 **
## cls_perc_eval 0.0052888 0.0015317 3.453 0.000607 ***
## cls_students 0.0004687 0.0003737 1.254 0.210384
## cls_levelupper 0.0606374 0.0575010 1.055 0.292200
## cls_creditsone credit 0.5061196 0.1149163 4.404 1.33e-05 ***
## bty_avg 0.0398629 0.0174780 2.281 0.023032 *
## pic_outfitnot formal -0.1083227 0.0721711 -1.501 0.134080
## pic_colorcolor -0.2190527 0.0711469 -3.079 0.002205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4974 on 449 degrees of freedom
## Multiple R-squared: 0.187, Adjusted R-squared: 0.1634
## F-statistic: 7.943 on 13 and 449 DF, p-value: 2.336e-14
There was a slight change in all the variables. Since the overally change appeared to be small, it suggests that the dropped data point may have been collinear with one or more of the other variables since the exclusion of it did not affect the model.
Answer:
# Remove cls_level m_full_2 <- lm(score ~ rank + ethnicity +
# gender + language + age + cls_perc_eval + cls_students +
# cls_credits + bty_avg + pic_outfit + pic_color, data =
# evals) summary(m_full_2) Remove cls_students m_full_3 <-
# lm(score ~ rank + ethnicity + gender + language + age +
# cls_perc_eval + cls_credits + bty_avg + pic_outfit +
# pic_color, data = evals) summary(m_full_3) Remove rank
# m_full_4 <- lm(score ~ ethnicity + gender + language + age
# + cls_perc_eval + cls_credits + bty_avg + pic_outfit +
# pic_color, data = evals) summary(m_full_4) Remove
# pic_outfit
m_full_5 <- lm(score ~ ethnicity + gender + language + age +
cls_perc_eval + cls_credits + bty_avg + pic_color, data = evals)
summary(m_full_5)
##
## Call:
## lm(formula = score ~ ethnicity + gender + language + age + cls_perc_eval +
## cls_credits + bty_avg + pic_color, data = evals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85320 -0.32394 0.09984 0.37930 0.93610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.771922 0.232053 16.255 < 2e-16 ***
## ethnicitynot minority 0.167872 0.075275 2.230 0.02623 *
## gendermale 0.207112 0.050135 4.131 4.30e-05 ***
## languagenon-english -0.206178 0.103639 -1.989 0.04726 *
## age -0.006046 0.002612 -2.315 0.02108 *
## cls_perc_eval 0.004656 0.001435 3.244 0.00127 **
## cls_creditsone credit 0.505306 0.104119 4.853 1.67e-06 ***
## bty_avg 0.051069 0.016934 3.016 0.00271 **
## pic_colorcolor -0.190579 0.067351 -2.830 0.00487 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4992 on 454 degrees of freedom
## Multiple R-squared: 0.1722, Adjusted R-squared: 0.1576
## F-statistic: 11.8 on 8 and 454 DF, p-value: 2.58e-15
Using backward selection with p-value as the criteria, I eliminated cls_level, cls_students, rank, and pic_outfit. The resulting model has all predictor variables that are significant within a value of 0.05. The equation is the following:
\[\widehat {evaluation\quad score\quad} =\quad 3.77\quad +\quad ethnicity\quad *\quad 0.17\quad +\quad gender\quad *\quad 0.21\quad \\+\quad language\quad *\quad -0.21\quad +\quad age\quad *\quad -0.006\quad +\quad cls\_ perc\_ eval\quad *\quad 0.005\quad \\+\quad cls\_ credits\quad *\quad 0.51\quad +\quad bty\_ avg\quad *\quad 0.05\quad +\quad pic\_ color\quad *\quad -0.19\]
Answer:
The assumptions to check are the following:
hist(m_full_5$residuals)
qqnorm(m_full_5$residuals)
qqline(m_full_5$residuals)
The residuals appear to have a left skew so we should proceed with caution.
ggplot(m_full_5, aes(y = abs(resid(m_full_5)), x = m_full_5$fitted.values)) +
geom_point() + geom_hline(yintercept = 0)
The variability appears to be relatively constant; however, it appears there may be somewhat less variability at the two ends of the fitted values (<3.5 and >4.75).
As with the previous response for this assumption, we can assume independence.
# Ethnicity
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$ethnicity)) +
geom_boxplot() + geom_point()
# Gender
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$gender)) +
geom_boxplot() + geom_point()
# Language
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$language)) +
geom_boxplot() + geom_point()
# Age
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$age)) +
geom_point() + geom_hline(yintercept = 0)
# cls_perc_eval
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$cls_perc_eval)) +
geom_point() + geom_hline(yintercept = 0)
# cls_credits
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$cls_credits)) +
geom_boxplot() + geom_point()
# bty_avg
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$bty_avg)) +
geom_point() + geom_hline(yintercept = 0)
# pic_color
ggplot(data = m_full_5, aes(y = resid(m_full_5), x = evals$pic_color)) +
geom_boxplot() + geom_point()
Given that all the categorical variables chosen have two levels, they definitely appear to have linera relationships. Also, the numerical variables do not appear to have any underlying patterns that would suggest the data follows a non-linear correlation.
Answer:
It does not appear that this issue would affect the conditions of linear regression. It is not required that each case have a different professor. It is possible that some professors that teach many types of classes, or teach one type of class many times, would cause some issues with the model; however, as long as there arent’ a significant number of these types of issues, the model should not be greatly impacted.
Answer:
The characteristics of a professor and course that would be associated with a high evaluation score are as follows:
Answer:
No, universities can vary widely in the types of professor’s they employ, the prestige (difficulty) of some departments vs. others, and the types of high level classes they offer. I would think that in order to expand this study to other universities, then there should be some type of sampling done to try and normalize the data.