# ----------------------
# Setup and Data Loading
# ----------------------
# Install and load required packages
if (!require("wooldridge")) install.packages("wooldridge")
## Loading required package: wooldridge
if (!require("dplyr")) install.packages("dplyr")
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
library(wooldridge)
library(dplyr)
library(ggplot2)

# Load the VOTE1 dataset
data("vote1")

Part (i) and (ii)

——————-

# Estimate model with interaction
model1 <- lm(voteA ~ prtystrA + expendA + expendB + I(expendA*expendB), data = vote1)

# Display comprehensive results
cat("Model 1 Results:\n")
## Model 1 Results:
print(summary(model1))
## 
## Call:
## lm(formula = voteA ~ prtystrA + expendA + expendB + I(expendA * 
##     expendB), data = vote1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.9999  -8.7632  -0.1726   8.2310  29.7325 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.212e+01  4.591e+00   6.995 5.99e-11 ***
## prtystrA              3.419e-01  8.799e-02   3.886 0.000146 ***
## expendA               3.828e-02  4.960e-03   7.718 1.00e-12 ***
## expendB              -3.172e-02  4.588e-03  -6.915 9.32e-11 ***
## I(expendA * expendB) -6.629e-06  7.186e-06  -0.923 0.357584    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.13 on 168 degrees of freedom
## Multiple R-squared:  0.5708, Adjusted R-squared:  0.5606 
## F-statistic: 55.86 on 4 and 168 DF,  p-value: < 2.2e-16
# Calculate and store coefficients for later use
coef <- coef(model1)

# Print partial effects formulas
cat("\nPartial Effects:\n")
## 
## Partial Effects:
cat("∂voteA/∂expendB =", coef["expendB"], "+", coef["I(expendA * expendB)"], "* expendA\n")
## ∂voteA/∂expendB = -0.03172384 + -6.629499e-06 * expendA
cat("∂voteA/∂expendA =", coef["expendA"], "+", coef["I(expendA * expendB)"], "* expendB\n")
## ∂voteA/∂expendA = 0.03828094 + -6.629499e-06 * expendB

——————-

Part (iii)

——————-

# Calculate mean of expendA
mean_expendA <- mean(vote1$expendA)
cat("\nPart (iii):\n")
## 
## Part (iii):
cat("Average expendA:", mean_expendA, "\n")
## Average expendA: 310.611
# Fix expendA at 300
fixed_expendA <- 300

# Calculate effect of $100,000 increase in expendB when expendA = 300
effect_B <- coef["expendB"] + coef["I(expendA * expendB)"] * fixed_expendA
cat("Effect of $100,000 increase in expendB when expendA = 300:", effect_B * 100, "\n")
## Effect of $100,000 increase in expendB when expendA = 300: -3.371269

——————-

Part (iv)

——————-

# Fix expendB at 100
fixed_expendB <- 100

# Calculate effect of ΔexpendA = 100
effect_A <- coef["expendA"] + coef["I(expendA * expendB)"] * fixed_expendB
cat("\nPart (iv):\n")
## 
## Part (iv):
cat("Effect of $100,000 increase in expendA when expendB = 100:", effect_A * 100, "\n")
## Effect of $100,000 increase in expendA when expendB = 100: 3.761799

——————-

Part (v)

——————-

# Calculate shareA
vote1$shareA <- vote1$expendA / (vote1$expendA + vote1$expendB)

# Estimate new model with shareA
model2 <- lm(voteA ~ prtystrA + expendA + expendB + shareA, data = vote1)

cat("\nPart (v) - Model with shareA:\n")
## 
## Part (v) - Model with shareA:
print(summary(model2))
## 
## Call:
## lm(formula = voteA ~ prtystrA + expendA + expendB + shareA, data = vote1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5859  -3.3764  -0.3036   3.2095  31.5656 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.195374   2.567856   7.086 3.63e-11 ***
## prtystrA     0.157272   0.049685   3.165  0.00184 ** 
## expendA     -0.006670   0.002833  -2.354  0.01971 *  
## expendB      0.004267   0.002607   1.637  0.10351    
## shareA      49.439436   2.530859  19.535  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.167 on 168 degrees of freedom
## Multiple R-squared:  0.8681, Adjusted R-squared:  0.865 
## F-statistic: 276.5 on 4 and 168 DF,  p-value: < 2.2e-16

——————-

Part (vi)

——————-

# Calculate partial effect at expendA = 300 and expendB = 0
expendA_eval <- 300
expendB_eval <- 0

# Calculate shareA at these values (adding small constant to avoid division by zero)
shareA_eval <- expendA_eval / (expendA_eval + expendB_eval + 0.000001)

# Calculate partial effect
beta_expendB <- coef(model2)["expendB"]
beta_shareA <- coef(model2)["shareA"]
partial_effect <- beta_expendB + beta_shareA * (-expendA_eval/(expendA_eval + expendB_eval + 0.000001)^2)

cat("\nPart (vi):\n")
## 
## Part (vi):
cat("Partial effect of expendB when expendA = 300 and expendB = 0:", partial_effect, "\n")
## Partial effect of expendB when expendA = 300 and expendB = 0: -0.1605312