Montiverdi Homework 5:
2026-07-12
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
Homework 5:
Problem 2.12 parts a-d. Use shortcuts in R whenever you can. Interpret results in the context of the problem. A: Fit a simple linear regression model to the data B: Test for significance of regression C: Plant management believes that an increase in average ambient temperature of 1 degree will increase average monthly steam consumption by 10,000 lb. Do the data support this statement? D: Construct a 99% prediction interval on steam usage in a month with average ambient temperature of 58 °
Briefly interpret the rest of the summary outputs (as many as you can) of the lm() function from 1.
Compare results using the bayesreg package. See corresponding R file for examples.
# Read data
data<-read.csv("data-prob-2-12.csv")
data## temp usage
## 1 21 185.79
## 2 24 214.47
## 3 32 288.03
## 4 47 424.84
## 5 50 454.68
## 6 59 539.03
## 7 68 621.55
## 8 74 675.06
## 9 62 562.03
## 10 50 452.93
## 11 41 369.95
## 12 30 273.98#Part A: Fit a simple linear regression model to the data (t approach)
#Assign x and y
x=data[,1]
y=data[,2]
n=length(x)
#Model and Summary
fit1=lm(y~x)
summary(fit1)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5629 -1.2581 -0.2550 0.8681 4.0581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.33209 1.67005 -3.792 0.00353 **
## x 9.20847 0.03382 272.255 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.946 on 10 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 7.412e+04 on 1 and 10 DF, p-value: < 2.2e-16confint(fit1)## 2.5 % 97.5 %
## (Intercept) -10.053181 -2.610993
## x 9.133106 9.283830#Plot data
plot(x,y)
abline(fit1, col = "red")
Interpretation: Fitted model: yhat = -6.33 + 9.21x
We are 95% confident the true slope (change in average steam usage per 1-degree increase in temperature) lies between 9.133 and 9.284 (thousand lb). Since this interval is entirely positive and narrow, it confirms a strong, estimated positive relationship between temperature and steam usage.
We are also 95% confident the true intercept lies between -10.053 and -2.611,though since a temperature of 0 degrees is outside the observed data range, this interval has little practical meaning.
Slope (9.21): for every 1-degree increase in average temperature, average monthly steam usage is estimated to increase by about 9.21 (thousand lb), on average.
Intercept (-6.33): the model’s estimated usage when temperature = 0 degrees
#Part B: Test for significance of regression (f approach)
#Anova
anova(fit1) ## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 280590 280590 74123 < 2.2e-16 ***
## Residuals 10 38 4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretation: The f-value is 74,124 and the p-value < 0.0001 which is lower than the standard alpha of 0.05. H0: B1 = 0 (no linear relationship) vs Ha: B1 does not = 0. In which case we will reject the null hypothesis because there is a significant linear relationship between x (temp) and y (usage).
#Part C: Plant management believes that an increase in average ambient temperature of 1 degree will #increase average monthly steam consumption by 10,000 lb. Does the data support this statement?
#Slope
B1H = coef(fit1)["x"]
#Standard error
se_B1H = summary(fit1)$coefficients["x", "Std. Error"]
#T-statistic
t_stat = (B1H - 10) / se_B1H
#Degrees of Freedom Residual
df = fit1$df.residual
#P-value
p_val = 2 * pt(-abs(t_stat), df)
t_stat; p_val## x
## -23.40222## x
## 4.597358e-10Interpretation: H0: beta1 = 10 (a 1-degree increase in temp increases usage by 10,000 lb) Ha: beta1 != 10 with t = -23.41, p < .0001, we reject H0 at alpha = 0.05. The data does NOT support management’s claim. The true slope is statistically significantly different from 10 (10,000 lb per degree). Since our estimated slope (B1H = 9.21) is significantly less than 10, the data suggest that a 1-degree increase in temperature is associated with a smaller increase in steam usage than management believes.
#Part D: Construct a 99% prediction interval on steam usage in a month with average ambient #temperature of 58 °
predict(fit1, data.frame(x = 58), level = 0.99, interval = "predict")## fit lwr upr
## 1 527.759 521.2237 534.2944Interpretation: We predict with 99% confidence that steam usage in a single month with average ambient temperature of 58 degrees will fall between 521.22 and 534.29 (thousand lb), with a point prediction (best estimate) of 527.76 (thousand lb).
#Checking fit for normalcy
qqnorm(resid(fit1))
qqline(resid(fit1))
Interpretation: If residuals fall approximately along the diagonal reference line, the normality assumption for inference (t-tests, prediction intervals) is reasonably satisfied.
#Question 2: Briefly interpret the rest of the summary outputs (as many as you can) of the lm() #function from 1.
summary(fit1)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5629 -1.2581 -0.2550 0.8681 4.0581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.33209 1.67005 -3.792 0.00353 **
## x 9.20847 0.03382 272.255 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.946 on 10 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 7.412e+04 on 1 and 10 DF, p-value: < 2.2e-16Interpretation: Residuals Min:-2.5629, 1Q:-1.2581, Median:-0.2550, 3Q:0.8681, Max:4.0581 The residuals are reasonably centered around 0 and fairly symmetric, though the max (4.0581) is somewhat larger than the min (-2.5629), suggesting a slight right skew a possible mild-outlier on the high side. Overall the spread looks reasonable for a well-fitting linear model.
Intercept significance (-6.332, t= -3.792, p = 0.00353): The intercept (-6.33) is statistically significant (p < 0.01), meaning it’s reliably different from 0. However, since a temperature of 0 degrees is far outside the observed data range, the intercept has little interpretation here.
Slope significance (9.208, t= 272.255, p < 2e-16): This confirms the same conclusion as Part B: temperature is an extremely significant predictor of steam usage (t = 272.26, p < 0.0001).
Residual standard error (1.946 on 10 degrees of freedom): The typical deviation of observed steam usage from the model’s predicted values is about 1.946 (thousand lb). This gives a sense of the model’s typical prediction error in the original units.
Multiple R-squared (0.9999): About 99.99% of the variability in steam usage is explained by temperature alone. This is an extremely strong fit due to the value being significantly close to 1.0 near perfect linear predictor of usage in this dataset.
Adjusted R-squared (0.9999): Nearly identical to R-squared, as expected with only one predictor showing accuracy of the assumption
F-statistic (7.412e+04, p < 2.2e-16): Reconfirms Part B’s ANOVA result – overall the regression model is highly statistically significant.
#Question 3:Compare results using the bayesreg package. See corresponding R file for examples.
#Using shortcuts and bayesreg package
#No p-values (no statistical significance testing) for testing effects of regressors
#Rename
names(data)[names(data) == "temp"] <- "x"
names(data)[names(data) == "usage"] <- "y"
names(data)## [1] "x" "y"#BAYESIAN regression
library(bayesreg)## Loading required package: pgdraw## Loading required package: doParallel## Loading required package: foreach## Loading required package: iterators## Loading required package: parallel#Least squares
predict(fit1, data,level=0.95, interval="conf")## fit lwr upr
## 1 187.0457 184.7524 189.3390
## 2 214.6711 212.5637 216.7786
## 3 288.3389 286.6775 290.0003
## 4 426.4659 425.2139 427.7179
## 5 454.0913 452.8124 455.3702
## 6 536.9675 535.4011 538.5339
## 7 619.8437 617.7964 621.8910
## 8 675.0945 672.6735 677.5155
## 9 564.5929 562.8810 566.3048
## 10 454.0913 452.8124 455.3702
## 11 371.2151 369.8968 372.5334
## 12 269.9219 268.1578 271.6861predict(fit1, data,level=0.95, interval="pred")## fit lwr upr
## 1 187.0457 182.1414 191.9501
## 2 214.6711 209.8509 219.4914
## 3 288.3389 283.6963 292.9815
## 4 426.4659 421.9536 430.9782
## 5 454.0913 449.5715 458.6112
## 6 536.9675 532.3581 541.5769
## 7 619.8437 615.0495 624.6380
## 8 675.0945 670.1292 680.0599
## 9 564.5929 559.9320 569.2538
## 10 454.0913 449.5715 458.6112
## 11 371.2151 366.6839 375.7462
## 12 269.9219 265.2416 274.6023#Horseshoe Prior
mod1<-bayesreg(y~.,data, prior="hs", n.samples = 1e3) # hs prior is horseshoe
summary(mod1)## ==========================================================================================
## | Bayesian Penalised Regression Estimation ver. 1.3 |
## | (c) Enes Makalic, Daniel F Schmidt. 2016-2024 |
## ==========================================================================================
## Bayesian Gaussian horseshoe regression Number of obs = 12
## Number of vars = 1
## MCMC Samples = 1000 std(Error) = 2.308
## MCMC Burnin = 1000 R-squared = 0.9997
## MCMC Thinning = 5 WAIC = 26.983
##
## -------------+-----------------------------------------------------------------------------
## Parameter | mean(Coef) std(Coef) [95% Cred. Interval] tStat Rank ESS
## -------------+-----------------------------------------------------------------------------
## x | 9.25420 0.03958 9.12794 9.28774 233.807 1 ** 1000
## (Intercept) | -6.30482 1.95614 -10.13007 -2.10671 . . .
## -------------+-----------------------------------------------------------------------------predict(mod1, data, bayes.avg=TRUE)## pred se(pred) CI 2.5% CI 97.5%
## [1,] 187.0747 1.2040278 184.7133 189.6555
## [2,] 214.7004 1.1058646 212.5604 217.0259
## [3,] 288.3687 0.8695638 286.7349 290.1097
## [4,] 426.4970 0.6491397 425.2271 427.7533
## [5,] 454.1226 0.6625917 452.8602 455.4186
## [6,] 536.9996 0.8131165 535.3584 538.6129
## [7,] 619.8765 1.0663398 617.6800 622.0406
## [8,] 675.1278 1.2630746 672.5949 677.6556
## [9,] 564.6252 0.8897191 562.7970 566.3490
## [10,] 454.1226 0.6625917 452.8602 455.4186
## [11,] 371.2457 0.6860676 369.9078 372.4856
## [12,] 269.9516 0.9241432 268.2038 271.8505#Plot the least squares and horseshoe estimates with confidence and credible intervals
plot(x, y, xlab="x", ylab="y", frame.plot=FALSE)
newx = data.frame(x=seq(from=min(x), to=max(x), length.out=20))
dat2 = newx
#95% confidence interval (least squares)
c.lim = as.data.frame(predict(fit1, newx, level=0.95, interval="conf"))
abline(fit1, col="red")
lines(cbind(newx, c.lim$lwr), col="red", lty="dashed")
lines(cbind(newx, c.lim$upr), col="red", lty="dashed")
#Horseshoe credible interval
yhat_hs <- predict(mod1, dat2, bayes.avg=TRUE)
lines(dat2$x, yhat_hs[,1], col="black", lwd=1)
lines(dat2$x, yhat_hs[,3], col="black", lwd=1, lty="dashed")
lines(dat2$x, yhat_hs[,4], col="black", lwd=1, lty="dashed")
legend("topleft", c("Least Squares", "Horseshoe"), lty=c(1,1),
col=c("red","black"), lwd=c(2.5,2.5), cex=0.7)
#predict at new X=58 (matches Part D's temperature of interest)
dat3 = data.frame(x=58)
predict(fit1, dat3)## 1
## 527.759predict(fit1, dat3, level=0.95, interval="confidence")## fit lwr upr
## 1 527.759 526.2368 529.2813predict(mod1, dat3, bayes.avg=TRUE)## pred se(pred) CI 2.5% CI 97.5%
## [1,] 527.791 0.7899015 526.1881 529.3274Interpretation: Both models used temperature as the predictor variable. Both models had nearly identical slope estimates (lm slope estimate 9.20847, Bayesian mean estimate 9.25421). Both models indicate the slope is estimated very precisely (lm Std. Error (Slope) = 0.03382, Bayesian Std(Coef) = 0.04039). Both models produced an extremely large t value (lm 272.255) and tStat (Bayesian 229.093), indicating a very strong relationship between temperature and steam usage. Both models conclude temperature is an important predictor, but they express the evidence differently.