library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.3.2## Warning: package 'ggplot2' was built under R version 4.3.2## Warning: package 'tidyr' was built under R version 4.3.2## Warning: package 'readr' was built under R version 4.3.2## Warning: package 'purrr' was built under R version 4.3.2## Warning: package 'dplyr' was built under R version 4.3.2## Warning: package 'stringr' was built under R version 4.3.2## Warning: package 'forcats' was built under R version 4.3.2## Warning: package 'lubridate' was built under R version 4.3.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── 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(purrr)
library(plotrix)## Warning: package 'plotrix' was built under R version 4.3.2library(readxl)
library(ggthemes)## Warning: package 'ggthemes' was built under R version 4.3.2library(ggrepel)## Warning: package 'ggrepel' was built under R version 4.3.2library(effsize)## Warning: package 'effsize' was built under R version 4.3.3library(pwrss)## Warning: package 'pwrss' was built under R version 4.3.3##
## Attaching package: 'pwrss'
##
## The following object is masked from 'package:stats':
##
## power.t.testThis week’s data dive asks to build a generalized linear model (GLM) for our data using one column that can be easily converted to a binary column, along with 1-4 explanatory variables. From here, we are asked to build a logistic regression model and interpret the coefficients, along with building a confidence interval for this as well.
For my data, this can be a bit challenging because there are few columns that could reasonably be converted into binary variables. I could choose semester as the column and go with fall/spring, or I could choose previous programming experience. Ultimately, we are interested in measuring the growth in computational attitudes/skills, so I would argue previous programming experience would be the best to convert to a binary for conveying relevant information. The response variable will be the initial computational skill, and while I could do this for all 11 of those items, in the interest of time I will only look at the initial state of students “Testing their code,” which is the column labeled “Initial7.” If they have previous experience, the value will be attributed to a “1”, and if not, that value will be a “0”.
So let’s go ahead and import the dataset:
data0 <- read_excel("~/IUPUI/By Semester/Spring 2024/R_Stats/F18-F23_Survey_Data_Clean_R.xlsx")
View(data0)In this dataset, the binary column will be column 5, which is labeled “previous experience.” Right now, the column responses are y/n for “yes” and “no,” respectively. To convert those into a binary column, we can add the following column:
data1 <-
data0 |>
mutate(Experience = case_when(
startsWith(`Previous Experience`, 'y') ~ 1,
startsWith(`Previous Experience`, 'n') ~ 0))
view(data1)Now to build the GLM:
m1 <-
glm(Experience ~ Initial7,
data = data1,
family = binomial(link = "logit"))
summary(m1)##
## Call:
## glm(formula = Experience ~ Initial7, family = binomial(link = "logit"),
## data = data1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.27945 0.21267 -6.016 1.78e-09 ***
## Initial7 0.07861 0.03662 2.147 0.0318 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 639.90 on 528 degrees of freedom
## Residual deviance: 635.24 on 527 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 639.24
##
## Number of Fisher Scoring iterations: 4The column “Estimate” in the model corresponds to the values of \(\hat{\beta_0}\) and \(\hat{\beta_1}\), respectively. Because the value of \(\hat{\beta_1}\) is greater than 0, this does mean that the probability of having previous programming experience increases as the value of the students’ initial state increases. The p-value associated with \(\beta_1\) is less than 0.05, which gives a good indication of the confidence for this coefficient, further providing evidence that interpretation of \(\beta_1\) is likely “correct.”
The other coefficient, \(\beta_0\) is the y-intercept. The p-value is below 0.05, and in fact it’s 7 orders of magnitude below 0.05, so it is safe to say that this value is non-zero. The value itself has some meaning too, as it has the ability to give us some information about what happens when the x-value, in this case the present skill, is zero. In most cases, this is a fairly straightforward interpretation, but in this case, when 0 is not a valid number (because the scale starts at 1), it’s a little more challenging. However, qualitatively, \(\beta_0\) can be interpreted as the odds that a student has previous programming experience based on an initial competency level of 0. As described in the notes for this week, this coefficient can still be useful when determining 50/50 odds. Using both coefficients, we can find the 50/50 odds by:
B_0 <- 1.27945
B_1 <- 0.07861
odds <- B_0/B_1
odds## [1] 16.27592It is unsurprising that this number is somewhat nonsensical. The 50/50 odds would occur when the initial student response is ~16. This cannot exist based on the scale, which highlights the limitations of applying this procedure to this type of data. It’s also worth mentioning that the \(\beta_1\) value is just marginally above 0, so the positive correlation is weak at best. This result is somewhat surprising, and worth examining further in the future for more of the other survey items (maybe in the week 11 data dive for example).
To build the confidence intervals for the coefficient \(\beta_1\), we just need to multiply the standard error given in the glm by ~1.96. This would translate to:
upper_CI <- B_1 + .03662*1.96
lower_CI <- B_1 - 0.03662*1.96
upper_CI## [1] 0.1503852lower_CI## [1] 0.0068348Essentially, this confidence interval just tells us that the “true” value of \(\beta_1\) likely falls between ~0.0068 and ~0.1504.