To start read in the data in and create some models. These are multiple regressions where we can control the order and number of parameters that are fed into the model. Generally this would imply that were selecting what is the strongest predictor of your output to come first in lm(all ~ strongest + 2nd.strondgest + .....+ last.strongest). With that said we might want to make some predictions as to what predictors will be the most influential.
After running this we can call for a complete summary of any of the models with summary(fit1) or what ever model you might want a summary of. Also note that the file path has double forward slashes as in \\. R does not like single ones, so when coping the file path where your data is stored you have add extra ones in. The following code can be easly changed and explored for different paramters in the provided .csv file.
my.data <- read.csv(file="E:\\Flash Drive\\ESCI 407\\Data\\PRISM_climate_data.csv", header=TRUE, sep=",")
library(MASS)
par(mfrow=c(1:2))#set the plotting window to two figures side by side
fit1 <- lm(all ~ mint + snow + ppt, data=my.data)
fit2 <- lm(all_tsme ~ mint + snow + ppt, data=my.data)
fit3 <- lm(all_abam ~ mint + snow + ppt, data=my.data)
summary(fit1) ##
## Call:
## lm(formula = all ~ mint + snow + ppt, data = my.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.104656 -0.040288 0.001378 0.036285 0.099236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.221e-01 1.559e-01 5.915 4.98e-06 ***
## mint 2.173e-02 1.683e-02 1.291 0.209
## snow -1.033e-04 2.198e-04 -0.470 0.643
## ppt -3.508e-05 3.473e-05 -1.010 0.323
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05822 on 23 degrees of freedom
## Multiple R-squared: 0.1609, Adjusted R-squared: 0.05145
## F-statistic: 1.47 on 3 and 23 DF, p-value: 0.2488Above you can call for the output of which ever multiple regression you would like with the summary(fit#) you choose. For this I think were mostly interested in the model coefficients which we can extract with fit1$coefficients.
fit1$coefficients ## (Intercept) mint snow ppt
## 9.220852e-01 2.173229e-02 -1.032519e-04 -3.508386e-05So based on \(all = (a)mint + (b)snow + (c)ppt\) we can say that the this particular model would be roughly…..
\(all = (0.922)mint - (0.022)snow - (0.00003)ppt\)
In the chunk of code below, we are letting R run a bunch of simulations to determine what it thinks is the most parsimonious model given all the predictors you’ve fed into it. It will then systematically remove the least predictive components one by one and spit out what it thinks is the best model. There are criteria by which it decides to throw out predictors but that can get complicated so will leave it at the default. Also the AIC score is basically a measure of how well each model fits the data relative the other models. The AIC sore wants to be low but again is only relative to the other AIC scores. Step-wise regression is typically only advised with cation. The final result at the bottom here gives the most “parsimonious” model as \(all = ppt + dewpt\) with the listed coefficients. Essentially R decided that the water year the dewpt were the strongest predictive factors and the others are adding more noise.
Lots of script below, but it keeps paring down to a final multiple regression.
fit1 <- lm(all ~ meant + snow + ppt + maxvpd + minvpd + dewpt
, data=my.data)
step <- stepAIC(fit1, direction="both")## Start: AIC=-146.04
## all ~ meant + snow + ppt + maxvpd + minvpd + dewpt
##
## Df Sum of Sq RSS AIC
## - meant 1 0.0002340 0.072212 -147.95
## - maxvpd 1 0.0004018 0.072380 -147.88
## - snow 1 0.0005332 0.072511 -147.84
## - minvpd 1 0.0009460 0.072924 -147.68
## - dewpt 1 0.0029007 0.074879 -146.97
## - ppt 1 0.0053600 0.077338 -146.10
## <none> 0.071978 -146.03
##
## Step: AIC=-147.95
## all ~ snow + ppt + maxvpd + minvpd + dewpt
##
## Df Sum of Sq RSS AIC
## - maxvpd 1 0.0003577 0.072570 -149.81
## - snow 1 0.0005615 0.072774 -149.74
## - minvpd 1 0.0007398 0.072952 -149.67
## - ppt 1 0.0051751 0.077387 -148.08
## <none> 0.072212 -147.95
## - dewpt 1 0.0107338 0.082946 -146.21
## + meant 1 0.0002340 0.071978 -146.03
##
## Step: AIC=-149.81
## all ~ snow + ppt + minvpd + dewpt
##
## Df Sum of Sq RSS AIC
## - snow 1 0.0007223 0.073292 -151.55
## - minvpd 1 0.0019247 0.074495 -151.11
## <none> 0.072570 -149.81
## - ppt 1 0.0056375 0.078207 -149.79
## - dewpt 1 0.0108248 0.083395 -148.06
## + maxvpd 1 0.0003577 0.072212 -147.95
## + meant 1 0.0001900 0.072380 -147.88
##
## Step: AIC=-151.55
## all ~ ppt + minvpd + dewpt
##
## Df Sum of Sq RSS AIC
## - minvpd 1 0.0021209 0.075413 -152.78
## <none> 0.073292 -151.55
## - ppt 1 0.0082522 0.081544 -150.67
## + snow 1 0.0007223 0.072570 -149.81
## + maxvpd 1 0.0005185 0.072774 -149.74
## + meant 1 0.0002986 0.072994 -149.66
## - dewpt 1 0.0124801 0.085772 -149.30
##
## Step: AIC=-152.78
## all ~ ppt + dewpt
##
## Df Sum of Sq RSS AIC
## <none> 0.075413 -152.78
## + minvpd 1 0.0021209 0.073292 -151.55
## - ppt 1 0.0095850 0.084998 -151.55
## + maxvpd 1 0.0019204 0.073493 -151.47
## + meant 1 0.0018782 0.073535 -151.46
## - dewpt 1 0.0106014 0.086014 -151.22
## + snow 1 0.0009184 0.074495 -151.11step$anova # display results ## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## all ~ meant + snow + ppt + maxvpd + minvpd + dewpt
##
## Final Model:
## all ~ ppt + dewpt
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 20 0.07197822 -146.0352
## 2 - meant 1 0.0002339726 21 0.07221219 -147.9475
## 3 - maxvpd 1 0.0003577500 22 0.07256994 -149.8141
## 4 - snow 1 0.0007222522 23 0.07329219 -151.5467
## 5 - minvpd 1 0.0021208512 24 0.07541305 -152.7765