wind <- read.csv("C:/Users/justt/Desktop/School/622/Homework/HW 2/data_Windmill.csv")
summary(wind)
## Velocity Output
## Min. : 2.450 Min. :0.123
## 1st Qu.: 3.950 1st Qu.:1.144
## Median : 6.000 Median :1.800
## Mean : 6.132 Mean :1.610
## 3rd Qu.: 8.150 3rd Qu.:2.166
## Max. :10.200 Max. :2.386
Y <- wind$Output
X <- wind$Velocity
quad_mod <- lm(Y~X + I(X^2), data = wind)
summary(quad_mod)
##
## Call:
## lm(formula = Y ~ X + I(X^2), data = wind)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26347 -0.02537 0.01264 0.03908 0.19903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.155898 0.174650 -6.618 1.18e-06 ***
## X 0.722936 0.061425 11.769 5.77e-11 ***
## I(X^2) -0.038121 0.004797 -7.947 6.59e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1227 on 22 degrees of freedom
## Multiple R-squared: 0.9676, Adjusted R-squared: 0.9646
## F-statistic: 328.3 on 2 and 22 DF, p-value: < 2.2e-16
plot(X, Y, xlab = "Velocity", main = "Quadratic", ylab = "Output", cex = 0.5)
lines(X, quad_mod$fitted.values, col="blue", lwd = 1)
b_sp <- lm(Y~bs(X, df = 4), data = wind)
summary(b_sp)
##
## Call:
## lm(formula = Y ~ bs(X, df = 4), data = wind)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16694 -0.06493 0.01365 0.06848 0.11538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.18646 0.07063 2.640 0.0157 *
## bs(X, df = 4)1 1.18434 0.15531 7.626 2.42e-07 ***
## bs(X, df = 4)2 1.57874 0.14060 11.229 4.36e-10 ***
## bs(X, df = 4)3 2.15893 0.15095 14.302 5.78e-12 ***
## bs(X, df = 4)4 2.10989 0.09407 22.430 1.20e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09676 on 20 degrees of freedom
## Multiple R-squared: 0.9817, Adjusted R-squared: 0.978
## F-statistic: 267.7 on 4 and 20 DF, p-value: < 2.2e-16
bsp_interval <- predict(b_sp, interval = "confidence", level = 0.99)
bsp_interval
## fit lwr upr
## 1 0.1864577 -0.01449874 0.3874141
## 2 0.4221021 0.28329480 0.5609094
## 3 0.5905826 0.47867879 0.7024865
## 4 0.7059901 0.60179775 0.8101825
## 5 0.9416193 0.83260905 1.0506296
## 6 1.0569047 0.94234145 1.1714679
## 7 1.2288118 1.11066078 1.3469628
## 8 1.2920829 1.17521044 1.4089555
## 9 1.4654803 1.35998979 1.5709708
## 10 1.5714619 1.47331772 1.6696061
## 11 1.6679435 1.56731392 1.7685730
## 12 1.7341634 1.62722791 1.8410988
## 13 1.7711014 1.66171542 1.8804874
## 14 1.8085503 1.69879563 1.9183050
## 15 1.8369613 1.72818576 1.9457368
## 16 1.9603534 1.85905310 2.0616537
## 17 2.0336677 1.93019605 2.1371394
## 18 2.1105235 1.99595708 2.2250900
## 19 2.1569932 2.03393302 2.2800535
## 20 2.2397156 2.11146810 2.3679630
## 21 2.2675948 2.14644449 2.3887452
## 22 2.2945499 2.18453279 2.4045671
## 23 2.2991143 2.18520514 2.4130234
## 24 2.3009263 2.15230667 2.4495460
## 25 2.2963479 2.10240137 2.4902944
plot(X, Y, xlab = "Velocity", main = "Regression Spline with df = 4", ylab = "Output", cex = .5)
lines(X, b_sp$fitted.values, col = "blue", lwd = 1)
lines(X, bsp_interval[,2], col = "red", lty = 2)
lines(X, bsp_interval[,3], col = "red", lty = 2)
anova(quad_mod, b_sp)
## Analysis of Variance Table
##
## Model 1: Y ~ X + I(X^2)
## Model 2: Y ~ bs(X, df = 4)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 22 0.33108
## 2 20 0.18724 2 0.14384 7.6822 0.003347 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The model with the smallest sum of squared errors is the b spline model has a lower RSS of 0.18724 than the quadratic model’s RSS of 0.33108.
There is a statistically significant difference between the predictive performance of the two models as the P-value is 0.003347, which is less than 0.05.
The model that is best at predicting Output is the b spline model because the RSS is lower than that of the quadratic model, and it also has an Adjusted R-squared of 0.978, 98.7% which is a pretty large indicator that this is good at predicting the Output variable.
b_sp2 <- lm(Y~bs(X, df = 10), data = wind)
anova(b_sp, b_sp2)
## Analysis of Variance Table
##
## Model 1: Y ~ bs(X, df = 4)
## Model 2: Y ~ bs(X, df = 10)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 20 0.18724
## 2 14 0.17660 6 0.010637 0.1405 0.9881
plot(X, Y, main = "B Spline with df = 10", xlab = "Times", ylab = "Acceleration")
lines(X, b_sp2$fitted.values, col = "blue")
b_sp2_ci <- predict(b_sp2, interval = "confidence", level = 0.95)
lines(X, b_sp2_ci[,2], col = "red", lty = 2)
lines(X, b_sp2_ci[,3], col = "red", lty = 2)
The P-value is 0.9881, larger than 0.05, which shows that there is not a statistically significant difference between the predictive performance of the two models.
It also appears that this model with the 10 degrees of freedom is an overfit to the model.
c_sp = lm(Y~ns(X, df = 6), data = wind)
anova(b_sp, c_sp)
## Analysis of Variance Table
##
## Model 1: Y ~ bs(X, df = 4)
## Model 2: Y ~ ns(X, df = 6)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 20 0.18724
## 2 18 0.18604 2 0.0011954 0.0578 0.944
csp_interval <- predict(c_sp, interval = "confidence", level = 0.85)
csp_interval
## fit lwr upr
## 1 0.1813098 0.06231097 0.3003086
## 2 0.4105225 0.33124873 0.4897963
## 3 0.5870989 0.52131612 0.6528817
## 4 0.7121712 0.64428747 0.7800550
## 5 0.9662633 0.88576389 1.0467628
## 6 1.0809850 1.00264919 1.1593208
## 7 1.2350628 1.16583929 1.3042863
## 8 1.2874445 1.21642980 1.3584592
## 9 1.4309811 1.33798159 1.5239806
## 10 1.5378558 1.45190860 1.6238029
## 11 1.6544810 1.58615440 1.7228076
## 12 1.7389925 1.66827392 1.8097112
## 13 1.7836864 1.71003464 1.8573381
## 14 1.8252945 1.75235358 1.8982354
## 15 1.8546607 1.78385263 1.9254689
## 16 1.9689915 1.89076504 2.0472179
## 17 2.0340475 1.94481569 2.1232793
## 18 2.1063603 2.02573978 2.1869809
## 19 2.1521916 2.07838469 2.2259985
## 20 2.2358039 2.14908546 2.3225223
## 21 2.2633765 2.17576106 2.3509919
## 22 2.2890786 2.22021602 2.3579413
## 23 2.2944607 2.22821502 2.3607063
## 24 2.3023578 2.21711962 2.3875961
## 25 2.3065215 2.19394116 2.4191019
plot(X, Y, xlab = "Velocity", main = "Natural Cubic Spline with df = 6", ylab = "Output", cex = .5)
lines(X, c_sp$fitted.values, col = "blue", lwd = 1)
lines(X, csp_interval[,2], col = "red", lty = 2)
lines(X, csp_interval[,3], col = "red", lty = 2)