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.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Homework 5:

  1. 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 °

  2. Briefly interpret the rest of the summary outputs (as many as you can) of the lm() function from 1.

  3. 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-16
confint(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 is between -6.33 and 9.21
#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 ' ' 1
#Interpretation: 
#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-10
#Interpretation:
#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.2944
#Interpretation:
#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-16
#Interpretation:
#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 pat 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.6861
predict(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.285
## MCMC Burnin    =   1000                                         R-squared       =   0.9997
## MCMC Thinning  =      5                                         WAIC            =   26.985
## 
## -------------+-----------------------------------------------------------------------------
##    Parameter |  mean(Coef)  std(Coef)    [95% Cred. Interval]      tStat    Rank        ESS
## -------------+-----------------------------------------------------------------------------
##            x |     9.25419    0.03907      9.13062    9.28574    236.843       1 **    1000
##  (Intercept) |    -6.32633    1.90204    -10.05381   -2.58985          .       .          .
## -------------+-----------------------------------------------------------------------------
predict(mod1, data, bayes.avg=TRUE)
##           pred  se(pred)  CI 2.5% CI 97.5%
##  [1,] 187.0336 1.1655929 184.6626 189.3837
##  [2,] 214.6565 1.0704179 212.4657 216.8296
##  [3,] 288.3174 0.8441829 286.6741 290.0415
##  [4,] 426.4316 0.6536076 425.1558 427.8086
##  [5,] 454.0545 0.6727234 452.7210 455.4804
##  [6,] 536.9230 0.8327679 535.1700 538.5613
##  [7,] 619.7916 1.0870945 617.4939 621.9883
##  [8,] 675.0373 1.2824557 672.4504 677.6119
##  [9,] 564.5459 0.9103730 562.6021 566.4041
## [10,] 454.0545 0.6727234 452.7210 455.4804
## [11,] 371.1859 0.6774493 369.8937 372.5044
## [12,] 269.9022 0.8958860 268.1171 271.7408
#Student-t Ridge Prior
rv.t <- bayesreg(y~.,data, model = "t", prior = "ridge", t.dof = 5, n.samples = 1e4)
summary(rv.t)
## ==========================================================================================
## |                   Bayesian Penalised Regression Estimation ver. 1.3                    |
## |                     (c) Enes Makalic, Daniel F Schmidt. 2016-2024                      |
## ==========================================================================================
## Bayesian Student-t (DOF = 5) ridge regression                   Number of obs   =       12
##                                                                 Number of vars  =        1
## MCMC Samples   =  10000                                         std(Error)      =   2.5555
## MCMC Burnin    =   1000                                         R-squared       =   0.9999
## MCMC Thinning  =      5                                         WAIC            =   27.224
## 
## -------------+-----------------------------------------------------------------------------
##    Parameter |  mean(Coef)  std(Coef)    [95% Cred. Interval]      tStat    Rank        ESS
## -------------+-----------------------------------------------------------------------------
##            x |     9.22342    0.03772      9.14034    9.29187    244.520       1 **    9990
##  (Intercept) |    -6.99100    1.84518    -10.54652   -3.09923          .       .          .
## -------------+-----------------------------------------------------------------------------
predict(rv.t, data, bayes.avg=TRUE)
##           pred  se(pred)  CI 2.5% CI 97.5%
##  [1,] 186.6090 1.1322667 184.4351 188.9911
##  [2,] 214.2661 1.0398247 212.2576 216.4473
##  [3,] 288.0185 0.8191446 286.4272 289.7472
##  [4,] 426.3042 0.6266410 425.1038 427.6090
##  [5,] 453.9613 0.6433102 452.7198 455.2709
##  [6,] 536.9327 0.7943681 535.3227 538.5025
##  [7,] 619.9042 1.0386010 617.7752 621.9421
##  [8,] 675.2184 1.2268636 672.7212 677.6276
##  [9,] 564.5899 0.8686965 562.8035 566.2992
## [10,] 453.9613 0.6433102 452.7198 455.2709
## [11,] 370.9899 0.6535853 369.7253 372.3499
## [12,] 269.5804 0.8697605 267.8910 271.4174
#Plot the different estimates with confidence and credible intervals
plot(x, y, xlab="x", ylab="y", frame.plot=FALSE)
#lines(sort(x),yH, col="red", lwd=2.5) #Least squares fit
 
#no duplicates in x-values
newx = data.frame(x=seq(from=min(x), to= max(x), length.out=20))
dat2 = newx  
#newx = data.frame(x=x)
 
#95% credible interval
c.lim=as.data.frame(predict(fit1, newx , level=0.95, interval="conf")) #
 
#Least squares
abline(fit1, col="red")
lines(cbind(newx,c.lim$lwr), col="red", lty="dashed")
lines(cbind(newx,c.lim$upr), col="red", lty="dashed")
 
dat2<-newx

#Guassian horseshoe 
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")
 
## Gaussian ridge
#yhat_G <- predict(rv.G, dat2, bayes.avg=TRUE)
#lines(dat2$x, yhat_G[,1], col="blue", lwd=1)
#lines(dat2$x, yhat_G[,3], col="blue", lwd=1, lty="dashed")
#lines(dat2$x, yhat_G[,4], col="blue", lwd=1, lty="dashed")
 
# Student-t ridge
yhat_t <- predict(rv.t, dat2, bayes.avg=TRUE)
lines(dat2$x, yhat_t[,1], col="green", lwd=1)
lines(dat2$x, yhat_t[,3], col="green", lwd=1, lty="dashed")
lines(dat2$x, yhat_t[,4], col="green", lwd=1, lty="dashed")
 
legend("topleft",c("Least Squares", "Horseshoe","Student-t (dof=5)"),lty=c(1,1,1),col=c("red", "black","green"),
       lwd=c(2.5,2.5,2.5), cex=0.7)

#predict at new X=58 (matches Part D's temperature of interest)
dat3 = data.frame(x=58)
 
#Least squares point prediction at x=58
predict(fit1, dat3)
##       1 
## 527.759
#Least squares confidence interval
predict(fit1, dat3, level=0.95, interval="confidence")
##       fit      lwr      upr
## 1 527.759 526.2368 529.2813
#Horseshoe prediction at x=58
predict(mod1, dat3, bayes.avg=TRUE)
##          pred  se(pred)  CI 2.5% CI 97.5%
## [1,] 527.7154 0.8090222 526.0385 529.3278
#Student-t ridge prediction at x=58
predict(rv.t, dat3, bayes.avg=TRUE)
##          pred  se(pred)  CI 2.5% CI 97.5%
## [1,] 527.7137 0.7716915 526.1526 529.2388
#predict(rv.G, dat3, bayes.avg=TRUE)
 
#Interpretation:
#Point predictions at x=58 are nearly identical across all three methods: 
#527.76 (least squares), 527.773 (horseshoe), 527.71 (Student-t ridge) --
#differences are small. This makes sense given the extremely strong 
#linear fit (R-squared = 0.9999), which leaves little room for the Bayesian 
#priors to shift the estimates away from least squares.
#
#Interval widths are also similar: least squares' 95% CI (width ~3.04) is
#close to the horseshoe (~3.258) and Student-t ridge (~3.158) credible intervals.
#The least-squares 95% prediction interval (width ~9.19) is much wider, since
#it accounts for individual observation variability, not just estimation
#uncertainty so it isn't directly comparable to the credible intervals above.
#Overall: least squares and both Bayesian methods agree closely here, since 
#the strong signal in the data leaves minimal room for the priors to matter.