Dataset integrity. The ZIP contains the original, unmodified CSV data. A leading CSV-generated row-index column is omitted from the regression variables.

Analysis design. The project uses the exact CSV files supplied with the assignment and base-R OLS. Conventional standard errors assume homoskedastic errors. “White robust” errors mean the original HC0 Eicker–Huber–White sandwich estimator. All tests are two-sided or joint tests at the 5% significance level. Model estimates are associational, not necessarily causal.

1 Part I: Computing exercises

1.1 (i) Descriptive statistics for equation (1)

The dependent variable is course score (percentage points). The six regressors are weekly study hours, high-school economics, entering college GPA, male, prior calculus, and good past attendance.

Part I — sample descriptive statistics
Variable N Mean SD Minimum Maximum
score 599 72.110 13.097 22.130 97.540
study 599 13.995 7.843 0.000 50.000
econhs 599 0.364 0.482 0.000 1.000
colgpa 599 2.806 0.539 1.104 4.000
male 599 0.511 0.500 0.000 1.000
calculus 599 0.676 0.468 0.000 1.000
attgood 599 0.571 0.495 0.000 1.000

There are 599 observations in the supplied Part I dataset, of which 599 are available for equation (1).

1.2 (ii) OLS coefficients and usual standard errors

We estimate

\[ \mathrm{score}_i = \beta_1 + \beta_2\mathrm{study}_i + \beta_3\mathrm{econhs}_i + \beta_4\mathrm{colgpa}_i + \beta_5\mathrm{male}_i + \beta_6\mathrm{calculus}_i + \beta_7\mathrm{attgood}_i + u_i. \]

OLS with conventional (non-robust) standard errors
Variable Estimate Std_Error t_stat p_value
5 male 3.688 0.877 4.207 <0.0001
6 calculus 4.647 0.925 5.026 <0.0001

Holding other included regressors constant, the male indicator has an estimated conditional score difference of 3.688 percentage points, comparing male to non-male students. The calculus indicator has an estimated conditional score difference of 4.647 percentage points, comparing students who took calculus with those who did not. Neither coefficient alone establishes a causal effect.

1.3 (iii) Tests of the calculus and good-attendance coefficients

Each null is \(H_0:\beta_j=0\) against a two-sided alternative. Under the classical model, \(t=\hat\beta_j/\widehat{\mathrm{SE}}(\hat\beta_j)\) uses 592 residual degrees of freedom.

Conventional individual coefficient tests
Variable Estimate Std_Error Test_statistic p_value Distribution Decision_5pct
calculus 4.647 0.925 5.026 <0.0001 t(592) Reject H0
attgood 0.186 0.865 0.215 0.8296 t(592) Do not reject H0

At 5%, we reject the zero-calculus-effect null and do not reject the zero-good-attendance-effect null, conditional on the maintained model.

1.4 (iv) Wald, likelihood-ratio and Lagrange-multiplier tests

The joint null is \(H_0:\beta_2=\beta_3=\beta_7=0\) (study, economics in high school, and good attendance), with three restrictions. The restricted model keeps intercept, college GPA, male, and calculus. Both models use the same sample. The reference critical value is \(\chi^2_{3,0.95}=\) 7.815.

We compute \(W=(R\hat\beta)'[R\widehat{\mathrm{Var}}(\hat\beta)R']^{-1}(R\hat\beta)\), \(LR=n\log(SSR_r/SSR_u)\), and the classical score statistic \(LM=nR^2_{\mathrm{aux}}\) using restricted-model residuals and unrestricted regressors. These are asymptotic chi-square comparisons (the Wald statistic is \(q\) times the corresponding classical F statistic).

Three joint tests using chi-square(3)
Test Statistic DF Critical_95pct p_value Decision_5pct
Wald 0.671 3 7.815 0.8799 Do not reject H0
LR 0.679 3 7.815 0.8781 Do not reject H0
LM 0.679 3 7.815 0.8782 Do not reject H0

At the requested 5% significance level, all three tests are interpreted against the stated chi-square reference distribution. Their finite-sample values need not agree exactly.

1.5 (v) White’s heteroskedasticity test

We regress \(\hat u_i^2\) on a constant, all six original regressors, their squares, and their pairwise interactions. Squared binary regressors duplicate the binary regressors; the actual degrees of freedom equal the rank of the auxiliary design minus one, rather than counting duplicated columns. The LM statistic is \(nR^2\).

White LM heteroskedasticity test
Statistic DF Critical_95pct p_value Candidate_Auxiliary_Terms Independent_Auxiliary_Terms
56.263 23 35.172 0.0001 27.000 23.000

The test detects evidence against homoskedasticity at 5% (\(p\) = 0.0001; reference \(\chi^2_{23}\)). Failure to reject would not prove errors are homoskedastic.

1.6 (vi) White HC0 standard errors and repeated individual tests

HC0 estimates \(\widehat V_{HC0}=(X'X)^{-1}X'\operatorname{diag}(\hat u_i^2)X(X'X)^{-1}\). The reported robust two-sided \(p\)-values use the large-sample standard-normal reference, in contrast with the classical Student \(t\) distribution above.

HC0 robust tests of calculus and good attendance
Variable Estimate Std_Error Test_statistic p_value Distribution Decision_5pct
calculus 4.647 0.943 4.928 <0.0001 Asymptotic N(0, 1) Reject H0
attgood 0.186 0.866 0.215 0.8297 Asymptotic N(0, 1) Do not reject H0

Using HC0 errors, the calculus null is rejected and the good-attendance null is not rejected at 5%. Changes relative to (iii) arise from the estimated covariance matrix and its large-sample reference.

1.7 (vii) Expanded regression and t-ratio rankings

This specification adds all other genuine variables in the supplied econmath file as regressors. The first CSV field is a row index, not a substantive variable. The ACT variables have missing values, so listwise deletion leaves 572 observations rather than the original 599; the expanded and baseline models should not be compared without acknowledging that change.

Here, the assignment’s “magnitude of t-ratio” is interpreted as absolute \(|\hat\beta/SE|\) (both positive and negative predictive associations can have large magnitude). HC0 uses the analogous absolute robust z-ratio.

All expanded-model regressors ranked by absolute ratio
Variable t_Usual Rank_Usual z_HC0 Rank_HC0
colgpa 11.641 1 11.225 1
calculus 3.677 2 3.496 2
male 3.064 3 3.051 3
actmth 2.547 4 2.508 4
work -2.372 5 -2.293 5
hsgpa 2.092 6 2.175 6
act 1.823 7 1.686 7
acteng -0.812 8 -0.794 8
age -0.779 9 -0.614 10
attexc 0.675 10 0.637 9
mothcoll -0.570 11 -0.580 11
attgood 0.468 12 0.433 12
fathcoll -0.239 13 -0.254 13
econhs -0.208 14 -0.213 14
study -0.095 15 -0.101 15
mathscr -0.030 16 -0.028 16

The top three variables using conventional errors are colgpa, calculus, male. Using HC0 errors they are colgpa, calculus, male. The top-three ordering is unchanged. A large absolute statistic measures evidence relative to sampling uncertainty, not economic importance, causal importance, or practical effect size. Comparisons can be influenced by measurement units, variable variation, and correlation across regressors.

1.8 (viii) Limitations of identifying the effect of calculus

The calculus indicator may be endogenous: students self-select into calculus based on preparation, ability, aspirations or advising, so omitted determinants may affect both course choice and score. The expanded model controls for age, work, high-school GPA, test scores, prior math-quiz performance, excellent attendance, and parental college education, but control inclusion does not itself create random assignment. Some factors may be affected by prior coursework or share measurement content with calculus; timing and a defensible causal design would be necessary before interpreting coefficients causally.

There is also multicollinearity between some academic measures (ACT components and composite), which can inflate standard errors. In the expanded sample, the auxiliary variance inflation factor for calculus is 1.19. Heteroskedasticity changes inference and motivates HC0 standard errors, but robust errors do not correct omitted-variable bias, selection bias, or measurement error. Finally, adding the ACT variables changes the complete-case sample; both specification and sample composition can influence the estimated calculus coefficient.

1.9 Optional regression diagnostics (generated from these data)

2 Part II: Alcohol abuse and employment

Empirical question. Is alcohol abuse associated with the probability of employment after accounting for individual socioeconomic, health and location characteristics? We estimate an OLS linear probability model (LPM), where employ = 1 means employed and abuse = 1 denotes alcohol abuse. The regression includes age and its square, education, marriage, family size, state unemployment rate, race indicator, health dummies, region dummies and metropolitan-location dummies. The omitted health, regional and location categories serve as references. Interview quarter indicators are omitted because they are constant across this subsample.

Sample profile (rates are proportions)
Sample_N Alcohol_abuse_rate Employment_rate In_workforce_rate
387.000 0.114 0.891 0.917
Employment and workforce participation by abuse status
abuse N employ inwf
0.000 343 0.895 0.921
1.000 44 0.864 0.886
Adjusted LPM coefficient for alcohol abuse
Term Coefficient HC0_SE p_value CI95_Lower CI95_Upper
Alcohol abuse indicator -0.024 0.053 0.6447 -0.128 0.079

The conditional association between alcohol abuse and employment is -2.43 percentage points (HC0 standard error 5.27 percentage points; p = 0.6447). At the 5% level we do not reject the hypothesis of a zero conditional association. The main model uses 387 observations and has \(R^2 =\) 0.209. A studentized residual-squared auxiliary diagnostic yields \(nR^2=\) 42.830 (\(df =\) 17, \(p =\) 0.0005). Hence heteroskedasticity-robust standard errors are reported. The small number of non-employed observations also limits precision.

Sensitivity analyses: HC0 inference for the abuse indicator
Outcome N Abuse_Estimate HC0_SE p_value
Employed in full sample 387 -0.024 0.053 0.6447
In the labor force 387 -0.034 0.048 0.4883
Employed if in labor force 355 0.009 0.030 0.7616

Econometric interpretation. The two sensitivity regressions change the outcome to labor-force participation and employment among labor-force participants. The latter changes the sample and conditions on participation, so it is not directly the same estimand as overall employment. The LPM may generate predictions outside [0,1] and generally has heteroskedastic disturbances. Most importantly, cross-sectional OLS does not resolve reverse causality, unobserved confounding, or measurement error in alcohol abuse: none of these coefficients should be described as the causal effect of abuse on employment. This same report includes the concise Part II results, with detailed results in the appendix.

3 Appendix A: Complete R analysis code

All executable code is embedded directly in this .Rmd file. The core analysis source is reproduced below for inspection, so no external R scripts are required.

# ECMT 6002/6702: Reproducible groupwork econometric computations.
# All computations use base R; report rendering uses rmarkdown + knitr.
# All calculations are embedded inside this single .Rmd file.

options(stringsAsFactors = FALSE, scipen = 8)

read_assignment_csv <- function(path) {
  if (!file.exists(path)) stop("Missing data file: ", path)
  d <- utils::read.csv(path, check.names = TRUE)
  # CSV export contains an observation-index column, not an explanatory variable.
  index_names <- intersect(names(d), c("Unnamed..0", "X", "X.1"))
  if (length(index_names) > 0L) {
    for (v in index_names) {
      if (is.numeric(d[[v]]) && identical(as.integer(d[[v]]), seq_len(nrow(d)))) {
        d[[v]] <- NULL
      }
    }
  }
  d
}

assert_variables <- function(data, variables, label) {
  missing <- setdiff(variables, names(data))
  if (length(missing)) stop(label, " is missing: ", paste(missing, collapse = ", "))
  if (any(!vapply(data[variables], is.numeric, logical(1)))) {
    stop(label, " requires numeric variables in: ", paste(variables, collapse = ", "))
  }
}

descriptive <- function(data, variables) {
  out <- lapply(variables, function(nm) {
    v <- data[[nm]]
    v <- v[is.finite(v)]
    data.frame(Variable = nm, N = length(v), Mean = mean(v),
               SD = stats::sd(v), Minimum = min(v), Maximum = max(v))
  })
  do.call(rbind, out)
}

usual_coefs <- function(model) {
  ctab <- stats::coef(summary(model))
  if (anyNA(stats::coef(model))) stop("Aliased coefficient: revise model specification")
  data.frame(Variable = rownames(ctab), Estimate = unname(ctab[, 1]),
             Std_Error = unname(ctab[, 2]), t_stat = unname(ctab[, 3]),
             p_value = unname(ctab[, 4]), row.names = NULL)
}

white_vcov_hc0 <- function(model) {
  # Eicker-Huber-White / White (HC0) covariance: (X'X)^-1 X'diag(e^2)X (X'X)^-1.
  x <- stats::model.matrix(model)
  if (qr(x)$rank < ncol(x)) stop("Cannot compute HC0: regression matrix is singular")
  e <- stats::residuals(model)
  bread <- solve(crossprod(x))
  meat <- crossprod(x, sweep(x, 1, as.numeric(e^2), `*`))
  ans <- bread %*% meat %*% bread
  dimnames(ans) <- list(colnames(x), colnames(x))
  ans
}

white_coefs <- function(model) {
  beta <- stats::coef(model)
  ses <- sqrt(diag(white_vcov_hc0(model)))
  stat <- beta / ses
  data.frame(Variable = names(beta), Estimate = unname(beta),
             Std_Error = unname(ses), z_stat = unname(stat),
             p_value = 2 * stats::pnorm(abs(stat), lower.tail = FALSE),
             CI_Lower = unname(beta - stats::qnorm(.975) * ses),
             CI_Upper = unname(beta + stats::qnorm(.975) * ses),
             row.names = NULL)
}

simple_tests <- function(model, variables, robust = FALSE) {
  coef_table <- if (robust) white_coefs(model) else usual_coefs(model)
  rows <- match(variables, coef_table$Variable)
  if (anyNA(rows)) stop("Missing requested term(s): ", paste(variables, collapse = ", "))
  tab <- coef_table[rows, , drop = FALSE]
  stat_column <- if (robust) "z_stat" else "t_stat"
  data.frame(Variable = tab$Variable, Estimate = tab$Estimate,
             Std_Error = tab$Std_Error, Test_statistic = tab[[stat_column]],
             p_value = tab$p_value,
             Distribution = if (robust) "Asymptotic N(0, 1)" else
               paste0("t(", stats::df.residual(model), ")"),
             Decision_5pct = ifelse(tab$p_value < .05, "Reject H0", "Do not reject H0"),
             row.names = NULL)
}

joint_wald_lr_lm <- function(full_model, restricted_model, hypothesis_terms) {
  # Both models MUST use identical observations and classical homoskedastic OLS.
  if (stats::nobs(full_model) != stats::nobs(restricted_model)) {
    stop("Restricted and unrestricted models use different samples")
  }
  q <- length(hypothesis_terms)
  b <- stats::coef(full_model)
  selector <- matrix(0, q, length(b),
                     dimnames = list(hypothesis_terms, names(b)))
  for (j in seq_along(hypothesis_terms)) selector[j, hypothesis_terms[j]] <- 1
  rbeta <- as.vector(selector %*% b)
  rvrt <- selector %*% stats::vcov(full_model) %*% t(selector)
  wald <- as.numeric(crossprod(rbeta, solve(rvrt, rbeta)))
  n <- stats::nobs(full_model)
  rss_f <- sum(stats::residuals(full_model)^2)
  rss_r <- sum(stats::residuals(restricted_model)^2)
  lr <- n * log(rss_r / rss_f)
  # Classical score/LM: regress restricted-model residuals on all unrestricted X.
  e_restricted <- stats::residuals(restricted_model)
  aux <- stats::lm.fit(stats::model.matrix(full_model), e_restricted)
  lm_stat <- n * (1 - sum(aux$residuals^2) / sum(e_restricted^2))
  values <- c(Wald = wald, LR = lr, LM = lm_stat)
  out <- data.frame(Test = names(values), Statistic = unname(values),
                    DF = q,
                    Critical_95pct = stats::qchisq(.95, df = q),
                    p_value = stats::pchisq(values, df = q, lower.tail = FALSE),
                    row.names = NULL)
  out$Decision_5pct <- ifelse(out$p_value < .05, "Reject H0", "Do not reject H0")
  out
}

white_heteroskedasticity <- function(model) {
  # White auxiliary regression with all original regressors, squares and products.
  # Redundant squares of binary variables are automatically excluded via QR rank.
  x <- stats::model.matrix(model)[, -1, drop = FALSE]
  z <- x
  for (j in seq_len(ncol(x))) {
    new_col <- matrix(x[, j]^2, ncol = 1,
                      dimnames = list(NULL, paste0("sq_", colnames(x)[j])))
    z <- cbind(z, new_col)
  }
  if (ncol(x) >= 2L) {
    for (j in seq_len(ncol(x) - 1L)) {
      for (k in (j + 1L):ncol(x)) {
        new_col <- matrix(x[, j] * x[, k], ncol = 1,
                          dimnames = list(NULL,
                            paste0("cross_", colnames(x)[j], "_", colnames(x)[k])))
        z <- cbind(z, new_col)
      }
    }
  }
  e2 <- as.numeric(stats::residuals(model)^2)
  aux <- stats::lm.fit(cbind("(Intercept)" = 1, z), e2)
  tss <- sum((e2 - mean(e2))^2)
  r2 <- 1 - sum(aux$residuals^2) / tss
  r2 <- max(0, min(1, r2))
  df <- aux$rank - 1L
  if (df <= 0) stop("White auxiliary regression has no identifiable predictors")
  lm <- length(e2) * r2
  list(statistic = lm, df = df, critical = stats::qchisq(.95, df),
       p_value = stats::pchisq(lm, df, lower.tail = FALSE),
       candidate_terms = ncol(z), independent_terms = df,
       n = length(e2), r_squared = r2)
}

auxiliary_heteroskedasticity <- function(model) {
  # Koenker-style studentized auxiliary n*R2 diagnostic on original regressors.
  e2 <- as.numeric(stats::residuals(model)^2)
  z <- stats::model.matrix(model)
  fit <- stats::lm.fit(z, e2)
  r2 <- max(0, 1 - sum(fit$residuals^2) / sum((e2 - mean(e2))^2))
  df <- fit$rank - 1L
  stat <- length(e2) * r2
  data.frame(Test = "Studentized auxiliary nR2 diagnostic",
             Statistic = stat, DF = df,
             p_value = stats::pchisq(stat, df, lower.tail = FALSE),
             Decision_5pct = ifelse(stats::pchisq(stat, df, lower.tail = FALSE) < .05,
                                    "Reject homoskedasticity", "Do not reject homoskedasticity"))
}

rank_by_t <- function(model) {
  classic <- usual_coefs(model)
  robust <- white_coefs(model)
  classic <- classic[classic$Variable != "(Intercept)", , drop = FALSE]
  robust <- robust[robust$Variable != "(Intercept)", , drop = FALSE]
  matched <- match(classic$Variable, robust$Variable)
  tab <- data.frame(Variable = classic$Variable,
                    t_Usual = classic$t_stat,
                    Rank_Usual = rank(-abs(classic$t_stat), ties.method = "min"),
                    z_HC0 = robust$z_stat[matched],
                    Rank_HC0 = rank(-abs(robust$z_stat[matched]), ties.method = "min"),
                    row.names = NULL)
  tab <- tab[order(tab$Rank_Usual), , drop = FALSE]
  rownames(tab) <- NULL
  tab
}

save_table <- function(data, filename) {
  utils::write.csv(data, file.path(tables_dir, filename), row.names = FALSE,
                   na = "", fileEncoding = "UTF-8")
}

format_p <- function(x) {
  ifelse(is.na(x), "NA", ifelse(x < .0001, "<0.0001", sprintf("%.4f", x)))
}

# --- INPUTS ------------------------------------------------------------------
# When using RStudio Knit, working directory is this .Rmd folder.
# When rendering from R Console, use the folder containing the Rmd.
project_root <- getwd()
tables_dir <- file.path(project_root, "output", "tables")
figures_dir <- file.path(project_root, "output", "figures")
images_dir <- file.path(project_root, "output", "result_images")
for (path in c(tables_dir, figures_dir, images_dir)) {
  dir.create(path, recursive = TRUE, showWarnings = FALSE)
}
econ <- read_assignment_csv("econmath_subsample.csv")
alc <- read_assignment_csv("alcohol_subsample.csv")
assert_variables(econ,
  c("score", "study", "econhs", "colgpa", "male", "calculus", "attgood"),
  "Economics dataset")
assert_variables(alc, c("abuse", "status", "employ", "inwf", "age", "agesq",
                      "educ", "married", "famsize", "unemrate"), "Alcohol dataset")

# --- PART I: ORIGINAL SPECIFICATION -----------------------------------------
p1_variables <- c("score", "study", "econhs", "colgpa", "male", "calculus", "attgood")
p1_desc <- descriptive(econ, p1_variables)
p1_formula <- score ~ study + econhs + colgpa + male + calculus + attgood
p1_ols <- stats::lm(p1_formula, data = econ, na.action = stats::na.omit)
p1_ols_all <- usual_coefs(p1_ols)
p1_ols_asked <- p1_ols_all[p1_ols_all$Variable %in% c("male", "calculus"), , drop = FALSE]
p1_individual <- simple_tests(p1_ols, c("calculus", "attgood"), robust = FALSE)

# Joint H0: coefficients of study, econhs, attgood all equal zero.
p1_restricted <- stats::lm(score ~ colgpa + male + calculus,
                           data = stats::model.frame(p1_ols))
p1_joint <- joint_wald_lr_lm(p1_ols, p1_restricted,
                             c("study", "econhs", "attgood"))
p1_white <- white_heteroskedasticity(p1_ols)
p1_hc0_all <- white_coefs(p1_ols)
p1_hc0_asked <- p1_hc0_all[p1_hc0_all$Variable %in% c("calculus", "attgood"),,
                             drop = FALSE]
p1_robust_tests <- simple_tests(p1_ols, c("calculus", "attgood"), robust = TRUE)

# --- PART I: EXPANDED SPECIFICATION -----------------------------------------
# Include EVERY other genuine explanatory variable, not the CSV row-number index.
p1_extended_predictors <- setdiff(names(econ), "score")
p1_extended_formula <- stats::reformulate(p1_extended_predictors, response = "score")
p1_extended <- stats::lm(p1_extended_formula, data = econ,
                           na.action = stats::na.omit)
p1_extended_usual <- usual_coefs(p1_extended)
p1_extended_hc0 <- white_coefs(p1_extended)
p1_ranking <- rank_by_t(p1_extended)
p1_top_usual <- head(p1_ranking$Variable, 3)
p1_top_hc0 <- head(p1_ranking$Variable[order(p1_ranking$Rank_HC0)], 3)
p1_expanded_sample <- stats::nobs(p1_extended)

# Auxiliary VIF for calculus (descriptive diagnostic, not causal identification).
p1_extended_frame <- stats::model.frame(p1_extended)
p1_other_predictors <- setdiff(p1_extended_predictors, "calculus")
p1_vif_calculus_aux <- stats::lm(
  stats::reformulate(p1_other_predictors, response = "calculus"),
  data = p1_extended_frame)
p1_calculus_vif <- 1 / (1 - summary(p1_vif_calculus_aux)$r.squared)

# --- PART II: ALCOHOL ABUSE AND EMPLOYMENT -----------------------------------
# Interview quarter has no within-sample variation in this particular subsample:
# qrt1 == 1, qrt2 == 0, qrt3 == 0 for everyone. Omitting them avoids singularity.
p2_constants <- names(alc)[vapply(alc, function(x) {
  vals <- unique(stats::na.omit(x))
  length(vals) <= 1L
}, logical(1))]
p2_valid_status <- all(alc$employ == as.integer(alc$status == 3), na.rm = TRUE) &&
  all(alc$inwf == as.integer(alc$status > 1), na.rm = TRUE)
if (!p2_valid_status) stop("employ/inwf indicators disagree with status coding")
p2_variables <- c("abuse", "employ", "inwf", "age", "educ", "unemrate")
p2_desc <- descriptive(alc, p2_variables)
p2_by_abuse <- aggregate(cbind(employ, inwf) ~ abuse, data = alc,
                         FUN = function(v) mean(v, na.rm = TRUE))
p2_by_abuse$N <- as.integer(table(alc$abuse)[as.character(p2_by_abuse$abuse)])
p2_by_abuse <- p2_by_abuse[, c("abuse", "N", "employ", "inwf")]

p2_controls <- c("age", "agesq", "educ", "married", "famsize", "unemrate",
                 "white", "exhealth", "vghealth", "goodhealth", "fairhealth",
                 "northeast", "midwest", "south", "centcity", "outercity")
p2_predictors <- c("abuse", p2_controls)
assert_variables(alc, c("employ", "inwf", p2_predictors), "Part II")
p2_formula <- stats::reformulate(p2_predictors, response = "employ")
p2_ols <- stats::lm(p2_formula, data = alc, na.action = stats::na.omit)
p2_usual <- usual_coefs(p2_ols)
p2_hc0 <- white_coefs(p2_ols)
p2_abuse_robust <- p2_hc0[p2_hc0$Variable == "abuse", , drop = FALSE]
p2_aux_test <- auxiliary_heteroskedasticity(p2_ols)

# Sensitivity: model labor-force participation, and employment conditional on participation.
p2_labor_model <- stats::lm(stats::reformulate(p2_predictors, response = "inwf"),
                            data = alc, na.action = stats::na.omit)
p2_labor_hc0 <- white_coefs(p2_labor_model)
p2_workforce <- subset(alc, inwf == 1)
p2_conditional_model <- stats::lm(p2_formula, data = p2_workforce,
                                  na.action = stats::na.omit)
p2_conditional_hc0 <- white_coefs(p2_conditional_model)
p2_robustness <- data.frame(
  Outcome = c("Employed in full sample", "In the labor force", "Employed if in labor force"),
  N = c(stats::nobs(p2_ols), stats::nobs(p2_labor_model),
        stats::nobs(p2_conditional_model)),
  Abuse_Estimate = c(p2_abuse_robust$Estimate,
                     p2_labor_hc0$Estimate[p2_labor_hc0$Variable == "abuse"],
                     p2_conditional_hc0$Estimate[p2_conditional_hc0$Variable == "abuse"]),
  HC0_SE = c(p2_abuse_robust$Std_Error,
             p2_labor_hc0$Std_Error[p2_labor_hc0$Variable == "abuse"],
             p2_conditional_hc0$Std_Error[p2_conditional_hc0$Variable == "abuse"]),
  p_value = c(p2_abuse_robust$p_value,
              p2_labor_hc0$p_value[p2_labor_hc0$Variable == "abuse"],
              p2_conditional_hc0$p_value[p2_conditional_hc0$Variable == "abuse"]))

# Export all numerical tables so they remain usable outside the HTML report.
save_table(p1_desc, "PartI_Descriptive.csv")
save_table(p1_ols_all, "PartI_OLS_Conventional.csv")
save_table(p1_individual, "PartI_Individual_Tests.csv")
save_table(p1_joint, "PartI_Wald_LR_LM.csv")
save_table(p1_hc0_all, "PartI_OLS_White_HC0.csv")
save_table(p1_robust_tests, "PartI_Robust_Tests.csv")
save_table(p1_extended_usual, "PartI_Expanded_OLS.csv")
save_table(p1_extended_hc0, "PartI_Expanded_White_HC0.csv")
save_table(p1_ranking, "PartI_Rankings.csv")
save_table(p2_desc, "PartII_Descriptive.csv")
save_table(p2_by_abuse, "PartII_By_Abuse.csv")
save_table(p2_usual, "PartII_OLS_Conventional.csv")
save_table(p2_hc0, "PartII_OLS_White_HC0.csv")
save_table(p2_robustness, "PartII_Robustness.csv")

# Output quick-check text files useful for actual RStudio screenshot capture.
utils::capture.output({
  cat("PART I: OLS ESTIMATION\n\n")
  print(summary(p1_ols))
  cat("\nPART I: JOINT CHI-SQUARE TESTS\n")
  print(p1_joint, row.names = FALSE)
  cat("\nPART I: WHITE TEST\n")
  print(p1_white)
  cat("\nPART I: EXPANDED MODEL t-RANKING\n")
  print(p1_ranking, row.names = FALSE)
}, file = file.path(project_root, "output", "PartI_Console_Results.txt"))
utils::capture.output({
  cat("PART II: EMPLOYMENT REGRESSION\n\n")
  print(summary(p2_ols))
  cat("\nPART II: WHITE HC0 INFERENCE\n")
  print(p2_hc0, row.names = FALSE)
  cat("\nPART II: ROBUSTNESS\n")
  print(p2_robustness, row.names = FALSE)
}, file = file.path(project_root, "output", "PartII_Console_Results.txt"))

4 Appendix B: Complete regression outputs and computed results

4.1 Part I: All classical OLS estimates

Part I full conventional coefficient table
Variable Estimate Std_Error t_stat p_value
(Intercept) 28.583 2.642 10.817 <0.0001
study -0.033 0.055 -0.603 0.5468
econhs -0.418 0.877 -0.477 0.6338
colgpa 13.902 0.803 17.319 <0.0001
male 3.688 0.877 4.207 <0.0001
calculus 4.647 0.925 5.026 <0.0001
attgood 0.186 0.865 0.215 0.8296

4.2 Part I: All HC0 robust estimates

Part I full HC0 coefficient table
Variable Estimate Std_Error z_stat p_value CI_Lower CI_Upper
(Intercept) 28.583 2.893 9.878 <0.0001 22.912 34.254
study -0.033 0.058 -0.565 0.5718 -0.147 0.081
econhs -0.418 0.858 -0.487 0.6263 -2.101 1.264
colgpa 13.902 0.824 16.870 <0.0001 12.287 15.517
male 3.688 0.877 4.205 <0.0001 1.969 5.407
calculus 4.647 0.943 4.928 <0.0001 2.799 6.496
attgood 0.186 0.866 0.215 0.8297 -1.512 1.884

4.3 Part I: Expanded-model estimates

Expanded model — conventional inference
Variable Estimate Std_Error t_stat p_value
(Intercept) 15.340 10.154 1.511 0.1314
age -0.360 0.462 -0.779 0.4361
work -0.109 0.046 -2.372 0.0181
study -0.005 0.055 -0.095 0.9241
econhs -0.182 0.875 -0.208 0.8355
colgpa 10.957 0.941 11.641 <0.0001
hsgpa 3.165 1.513 2.092 0.0369
acteng -0.175 0.215 -0.812 0.4169
actmth 0.465 0.183 2.547 0.0111
act 0.560 0.307 1.823 0.0688
mathscr -0.008 0.272 -0.030 0.9759
male 2.783 0.908 3.064 0.0023
calculus 3.511 0.955 3.677 0.0003
attexc 1.049 1.553 0.675 0.4999
attgood 0.649 1.388 0.468 0.6400
fathcoll -0.209 0.873 -0.239 0.8110
mothcoll -0.530 0.930 -0.570 0.5686
Expanded model — HC0 inference
Variable Estimate Std_Error z_stat p_value CI_Lower CI_Upper
(Intercept) 15.340 12.174 1.260 0.2076 -8.520 39.200
age -0.360 0.587 -0.614 0.5393 -1.510 0.790
work -0.109 0.048 -2.293 0.0218 -0.202 -0.016
study -0.005 0.052 -0.101 0.9199 -0.108 0.097
econhs -0.182 0.852 -0.213 0.8310 -1.852 1.489
colgpa 10.957 0.976 11.225 <0.0001 9.044 12.870
hsgpa 3.165 1.455 2.175 0.0296 0.313 6.016
acteng -0.175 0.220 -0.794 0.4270 -0.606 0.256
actmth 0.465 0.185 2.508 0.0122 0.102 0.829
act 0.560 0.332 1.686 0.0918 -0.091 1.211
mathscr -0.008 0.293 -0.028 0.9776 -0.583 0.566
male 2.783 0.912 3.051 0.0023 0.995 4.570
calculus 3.511 1.005 3.496 0.0005 1.542 5.480
attexc 1.049 1.645 0.637 0.5239 -2.176 4.274
attgood 0.649 1.500 0.433 0.6651 -2.291 3.590
fathcoll -0.209 0.821 -0.254 0.7994 -1.819 1.401
mothcoll -0.530 0.914 -0.580 0.5618 -2.322 1.261

4.4 Part II: All estimated coefficients and diagnostic

Part II conventional estimates
Variable Estimate Std_Error t_stat p_value
(Intercept) -0.177 0.309 -0.573 0.5670
abuse -0.024 0.046 -0.526 0.5994
age 0.017 0.015 1.161 0.2463
agesq -0.000 0.000 -1.203 0.2296
educ 0.015 0.005 2.843 0.0047
married 0.009 0.043 0.214 0.8305
famsize 0.017 0.011 1.531 0.1267
unemrate -0.001 0.014 -0.041 0.9671
white 0.031 0.046 0.671 0.5028
exhealth 0.471 0.083 5.685 <0.0001
vghealth 0.494 0.085 5.808 <0.0001
goodhealth 0.460 0.086 5.339 <0.0001
fairhealth 0.328 0.098 3.344 0.0009
northeast 0.100 0.056 1.785 0.0752
midwest 0.038 0.040 0.949 0.3435
south 0.059 0.042 1.411 0.1590
centcity -0.047 0.042 -1.116 0.2653
outercity 0.005 0.042 0.119 0.9051
Part II HC0 estimates
Variable Estimate Std_Error z_stat p_value CI_Lower CI_Upper
(Intercept) -0.177 0.348 -0.509 0.6109 -0.858 0.505
abuse -0.024 0.053 -0.461 0.6447 -0.128 0.079
age 0.017 0.016 1.026 0.3047 -0.015 0.049
agesq -0.000 0.000 -1.041 0.2977 -0.001 0.000
educ 0.015 0.006 2.676 0.0075 0.004 0.027
married 0.009 0.043 0.214 0.8306 -0.075 0.094
famsize 0.017 0.009 1.811 0.0702 -0.001 0.035
unemrate -0.001 0.011 -0.055 0.9561 -0.021 0.020
white 0.031 0.051 0.598 0.5496 -0.070 0.132
exhealth 0.471 0.126 3.729 0.0002 0.223 0.718
vghealth 0.494 0.127 3.882 0.0001 0.245 0.743
goodhealth 0.460 0.128 3.582 0.0003 0.208 0.712
fairhealth 0.328 0.153 2.146 0.0319 0.028 0.627
northeast 0.100 0.049 2.042 0.0412 0.004 0.196
midwest 0.038 0.045 0.846 0.3975 -0.050 0.125
south 0.059 0.046 1.294 0.1958 -0.030 0.148
centcity -0.047 0.051 -0.932 0.3514 -0.147 0.052
outercity 0.005 0.041 0.123 0.9023 -0.075 0.085
Part II studentized auxiliary heteroskedasticity diagnostic
Test Statistic DF p_value Decision_5pct
Studentized auxiliary nR2 diagnostic 42.830 17 0.0005 Reject homoskedasticity

4.5 Sample consistency, missingness, and replication

The status variable exactly agrees with derived employ (status = 3) and inwf (status > 1) indicators: Yes. Constant Part II variables: qrt1, qrt2, qrt3. Variables with missing Part I values: acteng, actmth, act. The report writes every numerical table to output/tables/ and full console-style results to output/PartI_Console_Results.txt and output/PartII_Console_Results.txt.

Screenshots: The brief requests actual screenshots for Part I. Open the output text files in RStudio or view the RStudio Console, capture real screenshots, and attach them to the submission if your instructor requires literal application screenshots. The report’s executable code and reproduced tables are not falsely presented as screenshots.

## **R version:** 4.5.3; **Platform:** aarch64-apple-darwin20; **Execution date:** 2026-10-23 13:48:56 CST.