R/Module 3 Assignment.R at main · OBOptics/R · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
##Load dataset
df<- bua6315_mod3_quickfix_data
View(df)

##Reformat Variabes
names(df) <- make.names(names(df))
View(df)


######QUESTION 1#########
# Load necessary packages
install.packages(c("stargazer"))
library(stargazer)

# Label models
model1 <- lm(Vehicles.Served ~ Garage.Bays + Population, data = df)
model2 <- lm(Vehicles.Served ~ Garage.Bays + Population + Access, data = df)
model3 <- lm(Vehicles.Served ~ Garage.Bays + Population + Access + Winter, data = df)


##Create Regression Model
stargazer(
  model1, model2, model3,
  type = "text",                # Use "html" or "latex" as needed
  title = "Table 1: Regression Results",
  dep.var.labels = "Vehicles.Served",
  covariate.labels = c("Garage.Bays", "Population", "Access", "Winter"),
  omit.stat = c("ser", "f"),    # We'll add them manually
  add.lines = list(
    c("Std. Error of Estimate", round(sqrt(deviance(model1)/df.residual(model1)), 3),
      round(sqrt(deviance(model2)/df.residual(model2)), 3),
      round(sqrt(deviance(model3)/df.residual(model3)), 3)),
    c("F-statistic (p‑value)",
      paste0(round(summary(model1)$fstatistic[1],2),
             " (p=", signif(pf(summary(model1)$fstatistic[1], summary(model1)$fstatistic[2], summary(model1)$fstatistic[3], lower.tail=FALSE),3), ")"),
      paste0(round(summary(model2)$fstatistic[1],2),
             " (p=", signif(pf(summary(model2)$fstatistic[1], summary(model2)$fstatistic[2], summary(model2)$fstatistic[3], lower.tail=FALSE),3), ")"),
      paste0(round(summary(model3)$fstatistic[1],2),
             " (p=", signif(pf(summary(model3)$fstatistic[1], summary(model3)$fstatistic[2], summary(model3)$fstatistic[3], lower.tail=FALSE),3), ")")
    )
  ),
  omit.empty = FALSE,
  notes = c("Standard errors in parentheses.",
            "Signif. codes: *p<0.1; **p<0.05; ***p<0.01")
)


######QUESTION 2#########

##Upload Libraries
library(broom)
library(dplyr)

# Construct a summary table
glance(model1) %>% select(adj.r.squared, sigma, AIC, BIC, p.value) -> m1
glance(model2) %>% select(adj.r.squared, sigma, AIC, BIC, p.value) -> m2
glance(model3) %>% select(adj.r.squared, sigma, AIC, BIC, p.value) -> m3

bind_rows(m1, m2, m3, .id = "model") %>%
  mutate(DeltaAIC = AIC - min(AIC),
         AIC_weight = exp((min(AIC) - AIC)/2) / sum(exp((min(AIC) - AIC)/2)))


######QUESTION 3#########

##Further Examine p-values of F test
# Fit your models
model1 <- lm(Vehicles.Served ~ Garage.Bays + Population, data = df)
model2 <- lm(Vehicles.Served ~ Garage.Bays + Population + Access, data = df)
model3 <- lm(Vehicles.Served ~ Garage.Bays + Population + Access + Winter, data = df)

# Custom function to extract only the F‑test p‑value
overall_p <- function(mod) {
  f <- summary(mod)$fstatistic
  p <- pf(f[1], f[2], f[3], lower.tail = FALSE)
  unname(p)
}

# Create a data frame with model names and F‑test p‑values
p_table <- data.frame(
  Model = c("Model 1", "Model 2", "Model 3"),
  Ftest_p_value = c(overall_p(model1),
                    overall_p(model2),
                    overall_p(model3))
)

# Optionally, format the p‑values
p_table$Ftest_p_value <- signif(p_table$Ftest_p_value, 4)

# View the table
print(p_table)


######QUESTION 4#########

##Examine Model 3
# Fit the model (e.g., your chosen Model 3)
model3 <- lm(Vehicles.Served ~ Garage.Bays + Population + Access + Winter, data = df)

# Extract F-test p-value
f_stat <- summary(model3)$fstatistic
p_f <- pf(f_stat[1], f_stat[2], f_stat[3], lower.tail = FALSE)
p_f  # prints the F-test p-value

# Get coefficient table with p-values
coefs <- summary(model3)$coefficients

# Extract the p-values column
p_values <- coefs[, "Pr(>|t|)"]

# Display
p_values

# Decision thresholds
alpha5 <- p_f < 0.05
alpha10 <- p_f < 0.10

# Joint significance
cat("Model joint significance: 5%?", alpha5, "| 10%?", alpha10, "\n\n")

# Check each predictor
sig5 <- p_values < 0.05
sig10 <- p_values < 0.10

data.frame(
  Predictor = names(p_values),
  p_value = round(p_values, 4),
  Significant_at_5pct = sig5,
  Significant_at_10pct = sig10
)

######QUESTION 5#########


new_obs <- data.frame(
  Garage.Bays     = 5,
  Population = 40000,
  Access     = 1,
  Winter     = 0
)


pred_val <- predict(model3, newdata = new_obs)
print(pred_val)

predict(model3, newdata = new_obs, interval = "prediction", level = 0.95)