# Statistics coursework sample: hypothesis tests and regression on real data
# R route for https://assignmenthelper.org/statistics-coursework-sample-r-excel/
#
# What this is: the R commands for every step on that page, runnable top to bottom. The page's
# figures were computed by the Python script published beside this file; this script was then run
# in R 4.5.3 (car 3.1-5, lmtest 0.9-40, sandwich 3.1-3) on 25 September 2026 and its output,
# statistics-coursework-student-por.R.out.txt, agrees with the Python output at every figure the
# page prints. Two things to know before you compare your numbers with ours:
#   1. t.test(G3 ~ internet) orders the groups alphabetically (no, then yes), so R prints
#      t = -3.658 and a negative confidence interval. The page reports yes minus no: t = 3.658,
#      95% CI 0.53 to 1.76. Same test, same size, opposite sign.
#   2. wilcox.test prints W = 29,582, the U statistic for the "no" group; the page's U = 45,616 is
#      the "yes" group's. They are two ways of counting the same thing, and the p-value is the same.
#
# Dataset: Student Performance, UCI Machine Learning Repository id 320, licence CC BY 4.0.
#          Cortez, P. (2008). https://doi.org/10.24432/C5TG7T
#          Download https://archive.ics.uci.edu/static/public/320/student+performance.zip and
#          unzip both layers; the file used is student-por.csv (649 rows). Put it in the working
#          directory, or change the path in read.csv below.
# Packages: car (leveneTest, vif), lmtest (bptest, resettest, coeftest, coefci), sandwich (vcovHC).
#          install.packages(c("car", "lmtest", "sandwich")) once.

library(car)
library(lmtest)
library(sandwich)

# 1. Read and describe ---------------------------------------------------------------------------
d <- read.csv("student-por.csv", sep = ";")
d$internet <- factor(d$internet)     # so leveneTest does not warn about coercing the group
str(d)
summary(d$G3)
sd(d$G3)
table(d$G3)                      # the 15 zeros are the teaching point
table(d$studytime)
table(d$failures)
table(d$internet)

# 2. Two-sample t-test: G3 by internet (Welch is the default in t.test) --------------------------
leveneTest(G3 ~ internet, data = d)          # centre = median by default (Brown-Forsythe)
t.test(G3 ~ internet, data = d)              # Welch; report t, df, p and the 95% CI of the difference
t.test(G3 ~ internet, data = d, var.equal = TRUE)   # pooled, for comparison with Excel's type 2
# Cohen's d with the pooled SD (no base-R function; three lines)
g <- split(d$G3, d$internet)
s_pooled <- sqrt(((length(g$yes) - 1) * var(g$yes) + (length(g$no) - 1) * var(g$no)) / (length(g$yes) + length(g$no) - 2))
(mean(g$yes) - mean(g$no)) / s_pooled
shapiro.test(g$yes); shapiro.test(g$no)      # both groups non-normal because of the zero grades
wilcox.test(G3 ~ internet, data = d)         # Mann-Whitney, the rank-based check (see note 2 above)

# 3. One-way ANOVA: G3 by study-time band, then Tukey HSD and eta squared ------------------------
d$studytime_f <- factor(d$studytime, levels = 1:4, labels = c("under 2 h", "2 to 5 h", "5 to 10 h", "over 10 h"))
aggregate(G3 ~ studytime_f, data = d, FUN = function(x) c(n = length(x), mean = mean(x), sd = sd(x)))
leveneTest(G3 ~ studytime_f, data = d)
fit_aov <- aov(G3 ~ studytime_f, data = d)
summary(fit_aov)                             # F(3, 645), p
TukeyHSD(fit_aov)                            # which bands differ
ss <- summary(fit_aov)[[1]][["Sum Sq"]]
ss[1] / sum(ss)                              # eta squared = SS between / SS total
oneway.test(G3 ~ studytime_f, data = d)      # Welch ANOVA, the fallback if Levene rejects
kruskal.test(G3 ~ studytime_f, data = d)     # the rank-based check

# 4. Chi-square test of independence: school x higher --------------------------------------------
tab <- table(d$school, d$higher)
tab
prop.table(tab, 1)                           # row percentages
chisq.test(tab)                              # Yates continuity correction on a 2 x 2 by default
chisq.test(tab, correct = FALSE)             # what Excel's CHISQ.TEST computes
chisq.test(tab)$expected                     # expected counts (all must be at least 5)
sqrt(chisq.test(tab, correct = FALSE)$statistic / sum(tab))   # Cramer's V for a 2 x 2, from the uncorrected statistic

# 5. Multiple linear regression, model fixed before looking at the output ------------------------
fit <- lm(G3 ~ studytime + failures + absences + higher + school + sex + Medu, data = d)
summary(fit)                                 # coefficients, R squared, adjusted R squared, F, residual SE
confint(fit)                                 # 95% confidence intervals
fit_raw <- lm(G3 ~ studytime, data = d)      # study time on its own, for the raw-versus-adjusted comparison
summary(fit_raw)$coefficients; confint(fit_raw); summary(fit_raw)$r.squared
summary(update(fit, . ~ . + internet))$coefficients   # internet added: ties the t-test to the model

# 6. Diagnostics --------------------------------------------------------------------------------
par(mfrow = c(2, 2)); plot(fit); par(mfrow = c(1, 1))   # residuals vs fitted, Q-Q, scale-location, leverage
shapiro.test(residuals(fit))
bptest(fit)                                  # Breusch-Pagan, Koenker studentised by default
resettest(fit)                               # Ramsey RESET, squared and cubed fitted values by default
vif(fit)
cooks.distance(fit)[cooks.distance(fit) > 4 / nrow(d)]
max(cooks.distance(fit))
z <- subset(d, G3 == 0)                      # the 15 zero-grade students
nrow(z); sum(z$absences == 0); table(z$school); sum(z$G2 == 0); range(z$G1)
coeftest(fit, vcov = vcovHC(fit, type = "HC3"))          # re-run 1: HC3 robust standard errors
coefci(fit, vcov = vcovHC(fit, type = "HC3"))            # and their confidence intervals
fit_nz <- update(fit, data = subset(d, G3 > 0))          # re-run 2: without the zero-grade students
summary(fit_nz)
shapiro.test(residuals(fit_nz)); bptest(fit_nz); resettest(fit_nz)   # the checks repeated on re-run 2

# 6b. Which term bends the line? ---------------------------------------------------------------
# Free each suspect from its linear coding; the nested F says whether categories fit better, and
# RESET is repeated on the re-coded model.
fit_st <- update(fit, . ~ . - studytime + factor(studytime))
anova(fit, fit_st)                           # equal steps for study time: not rejected
resettest(fit_st)                            # curvature unchanged
fit_fc <- update(fit, . ~ . - failures + factor(failures))
anova(fit, fit_fc)                           # the count coding for failures: rejected
resettest(fit_fc)                            # curvature gone
fit_sq <- update(fit, . ~ . + I(failures^2))
summary(fit_sq)$coefficients["I(failures^2)", ]; resettest(fit_sq)
aggregate(G3 ~ failures, data = d, FUN = mean)   # the whole drop comes with the first failure
summary(fit_fc); confint(fit_fc)             # re-run 3: failures as categories

# 7. Leakage: never do this in the real coursework; shown only to see why -----------------------
cor(d[, c("G1", "G2", "G3")])
fit_leak <- lm(G3 ~ studytime + failures + absences + higher + school + sex + Medu + G1 + G2, data = d)
summary(fit_leak)                            # adjusted R squared jumps; the other coefficients collapse