Lab 4: Multivariable Linear Regression

library(tidyverse)
library(openintro)
library(GGally)

## Warning: package 'GGally' was built under R version 4.3.3

library(broom)

Exercise 1

Being that there are no controls or experimental groups, this, in fact, is an observational study. Further, this is only an observational study; thus, there cannot be causation between the explanatory and response variables. Rather, there can only be correlation. In my opinion, it is possible to answer this question. From what I understand, we can surmise the instructor’s beauty has a positive correlation to the course evaluation. Note: depending on how you view this, it may also been seen as a negative correlation.

glimpse(evals)

## Rows: 463
## Columns: 23
## $ course_id     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
## $ prof_id       <int> 1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5,…
## $ score         <dbl> 4.7, 4.1, 3.9, 4.8, 4.6, 4.3, 2.8, 4.1, 3.4, 4.5, 3.8, 4…
## $ rank          <fct> tenure track, tenure track, tenure track, tenure track, …
## $ ethnicity     <fct> minority, minority, minority, minority, not minority, no…
## $ gender        <fct> female, female, female, female, male, male, male, male, …
## $ language      <fct> english, english, english, english, english, english, en…
## $ age           <int> 36, 36, 36, 36, 59, 59, 59, 51, 51, 40, 40, 40, 40, 40, …
## $ cls_perc_eval <dbl> 55.81395, 68.80000, 60.80000, 62.60163, 85.00000, 87.500…
## $ cls_did_eval  <int> 24, 86, 76, 77, 17, 35, 39, 55, 111, 40, 24, 24, 17, 14,…
## $ cls_students  <int> 43, 125, 125, 123, 20, 40, 44, 55, 195, 46, 27, 25, 20, …
## $ cls_level     <fct> upper, upper, upper, upper, upper, upper, upper, upper, …
## $ cls_profs     <fct> single, single, single, single, multiple, multiple, mult…
## $ cls_credits   <fct> multi credit, multi credit, multi credit, multi credit, …
## $ bty_f1lower   <int> 5, 5, 5, 5, 4, 4, 4, 5, 5, 2, 2, 2, 2, 2, 2, 2, 2, 7, 7,…
## $ bty_f1upper   <int> 7, 7, 7, 7, 4, 4, 4, 2, 2, 5, 5, 5, 5, 5, 5, 5, 5, 9, 9,…
## $ bty_f2upper   <int> 6, 6, 6, 6, 2, 2, 2, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 9, 9,…
## $ bty_m1lower   <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 7, 7,…
## $ bty_m1upper   <int> 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 6, 6,…
## $ bty_m2upper   <int> 6, 6, 6, 6, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 6, 6,…
## $ bty_avg       <dbl> 5.000, 5.000, 5.000, 5.000, 3.000, 3.000, 3.000, 3.333, …
## $ pic_outfit    <fct> not formal, not formal, not formal, not formal, not form…
## $ pic_color     <fct> color, color, color, color, color, color, color, color, …

?evals

## starting httpd help server ... done

library(ggplot2)

Exercise 2

Based on the distribution, it is skewed to the left. Students tend to have more positive evaluations than negative for their instructor. Broadly speaking, I would anticipate this. It would be more appropriate to assume a normal distribution where most instructors would be evaluated as “average” and fewer instructors would be evaluated in either extreme, unsatisfactory or excellent. Perhaps the students in this course did quite well.

hist(evals$score)

Exercise 3

Boxplot: there does not appear to be a relationship between the instructor’s age and beauty score. I think the inverse would make more sense. For example, the younger the instructor, the higher the beauty score. However, based on the median scores, this does not appear to be so.

Scatterplot: there appear to be more observations than the approximiateapproximate number of points on the scatterplot.

ggplot(data = evals, aes(x = bty_avg, y = score)) + geom_point()

boxplot(evals$bty_avg ~ evals$age)

plot (evals$score ~ evals$bty_avg)

nrow(evals)

## [1] 463

Exercise 4

The original scatterplot was not able to show the relationship between beauty average and score for the teacher due to the multiple ties–overlapping scores–that is just represented by a single circle on the scatterplot.

ggplot(data = evals, aes(x = bty_avg, y = score)) + geom_jitter()

ggplot(evals, aes(bty_avg, score)) + geom_point(position = position_jitter(w = 0.3, h = 0.3)) + ylab("score") + xlab("beauty average")

Exercise 5

Yes, “bty_avg” does, in fact, appear to be a practically significant predictor. Specifically, a significant predictor of evaluation socre with p-value close of 0. However, it may not be a practically significant predictor of evaluation score, though, since for every 1 point increase in “bty_avg,” the model can only predict an increase of 0.06664 which does not quite change the evaluation score.

Note equation: y = 3.88034 + 0.06664 * bty/avg

m_bty <- lm(evals$score ~ evals$bty_avg)
plot(jitter(evals$score,factor=1.2) ~ jitter(evals$bty_avg,factor=1.2))
abline(m_bty)

ggplot(data = evals, aes(x = bty_avg, y = score)) + geom_jitter() + geom_smooth(method = "lm")

## `geom_smooth()` using formula = 'y ~ x'

ggplot(data = evals, aes(x = bty_avg, y = score)) + geom_jitter() + geom_smooth(method = "lm", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

cor(evals$score, evals$bty_avg)

## [1] 0.1871424

summary(m_bty)

## 
## Call:
## lm(formula = evals$score ~ evals$bty_avg)
## 
## 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 ***
## evals$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

Exercise 6

ggplot(data = evals, aes(x = bty_f1lower, y = bty_avg)) + geom_point()

evals %>% summarise(cor(bty_avg, bty_f1lower))

## # A tibble: 1 × 1
##   `cor(bty_avg, bty_f1lower)`
##                         <dbl>
## 1                       0.844

plot(evals[,13:19])

evals %>% select(contains("bty"))%>% ggpairs()

m_bty_gen <- lm(score ~ bty_avg + gender, data = evals) 
tidy(m_bty_gen)

## # A tibble: 3 × 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)   3.75      0.0847     44.3  6.23e-168
## 2 bty_avg       0.0742    0.0163      4.56 6.48e-  6
## 3 gendermale    0.172     0.0502      3.43 6.52e-  4

names(m_bty_gen)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

plot_ss <- function(x, y, showSquares = FALSE, leastSquares = FALSE){
  plot(y~x, asp = 1)# xlab = paste(substitute(x)), ylab = paste(substitute(y)))
  
  if(leastSquares){
    m1 <- lm(y~x)
    y.hat <- m1$fit
  } else{
    cat("Click two points to make a line.")
    pt1 <- locator(1)
    points(pt1$x, pt1$y, pch = 4)
    pt2 <- locator(1)
    points(pt2$x, pt2$y, pch = 4)
    pts <- data.frame("x" = c(pt1$x, pt2$x),"y" = c(pt1$y, pt2$y))
    m1 <- lm(y ~ x, data = pts)
    y.hat <- predict(m1, newdata = data.frame(x))
  }
  r <- y - y.hat
  abline(m1)

  oSide <- x - r
  LLim <- par()$usr[1]
  RLim <- par()$usr[2]
  oSide[oSide < LLim | oSide > RLim] <- c(x + r)[oSide < LLim | oSide > RLim] # move boxes to avoid margins

  n <- length(y.hat)
  for(i in 1:n){
    lines(rep(x[i], 2), c(y[i], y.hat[i]), lty = 2, col = "blue")
    if(showSquares){
    lines(rep(oSide[i], 2), c(y[i], y.hat[i]), lty = 3, col = "orange")
    lines(c(oSide[i], x[i]), rep(y.hat[i],2), lty = 3, col = "orange")
    lines(c(oSide[i], x[i]), rep(y[i],2), lty = 3, col = "orange")
    }
  }

  SS <- round(sum(r^2), 3)
  cat("\r                                ")
  print(m1)
  cat("Sum of Squares: ", SS)
}
plot_ss(x = evals$bty_avg, y = evals$score, showSquares = TRUE)

## Click two points to make a line.
                                
## Call:
## lm(formula = y ~ x, data = pts)
## 
## Coefficients:
## (Intercept)            x  
##     3.88034      0.06664  
## 
## Sum of Squares:  131.868

Exercise 7

Normal Q-Q Plot - SIM: the residuals of the model are not normal as these values for the higher quantiles are less than what a normal distribution would predict.

m_bty_gen$fitted.values: there are some outliers; however, overall, most of the residual values are close to the fitted values.

c(1:nrow(evals)): The condition is met. Based on the sequence when this was gathered, the residuals show that they were randomly gathered.

evals\(gender & evals\)bty_avg: there is a linear relationship between gender and evaluation score. The median score and variability for males and females are similar, at least in the evaluation scores. As was shown in the former exercises, there is a linear relationship between beauty average and teaching evaluation score.

suppressWarnings(library(StMoSim))

qqnormSim(m_bty_gen$residuals)

plot(abs(m_bty_gen$residuals) ~ m_bty_gen$fitted.values)

plot(m_bty_gen$residuals ~ c(1:nrow(evals)))

plot(evals$score ~ evals$gender)

plot(evals$score ~ evals$bty_avg)

Exercise 8

Yes, this is still a significant predictor. Gender made beauty average more significant as the p-value computed is smaller now when compared to a model where beauty average was the only variable.

Exercise 9

Equation: scoreˆ=β^0+β1×bty_avg+β^2×(1)=β0+β^{1×bty_avg+β}2

Note: males tend to have a higher course evaluation socre than females for professors who get the same rating.

Exercise 10

This code should also be used: multiLines(m_bty_rank), but it is not working here.

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

names(m_bty_rank)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

Exercise 11

The variable to have the least association with the professor’s evaluation score might be the “number of professors,” which is cls_profs.

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

Exercise 12

Based on this, it is the “number of professors” which has the least assosciation to “scores.” Further, it has the highest p-value in the model.

plot(evals$score ~ evals$cls_profs)

Exercise 13

Based on my interpretation, the ethnicity p-value of about 0.11 means that it has a rather weak relationship to scores and it may be dropped as part of the model.

Exercise 14

Yes, there was a slight change in the coefficients and significance of the other explanatory variables when “cls_profs” was removed. All of the values are now slightly lower. Which means, they are more significant now to the level than they were previously.

m_back <- 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_back)

## 
## 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

Exercise 15

Linear model: scoreˆβ^{5+×class_perceval+β}6×class_credits_one+β^{7×bty_avg+β}8×picture_color_colored=β^0+β1×ethnicity_not_minority+β^{2×gender_male+β}3×language_non−englist+β^4×age+

m_back2 <- lm(score ~ ethnicity + gender + language + age + cls_perc_eval + 
    cls_credits + bty_avg + pic_color, data = evals)
summary(m_back2)

## 
## 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

Exercise 16

First code: the residuals of the model appear to be nearly normal.

Second code: the absolute values of residuals against the fitted values. Based on my interpretation, the variability of the residuals are nearly constant.

Third code: the residuals are independent.

Fourth code: each variable is linearly related to the outcome. The variables are linearly related to the score, some more than others.

qqnormSim(m_back2$residuals)

plot(abs(m_back2$residuals) ~ m_back2$fitted.values)

plot(m_back2$residuals ~ c(1:nrow(evals)))

plot(evals$score ~ evals$ethnicity)

plot(evals$score ~ evals$gender)

plot(evals$score ~ evals$language)

plot(evals$score ~ evals$age)

plot(evals$score ~ evals$cls_perc_eval)

plot(evals$score ~ evals$cls_credits)

plot(evals$score ~ evals$bty_avg)

plot(evals$score ~ evals$pic_color)

Exercise 17

No, the class courses are independent of each other so the evaluation scores from one course is independent of the other even if the course is being taught by the same instructor.

Exercise 18

The instructor is not a minority and male. Likely graduated from an institution in the USA and teaches a one credit course. Must also have a high beauty average score from the students and the instructor’s class photo should not be in color. Must also be relatively young; and a somewhat high percentage of the class must have completed the evaluation.

Exercise 19

I would say no to this. The sample size of 6 it too small. Moreover, some of the predictor variables are subjective and may vary with culture. For example, beauty may be seen as subjective; although, there are certainly objective aspects to it. So, some of the non-objective predictor variables would be difficult to determine.