#create the cross validation setupcross_val <-trainControl(method ='cv', number =10) #create the different models 1-15 degree to fid optmum degreeupper_lim_deg =15#theese objects will store ersults and the models for the anova copmarisonresults_list <-vector("list", upper_lim_deg)models_list <-vector("list", upper_lim_deg)#Train the modelsfor(degree in1:upper_lim_deg) {#create the template of the formula functions <-as.formula(paste0('wage~poly(age, degree=',degree,')')) #fit the models poly_model<-train( functions, data=Wage,method ='lm',trControl=cross_val )#collect all of the pieces together in the resulting dataframe results_list[[degree]]<-data.frame(Degree= degree,RMSE=poly_model$results$RMSE,R_SQR= poly_model$results$Rsquared,MAE= poly_model$results$MAE )#store al o fthe models into a dataframe of the results models_list[[degree]]<- poly_model$finalModel}#responses for the metrics on the models given response_totals_df<-do.call(rbind, results_list)print(response_totals_df)
#compare agains teh anova o the modelsfor(degree in1:(upper_lim_deg-1)){ anova_comp <-anova(models_list[[degree]], models_list[[degree+1]]) pval=anova_comp$"Pr(>F)"[2]cat('Anova compare degree:',degree,' and ', degree+1,' has p val:', pval,'\n')}
Anova compare degree: 1 and 2 has p val: 3.07742e-32
Anova compare degree: 2 and 3 has p val: 0.001687063
Anova compare degree: 3 and 4 has p val: 0.05103865
Anova compare degree: 4 and 5 has p val: 0.369682
Anova compare degree: 5 and 6 has p val: 0.1162015
Anova compare degree: 6 and 7 has p val: 0.205311
Anova compare degree: 7 and 8 has p val: 0.7779522
Anova compare degree: 8 and 9 has p val: 0.03596326
Anova compare degree: 9 and 10 has p val: 0.9675292
Anova compare degree: 10 and 11 has p val: 0.7990375
Anova compare degree: 11 and 12 has p val: 0.9479033
Anova compare degree: 12 and 13 has p val: 0.7301651
Anova compare degree: 13 and 14 has p val: 0.6964635
Anova compare degree: 14 and 15 has p val: 0.4801506
It is a unique results where the anova results in a degree of 3 with the pvalue = .051 when running an anova comparison to degree = 4. However with cross validation using 10 fold - we get that degree 8 has the lowest RMSE. The anova sheds insight that the comparison of degree 8 and 9 is significant, however this can be due to random noise, which can influence the cv results. By the rule of parsimony we will go with the simplest model, therefore degree 3.
ggplot(data= Wage, aes(x= age, y= wage)) +geom_point(color="black", alpha=0.3, size =1.2) +stat_smooth(method="lm", formula= y ~poly(x, 3), se =TRUE, color="red",linewidth=1.2 ) +labs(title="Deg3 Polynomial Regression Fit: Wage vs. Age",x="Age (Years)",y="Wage" ) +theme_minimal()
The optimum number of cuts was chosen to be 5. It is true that degree poly 8 has the lowest RMSE, however the rate of return after each subsequent cut is not justified to get that many number of cuts. When compared to the polynomial model of degree 3 (RMSE=39.87239) the 5 cut yields an RMSE of (40.21877), only slightly being out beat by the poly model.
#set up cross alidation and cut params cross_v<-trainControl(method='cv',n=10)min =2max =15results<-vector('list', 14)#fit and create the model for(k in min:max){#create step formula step_function<-as.formula(paste0('wage~cut(age,',k,',)'))#trainining step_model <-train( step_function,data=Wage,method='lm',trControl=cross_v )#save to dataframe results[[k-1]] <-data.frame(Cuts= k,RMSE= step_model$results$RMSE,R_SQR= step_model$results$Rsquared,MAE = step_model$results$MAE )}step_totals_df<-do.call(rbind,results)print(step_totals_df)
library(leaps)set.seed(1)attach(College)train <-sample(length(Outstate), length(Outstate) /2)test <--trainCollege.train <- College[train, ]College.test <- College[test, ]fit <-regsubsets(Outstate ~ ., data = College.train, nvmax =17, method ="forward")fit.summary <-summary(fit)par(mfrow =c(1, 3))plot(fit.summary$cp, xlab ="Number of variables", ylab ="Cp", type ="l")min.cp <-min(fit.summary$cp)std.cp <-sd(fit.summary$cp)abline(h = min.cp +0.2* std.cp, col ="red", lty =2)abline(h = min.cp -0.2* std.cp, col ="red", lty =2)plot(fit.summary$bic, xlab ="Number of variables", ylab ="BIC", type='l')min.bic <-min(fit.summary$bic)std.bic <-sd(fit.summary$bic)abline(h = min.bic +0.2* std.bic, col ="red", lty =2)abline(h = min.bic -0.2* std.bic, col ="red", lty =2)plot(fit.summary$adjr2, xlab ="Number of variables", ylab ="Adjusted R2", type ="l", ylim =c(0.4, 0.84))max.adjr2 <-max(fit.summary$adjr2)std.adjr2 <-sd(fit.summary$adjr2)abline(h = max.adjr2 +0.2* std.adjr2, col ="red", lty =2)abline(h = max.adjr2 -0.2* std.adjr2, col ="red", lty =2)
#from table above 6 is the minimum number of optimum variableslm1 <-regsubsets(Outstate ~ ., data = College, method ="forward")coeffs <-coef(fit, id =6)names(coeffs)
From the graphs above, Board, alumni, and rate all show a linear pattern, meanwhile the Expend variables shows a very strong indication of being nonlinear with the curve. The variable phD is a very weak nonlinear pattern with a slight curve upward.
predict.regsubsets <-function(object, newdata, id) { form <-as.formula(object$call[[2]]) mat <-model.matrix(form, newdata) coefi <-coef(object, id = id) xvars <-names(coefi) mat[, xvars] %*% coefi}preds <-predict(fit, College.test, id=6)err <-mean((College.test$Outstate - preds)^2)cat('\nError Term:',err,'\n')
Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(PhD,
df = 2) + s(perc.alumni, df = 2) + s(Expend, df = 5) + s(Grad.Rate,
df = 2), data = College.train)
Deviance Residuals:
Min 1Q Median 3Q Max
-7402.89 -1114.45 -12.67 1282.69 7470.60
(Dispersion Parameter for gaussian family taken to be 3711182)
Null Deviance: 6989966760 on 387 degrees of freedom
Residual Deviance: 1384271126 on 373 degrees of freedom
AIC: 6987.021
Number of Local Scoring Iterations: NA
Anova for Parametric Effects
Df Sum Sq Mean Sq F value Pr(>F)
Private 1 1778718277 1778718277 479.286 < 2.2e-16 ***
s(Room.Board, df = 2) 1 1577115244 1577115244 424.963 < 2.2e-16 ***
s(PhD, df = 2) 1 322431195 322431195 86.881 < 2.2e-16 ***
s(perc.alumni, df = 2) 1 336869281 336869281 90.771 < 2.2e-16 ***
s(Expend, df = 5) 1 530538753 530538753 142.957 < 2.2e-16 ***
s(Grad.Rate, df = 2) 1 86504998 86504998 23.309 2.016e-06 ***
Residuals 373 1384271126 3711182
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova for Nonparametric Effects
Npar Df Npar F Pr(F)
(Intercept)
Private
s(Room.Board, df = 2) 1 1.9157 0.1672
s(PhD, df = 2) 1 0.9699 0.3253
s(perc.alumni, df = 2) 1 0.1859 0.6666
s(Expend, df = 5) 4 20.5075 2.665e-15 ***
s(Grad.Rate, df = 2) 1 0.5702 0.4506
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
From the Anova parametrics table all pvalues being below .05 means all vairables should be included in the model since they are all significant. Fro the non-paramteric table, only Expend has a pvalue below .05 so we reject the null hypothesis that, Expend is a strictly linear relationship, and all other variables fail to reject null hypothesis.
