#------------------------------------------------------# #------- Data 2 Decisions 2018 Exam #2 Solution -------# #------------------------------------------------------# # Packages needed for the multigraphsummary routine: # install.packages("MASS") # install.packages("e1071") # install.packages("car") # install.packages("lmtest") # install.packages("lawstat") multigraphsummary <- function(model){ # This routine produces a multigraph summary of a regression model, plus # various test of the OLS assumptions. # First, let's create data frame from the model object N = nobs(model) p = length(model$coefficients) # Grab the response (y) and predictor (x) values and put them into a data frame MyData = matrix(nrow = N, ncol = p) MyData = data.frame(MyData) MyData[,1] = model$model[,1] names(MyData)[1] = paste(attributes(model$terms)$variable[2]) for (i in 2:p){ MyData[,i] = model$model[,i] names(MyData)[i] = attributes(model$terms)$term.labels[i-1] } y = MyData[,1] x = MyData[,2] cat("------------------------", "\n") cat("Multigraph Regression Summary and Diagnostics", "\n") # turn off significance stars options(show.signif.stars=F) #------------------------------------------------------------------------- print(summary(model)) print(confint(model, level=0.95)) # Akiake Information Criterion and Schwarz's Bayesian Criterion cat("\n") cat("AIC =", AIC(model, k=2)," BIC =", AIC(model, k=log(N)), "\n") # PRESS: predicted residual error sum of squares pr <- residuals(model)/(1 - lm.influence(model)$hat) PRESS <- sum(pr^2) my.anova <- anova(model) tss <- sum(my.anova$"Sum Sq") # predictive R^2 pred.r.squared <- 1 - PRESS/(tss) cat("PRESS =", PRESS, " Predictive R^2 =", pred.r.squared, "\n") cat("------------------------", "\n") #------------------------------------------------------------------------- # Check for multicollinearity of p > 2 if (p>2) { cat("Checking for Multicollinearity", "\n") cat("Correlation Matrix:", "\n") print(cor(MyData)) cat("Eigenvalues:", "\n") print(eigen(cor(MyData))$values) # Condition Number: Ratio of max to min Eigen values of the correlation matrix cat("Condition Number (is it large, > 100 or so?):", "\n") print(kappa(cor(MyData), exact = TRUE)) library(car) cat("Variance Inflation Factors (>4=start worrying; >10=do something?):", "\n") print(vif(model)) cat("Is the mean VIF much bigger than 1?", "\n") cat("mean VIF = ", mean(vif(model)), "\n") cat("------------------------", "\n") } #------------------------------------------------------------------------- # Generate some basic plots of the data if (p>2) { plot(MyData) } # set up so that four graphs appear on each page par(mfrow=c(2,2)) if (p==2) { # plot the data and the best fit model ylabel = paste(attributes(model$terms)$variable[2]) xlabel = attributes(model$terms)$term.labels[1] plot(y ~ x, cex = 1.3, main = "Data and Model", xlab=xlabel, ylab=ylabel) # now add the model line and CI to the plot # create a new data frame that is a simple array of x-values newx = seq(min(x),max(x),length.out = 101) new <- data.frame(x = newx) names(new) = xlabel y_hat = predict(model, newdata=new, interval="confidence", level = 0.95, type="response") lines(newx, y_hat[,1], col = "red") # Note: this shows the confidence intervals, not the prediction intervals lines(newx, y_hat[,2], col = "blue", lty=2) lines(newx, y_hat[,3], col = "blue", lty=2) # add a sub-title to the graph that shows the fit R^2 mylabel = paste("R^2 = ",format(summary(model)$r.squared, digits = 3)) mtext(mylabel) } #plot the predicted y-values versus the measured y-values # install.packages("pls") plot(fitted(model) ~ y, cex = 1.3, main = "Predicted Response", xlab="Measured Y", ylab="Predicted Y") abline(a=0,b=1,col='red') # a y=x line # plot the residuals library(MASS) #for the studres function esr = studres(model) # externally studentized residuals plot(esr~predict(model), main = "Residuals", xlab="Predicted Response") abline(h=0,col='red') #------------------------------------------------------------------------- # Testing for Outliers #create a box and whisker plot #Box shows Q1 to Q3, with line at median. IQR = Q3 - Q1 #The top whisker denotes the min(maximum value, Q3 + 1.5*IQR) #The bottom whisker denotes the max(minimum value, Q1 - 1.5*IQR) library(e1071) boxplot(esr, ylab = "esr", main = "Box & Wiskers") cat("Testing for Outliers:", "\n") alpha = 0.01 tinv = qt(alpha/2/N,N-p-1) Grubbs_crit = (N-1)/sqrt(N)*sqrt(tinv^2/(N-3+tinv^2)) cat("Grubbs' Critical Value (alpha = 0.01) =", Grubbs_crit, "\n") # Now test for one outlier library(car) # Bonferonni p-value (Grubbs test) for most extreme observation cat("Grubbs' test for one outlier using outlierTest from the car package:", "\n") print(outlierTest(model)) cat("------------------------", "\n") #------------------------------------------------------------------------- # Testing for Normality # histogram of externally studentized residuals hist(studres(model), freq=FALSE, main="Distribution of esr", col="grey", xlab = "esr") #add a standard normal curve on top of the historgram xfit<-seq(min(studres(model)),max(studres(model)),length=40) yfit<-dnorm(xfit) lines(xfit, yfit) # Normal Probability Plot qqout = qqnorm(esr) qqline(esr, col = "red") # calculate the correlation coefficient from the qqplot qq_correlation = cor(qqout$x,qqout$y) mylabel = paste("r = ",format(qq_correlation, digits = 4)) text(-1.5, 2,labels = mylabel) #Testing for normality cat("Testing Externally Studentized Residuals (esr) for Normality:", "\n") library(e1071) #this generates G1, the unbiased estimator for skewness G1 = skewness(esr, na.rm = TRUE, type = 2) #type 1 = biased estimator, type 2 = unbiased estimator #Now let's calculate the standard error of the unbiased skewness estimator SE_G1 = sqrt(6*N*(N-1)/(N-2)/(N+1)/(N+3)) z1 = G1/SE_G1 if (z1 > 0) { pvalue = 2*(1-pnorm(z1)) } else { pvalue = 2*pnorm(z1) } cat("(Small p-values means we can reject the assumption of esr normality)", "\n") cat("Skewness Z = ",format(z1, digits = 4), ", p-value = ", pvalue, "\n") #this generates G2, the unbiased estimator for kurtosis G2 = kurtosis(esr, na.rm = TRUE, type = 2) #type 1 = biased estimator, type 2 = unbiased estimator #Now let's calculate the standard error of the unbiased skewness estimator SE_G2 = 2*SE_G1*sqrt((N*N-1)/(N-3)/(N+5)) z2 = G2/SE_G2 if (z2 > 0) { pvalue = 2*(1-pnorm(z2)) } else { pvalue = 2*pnorm(z2) } cat("Kurtosis Z = ",format(z2, digits = 4), ", p-value = ",pvalue, "\n") #performs the Shapiro-Wilk test for normality print(shapiro.test(esr)) #reject the null hypothesis that the data is normally distributed if the #p-value is less that your significance level (e.g., 0.05) cat("------------------------", "\n") #------------------------------------------------------------------------- # Testing for Influence cat("Testing Externally Studentized Residuals (esr) for Influence:", "\n") # plotting options specifically for lm objects # plot(model, which=4, cook.levels=cutoff) # 1 = residuals versus fitted # 2 = QQ plot # 3 = Scale-Location # 4 = Cook's Distance # 5 = Williams-like Graph # Let's create a Williams Graph norm_leverage = (length(esr)/2)*hatvalues(model) plot(abs(esr) ~ norm_leverage, xlab = "leverage*n/p", ylab = "|esr|", main ="Williams Graph") # add the critical lines to the Williams Graph abline(h=Grubbs_crit, col="red") # draw a red horizontal line at Gcrit abline(v=2, col="red") # draw a red vertical line at 2 D = cooks.distance(model) #this provides the cook's Distance for each data points cutoff <- 4/summary(model)$df[2] cat("Cook's Distance cut-off (4/df) = ", cutoff, "\n") cat("Maximum Cook's Distance = ", max(D), "\n") cat(" which occurs at index = ", which.max(D), "\n") # we can also calculate DFFITS and DFBETA plot(D ~ fitted(model), ylab = "Cook's D", main = "Cook's Distance", xlab="Predicted Response") # For Cook's D plot, identify D values > 4/DF abline(h=cutoff, col="red") xrange = max(fitted(model)) - min(fitted(model)) x_text = max(fitted(model)) - 0.1*xrange y_text = cutoff + 0.05*max(D) text(x_text, y_text, "D = 4/df", col = "red") dff = dffits(model) cat("Maximum DFFITS = ", max(abs(dff)), "\n") cat(" which occurs at index = ", which.max(abs(dff)), "\n") plot(dff ~ fitted(model), ylab = "DFFITS", main = "DFFITS", xlab="Predicted Response") dfb = dfbeta(model)[,'(Intercept)'] cat("Intercept Maximum DFBETA = ", max(abs(dfb)), "\n") cat(" which occurs at index = ", which.max(abs(dfb)), "\n") plot(dfb ~ fitted(model), ylab = "DFBETA", main = "Intercept DFBETA", xlab="Predicted Response") for (i in 2:p){ dfb = dfbeta(model)[,i] label = cat(attributes(model$terms)$term.labels[i-1], "Maximum DFBETA") cat(label, "= ", max(abs(dfb)), "\n") cat(" which occurs at index = ", which.max(abs(dfb)), "\n") mainlabel = paste(attributes(model$terms)$term.labels[i-1], "DFBETA") plot(dfb ~ fitted(model), ylab = "DFBETA", main = mainlabel, xlab="Predicted Response") } cat("------------------------", "\n") #------------------------------------------------------------------------- # Testing Externally Studentized Residuals (esr) for Homoscedasticity cat("Testing for Homoscedasticity:", "\n") cat("(Data sorted by y-hat and split in half)", "\n") data_sorted = data.frame(cbind(fitted(model),studres(model))) attach(data_sorted) data_sorted = data_sorted[order(X1),] detach(data_sorted) esr_sorted = data_sorted[,2] cat("(Small p-value indicates heteroscedasticity)", "\n", "\n") # Breusch-Pagan test for nonconstant variance library(lmtest) cat("Breusch-Pagan test from bptest, lmtest package:", "\n") print(bptest(model)) #from lmtest package, defaults to to testing against explanatory variables # split the data in half N = length(esr_sorted) group = rep(2,N) group[1:N/2] = 1 cat("Barlett test from bartlett.test, stats package:", "\n") print(bartlett.test(esr_sorted, as.factor(group))) # If the p-value is less than alpha, we reject the null hypothesis that # the variances of each group are equal. library(lawstat) cat("Brown-Forsythe test from levene.test, lawstat package:", "\n") print(levene.test(esr_sorted, as.factor(group))) # Note that this test is formualted using an F statistic, rather than the # t-statistic that was discussed in class. They work out the same, giving # the same p-value. cat("------------------------", "\n") #------------------------------------------------------------------------- # Testing for Homoscedasticity rho = acf(model$residuals, main = "Autocorrelation Function") r = rho[1]$acf[1,,] # manual creation of a lag plot n = length(model$residuals) yprime = rep(0,n-1) xprime = rep(0,n-1) for (i in 1:n-1){ yprime[i] = model$residuals[i+1] xprime[i] = model$residuals[i] } plot(yprime ~ xprime, main = "Lag 1 Plot", xlab = "residual lag 1", ylab = "residual") cat("Testing for Autocorrelated Residuals:", "\n") cat("(only useful if data is in some natural order, such as time)", "\n") #Durbin-Watson Test library(lmtest) print(dwtest(model)) cat("Estimated AR(1) parameter r = ", r, "\n") cat("------------------------", "\n") return } #Importing data: use .csv format if at all possible WineQuality <- read.csv('WineQuality.csv') # best main factor model model1 = lm(quality~volatile.acidity + chlorides + free.sulfur.dioxide + total.sulfur.dioxide + pH + sulphates + alcohol, data = WineQuality) multigraphsummary(model1) print(summary(model1)) cat("Model 1 AIC =", AIC(model1, k=2), "\n") # best main factor + interaction model (^2) model2 = lm(quality~fixed.acidity + volatile.acidity + citric.acid + chlorides + free.sulfur.dioxide + total.sulfur.dioxide + pH + sulphates + alcohol + fixed.acidity:citric.acid + fixed.acidity:chlorides + fixed.acidity:free.sulfur.dioxide + fixed.acidity:pH + fixed.acidity:sulphates + fixed.acidity:alcohol + volatile.acidity:citric.acid + volatile.acidity:free.sulfur.dioxide + volatile.acidity:total.sulfur.dioxide + volatile.acidity:alcohol + citric.acid:total.sulfur.dioxide + citric.acid:density + citric.acid:pH + citric.acid:alcohol + residual.sugar:total.sulfur.dioxide + chlorides:free.sulfur.dioxide + chlorides:density + free.sulfur.dioxide:total.sulfur.dioxide + free.sulfur.dioxide:density + free.sulfur.dioxide:sulphates + free.sulfur.dioxide:alcohol + total.sulfur.dioxide:density + total.sulfur.dioxide:alcohol + density:sulphates + pH:sulphates, data = WineQuality) multigraphsummary(model2) print(summary(model2)) cat("Model 2 AIC =", AIC(model2, k=2), "\n") # best main factor + two interaction model (^3) model3 = lm(quality ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + total.sulfur.dioxide + pH + sulphates + alcohol + fixed.acidity:volatile.acidity + fixed.acidity:citric.acid + fixed.acidity:residual.sugar + fixed.acidity:chlorides + fixed.acidity:free.sulfur.dioxide + fixed.acidity:total.sulfur.dioxide + fixed.acidity:density + fixed.acidity:pH + fixed.acidity:alcohol + volatile.acidity:citric.acid + volatile.acidity:residual.sugar + volatile.acidity:chlorides + volatile.acidity:free.sulfur.dioxide + volatile.acidity:total.sulfur.dioxide + volatile.acidity:density + volatile.acidity:sulphates + volatile.acidity:alcohol + citric.acid:residual.sugar + citric.acid:chlorides + citric.acid:free.sulfur.dioxide + citric.acid:total.sulfur.dioxide + citric.acid:density + citric.acid:pH + citric.acid:sulphates + citric.acid:alcohol + residual.sugar:chlorides + residual.sugar:free.sulfur.dioxide + residual.sugar:total.sulfur.dioxide + residual.sugar:density + residual.sugar:pH + residual.sugar:sulphates + residual.sugar:alcohol + chlorides:free.sulfur.dioxide + chlorides:density + chlorides:pH + chlorides:sulphates + chlorides:alcohol + free.sulfur.dioxide:total.sulfur.dioxide + free.sulfur.dioxide:density + free.sulfur.dioxide:pH + free.sulfur.dioxide:sulphates + free.sulfur.dioxide:alcohol + total.sulfur.dioxide:density + total.sulfur.dioxide:pH + total.sulfur.dioxide:sulphates + total.sulfur.dioxide:alcohol + density:pH + density:sulphates + density:alcohol + pH:sulphates + fixed.acidity:citric.acid:density + fixed.acidity:citric.acid:pH + fixed.acidity:residual.sugar:free.sulfur.dioxide + fixed.acidity:residual.sugar:total.sulfur.dioxide + fixed.acidity:residual.sugar:density + fixed.acidity:chlorides:density + fixed.acidity:chlorides:sulphates + fixed.acidity:chlorides:alcohol + fixed.acidity:free.sulfur.dioxide:total.sulfur.dioxide + fixed.acidity:free.sulfur.dioxide:density + fixed.acidity:free.sulfur.dioxide:pH + fixed.acidity:free.sulfur.dioxide:sulphates + fixed.acidity:total.sulfur.dioxide:density + fixed.acidity:total.sulfur.dioxide:alcohol + fixed.acidity:density:pH + fixed.acidity:density:alcohol + volatile.acidity:citric.acid:residual.sugar + volatile.acidity:citric.acid:chlorides + volatile.acidity:citric.acid:total.sulfur.dioxide + volatile.acidity:citric.acid:alcohol + volatile.acidity:residual.sugar:chlorides + volatile.acidity:residual.sugar:alcohol + volatile.acidity:chlorides:pH + volatile.acidity:chlorides:alcohol + volatile.acidity:free.sulfur.dioxide:pH + volatile.acidity:free.sulfur.dioxide:alcohol + volatile.acidity:density:sulphates + volatile.acidity:pH:sulphates + volatile.acidity:sulphates:alcohol + citric.acid:residual.sugar:chlorides + citric.acid:residual.sugar:free.sulfur.dioxide + citric.acid:residual.sugar:total.sulfur.dioxide + citric.acid:chlorides:free.sulfur.dioxide + citric.acid:free.sulfur.dioxide:density + citric.acid:total.sulfur.dioxide:density + citric.acid:total.sulfur.dioxide:pH + citric.acid:total.sulfur.dioxide:alcohol + citric.acid:density:pH + citric.acid:density:alcohol + citric.acid:pH:sulphates + citric.acid:pH:alcohol + residual.sugar:chlorides:density + residual.sugar:chlorides:alcohol + residual.sugar:free.sulfur.dioxide:density + residual.sugar:free.sulfur.dioxide:alcohol + residual.sugar:total.sulfur.dioxide:density + residual.sugar:total.sulfur.dioxide:alcohol + residual.sugar:density:pH + residual.sugar:sulphates:alcohol + chlorides:free.sulfur.dioxide:total.sulfur.dioxide + chlorides:density:pH + chlorides:density:sulphates + chlorides:density:alcohol + chlorides:pH:sulphates + free.sulfur.dioxide:total.sulfur.dioxide:density + free.sulfur.dioxide:density:pH + free.sulfur.dioxide:density:sulphates + free.sulfur.dioxide:density:alcohol + free.sulfur.dioxide:pH:alcohol + free.sulfur.dioxide:sulphates:alcohol + total.sulfur.dioxide:density:pH + total.sulfur.dioxide:density:sulphates + total.sulfur.dioxide:density:alcohol + total.sulfur.dioxide:pH:sulphates + total.sulfur.dioxide:pH:alcohol + total.sulfur.dioxide:sulphates:alcohol + density:pH:sulphates + fixed.acidity:residual.sugar:chlorides, data = WineQuality) multigraphsummary(model3) print(summary(model3)) cat("Model 3 AIC =", AIC(model3, k=2), "\n") # best main factor + two interaction model (^3) model4 = lm(quality ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + total.sulfur.dioxide + density + pH + sulphates + alcohol + fixed.acidity:volatile.acidity + fixed.acidity:citric.acid + fixed.acidity:residual.sugar + fixed.acidity:chlorides + fixed.acidity:free.sulfur.dioxide + fixed.acidity:total.sulfur.dioxide + fixed.acidity:density + fixed.acidity:pH + fixed.acidity:sulphates + fixed.acidity:alcohol + volatile.acidity:citric.acid + volatile.acidity:residual.sugar + volatile.acidity:chlorides + volatile.acidity:free.sulfur.dioxide + volatile.acidity:total.sulfur.dioxide + volatile.acidity:density + volatile.acidity:pH + volatile.acidity:sulphates + volatile.acidity:alcohol + citric.acid:residual.sugar + citric.acid:chlorides + citric.acid:free.sulfur.dioxide + citric.acid:total.sulfur.dioxide + citric.acid:density + citric.acid:pH + citric.acid:sulphates + citric.acid:alcohol + residual.sugar:chlorides + residual.sugar:free.sulfur.dioxide + residual.sugar:total.sulfur.dioxide + residual.sugar:density + residual.sugar:pH + residual.sugar:sulphates + residual.sugar:alcohol + chlorides:free.sulfur.dioxide + chlorides:total.sulfur.dioxide + chlorides:density + chlorides:pH + chlorides:sulphates + chlorides:alcohol + free.sulfur.dioxide:total.sulfur.dioxide + free.sulfur.dioxide:density + free.sulfur.dioxide:pH + free.sulfur.dioxide:sulphates + free.sulfur.dioxide:alcohol + total.sulfur.dioxide:density + total.sulfur.dioxide:pH + total.sulfur.dioxide:sulphates + total.sulfur.dioxide:alcohol + density:pH + density:sulphates + density:alcohol + pH:sulphates + pH:alcohol + sulphates:alcohol + fixed.acidity:volatile.acidity:citric.acid + fixed.acidity:volatile.acidity:residual.sugar + fixed.acidity:volatile.acidity:chlorides + fixed.acidity:volatile.acidity:free.sulfur.dioxide + fixed.acidity:volatile.acidity:total.sulfur.dioxide + fixed.acidity:volatile.acidity:density + fixed.acidity:volatile.acidity:pH + fixed.acidity:volatile.acidity:sulphates + fixed.acidity:volatile.acidity:alcohol + fixed.acidity:citric.acid:residual.sugar + fixed.acidity:citric.acid:chlorides + fixed.acidity:citric.acid:free.sulfur.dioxide + fixed.acidity:citric.acid:total.sulfur.dioxide + fixed.acidity:citric.acid:density + fixed.acidity:citric.acid:pH + fixed.acidity:citric.acid:sulphates + fixed.acidity:citric.acid:alcohol + fixed.acidity:residual.sugar:chlorides + fixed.acidity:residual.sugar:free.sulfur.dioxide + fixed.acidity:residual.sugar:total.sulfur.dioxide + fixed.acidity:residual.sugar:density + fixed.acidity:residual.sugar:pH + fixed.acidity:residual.sugar:sulphates + fixed.acidity:residual.sugar:alcohol + fixed.acidity:chlorides:free.sulfur.dioxide + fixed.acidity:chlorides:total.sulfur.dioxide + fixed.acidity:chlorides:density + fixed.acidity:chlorides:pH + fixed.acidity:chlorides:sulphates + fixed.acidity:chlorides:alcohol + fixed.acidity:free.sulfur.dioxide:total.sulfur.dioxide + fixed.acidity:free.sulfur.dioxide:density + fixed.acidity:free.sulfur.dioxide:pH + fixed.acidity:free.sulfur.dioxide:sulphates + fixed.acidity:free.sulfur.dioxide:alcohol + fixed.acidity:total.sulfur.dioxide:density + fixed.acidity:total.sulfur.dioxide:pH + fixed.acidity:total.sulfur.dioxide:sulphates + fixed.acidity:total.sulfur.dioxide:alcohol + fixed.acidity:density:pH + fixed.acidity:density:sulphates + fixed.acidity:density:alcohol + fixed.acidity:pH:sulphates + fixed.acidity:pH:alcohol + fixed.acidity:sulphates:alcohol + volatile.acidity:citric.acid:residual.sugar + volatile.acidity:citric.acid:chlorides + volatile.acidity:citric.acid:free.sulfur.dioxide + volatile.acidity:citric.acid:total.sulfur.dioxide + volatile.acidity:citric.acid:density + volatile.acidity:citric.acid:pH + volatile.acidity:citric.acid:sulphates + volatile.acidity:citric.acid:alcohol + volatile.acidity:residual.sugar:chlorides + volatile.acidity:residual.sugar:free.sulfur.dioxide + volatile.acidity:residual.sugar:total.sulfur.dioxide + volatile.acidity:residual.sugar:density + volatile.acidity:residual.sugar:pH + volatile.acidity:residual.sugar:sulphates + volatile.acidity:residual.sugar:alcohol + volatile.acidity:chlorides:free.sulfur.dioxide + volatile.acidity:chlorides:total.sulfur.dioxide + volatile.acidity:chlorides:density + volatile.acidity:chlorides:pH + volatile.acidity:chlorides:sulphates + volatile.acidity:chlorides:alcohol + volatile.acidity:free.sulfur.dioxide:density + volatile.acidity:free.sulfur.dioxide:pH + volatile.acidity:free.sulfur.dioxide:sulphates + volatile.acidity:free.sulfur.dioxide:alcohol + volatile.acidity:total.sulfur.dioxide:density + volatile.acidity:total.sulfur.dioxide:pH + volatile.acidity:total.sulfur.dioxide:sulphates + volatile.acidity:total.sulfur.dioxide:alcohol + volatile.acidity:density:pH + volatile.acidity:density:sulphates + volatile.acidity:density:alcohol + volatile.acidity:pH:sulphates + volatile.acidity:pH:alcohol + volatile.acidity:sulphates:alcohol + citric.acid:residual.sugar:chlorides + citric.acid:residual.sugar:free.sulfur.dioxide + citric.acid:residual.sugar:total.sulfur.dioxide + citric.acid:residual.sugar:density + citric.acid:residual.sugar:pH + citric.acid:residual.sugar:sulphates + citric.acid:residual.sugar:alcohol + citric.acid:chlorides:free.sulfur.dioxide + citric.acid:chlorides:total.sulfur.dioxide + citric.acid:chlorides:density + citric.acid:chlorides:pH + citric.acid:chlorides:sulphates + citric.acid:chlorides:alcohol + citric.acid:free.sulfur.dioxide:total.sulfur.dioxide + citric.acid:free.sulfur.dioxide:density + citric.acid:free.sulfur.dioxide:pH + citric.acid:free.sulfur.dioxide:sulphates + citric.acid:free.sulfur.dioxide:alcohol + citric.acid:total.sulfur.dioxide:density + citric.acid:total.sulfur.dioxide:pH + citric.acid:total.sulfur.dioxide:sulphates + citric.acid:total.sulfur.dioxide:alcohol + citric.acid:density:pH + citric.acid:density:sulphates + citric.acid:density:alcohol + citric.acid:pH:sulphates + citric.acid:pH:alcohol + citric.acid:sulphates:alcohol + residual.sugar:chlorides:free.sulfur.dioxide + residual.sugar:chlorides:total.sulfur.dioxide + residual.sugar:chlorides:density + residual.sugar:chlorides:pH + residual.sugar:chlorides:sulphates + residual.sugar:chlorides:alcohol + residual.sugar:free.sulfur.dioxide:total.sulfur.dioxide + residual.sugar:free.sulfur.dioxide:density + residual.sugar:free.sulfur.dioxide:pH + residual.sugar:free.sulfur.dioxide:sulphates + residual.sugar:free.sulfur.dioxide:alcohol + residual.sugar:total.sulfur.dioxide:density + residual.sugar:total.sulfur.dioxide:pH + residual.sugar:total.sulfur.dioxide:sulphates + residual.sugar:total.sulfur.dioxide:alcohol + residual.sugar:density:pH + residual.sugar:density:sulphates + residual.sugar:density:alcohol + residual.sugar:pH:sulphates + residual.sugar:pH:alcohol + residual.sugar:sulphates:alcohol + chlorides:free.sulfur.dioxide:total.sulfur.dioxide + chlorides:free.sulfur.dioxide:density + chlorides:free.sulfur.dioxide:pH + chlorides:free.sulfur.dioxide:sulphates + chlorides:free.sulfur.dioxide:alcohol + chlorides:total.sulfur.dioxide:density + chlorides:total.sulfur.dioxide:pH + chlorides:total.sulfur.dioxide:sulphates + chlorides:total.sulfur.dioxide:alcohol + chlorides:density:pH + chlorides:density:sulphates + chlorides:density:alcohol + chlorides:pH:sulphates + chlorides:pH:alcohol + chlorides:sulphates:alcohol + free.sulfur.dioxide:total.sulfur.dioxide:density + free.sulfur.dioxide:total.sulfur.dioxide:pH + free.sulfur.dioxide:total.sulfur.dioxide:sulphates + free.sulfur.dioxide:density:pH + free.sulfur.dioxide:density:sulphates + free.sulfur.dioxide:density:alcohol + free.sulfur.dioxide:pH:sulphates + free.sulfur.dioxide:pH:alcohol + free.sulfur.dioxide:sulphates:alcohol + total.sulfur.dioxide:density:pH + total.sulfur.dioxide:density:sulphates + total.sulfur.dioxide:density:alcohol + total.sulfur.dioxide:pH:sulphates + total.sulfur.dioxide:pH:alcohol + total.sulfur.dioxide:sulphates:alcohol + density:pH:sulphates + density:pH:alcohol + density:sulphates:alcohol + pH:sulphates:alcohol + fixed.acidity:volatile.acidity:citric.acid:chlorides + fixed.acidity:volatile.acidity:citric.acid:free.sulfur.dioxide + fixed.acidity:volatile.acidity:citric.acid:density + fixed.acidity:volatile.acidity:citric.acid:pH + fixed.acidity:volatile.acidity:citric.acid:sulphates + fixed.acidity:volatile.acidity:citric.acid:alcohol + fixed.acidity:volatile.acidity:residual.sugar:chlorides + fixed.acidity:volatile.acidity:residual.sugar:free.sulfur.dioxide + fixed.acidity:volatile.acidity:residual.sugar:total.sulfur.dioxide + fixed.acidity:volatile.acidity:residual.sugar:density + fixed.acidity:volatile.acidity:residual.sugar:sulphates + fixed.acidity:volatile.acidity:residual.sugar:alcohol + fixed.acidity:volatile.acidity:chlorides:total.sulfur.dioxide + fixed.acidity:volatile.acidity:chlorides:alcohol + fixed.acidity:volatile.acidity:free.sulfur.dioxide:density + fixed.acidity:volatile.acidity:free.sulfur.dioxide:pH + fixed.acidity:volatile.acidity:total.sulfur.dioxide:alcohol + fixed.acidity:volatile.acidity:density:pH + fixed.acidity:volatile.acidity:density:alcohol + fixed.acidity:volatile.acidity:pH:alcohol + fixed.acidity:volatile.acidity:sulphates:alcohol + fixed.acidity:citric.acid:residual.sugar:chlorides + fixed.acidity:citric.acid:residual.sugar:free.sulfur.dioxide + fixed.acidity:citric.acid:chlorides:density + fixed.acidity:citric.acid:chlorides:pH + fixed.acidity:citric.acid:chlorides:sulphates + fixed.acidity:citric.acid:chlorides:alcohol + fixed.acidity:citric.acid:free.sulfur.dioxide:total.sulfur.dioxide + fixed.acidity:citric.acid:total.sulfur.dioxide:density + fixed.acidity:citric.acid:total.sulfur.dioxide:sulphates + fixed.acidity:citric.acid:density:sulphates + fixed.acidity:citric.acid:density:alcohol + fixed.acidity:citric.acid:sulphates:alcohol + fixed.acidity:residual.sugar:chlorides:free.sulfur.dioxide + fixed.acidity:residual.sugar:chlorides:density + fixed.acidity:residual.sugar:free.sulfur.dioxide:pH + fixed.acidity:residual.sugar:free.sulfur.dioxide:sulphates + fixed.acidity:residual.sugar:free.sulfur.dioxide:alcohol + fixed.acidity:residual.sugar:density:sulphates + fixed.acidity:residual.sugar:sulphates:alcohol + fixed.acidity:chlorides:free.sulfur.dioxide:sulphates + fixed.acidity:chlorides:pH:alcohol + fixed.acidity:free.sulfur.dioxide:total.sulfur.dioxide:sulphates + fixed.acidity:free.sulfur.dioxide:density:pH + fixed.acidity:free.sulfur.dioxide:density:sulphates + fixed.acidity:free.sulfur.dioxide:pH:alcohol + fixed.acidity:free.sulfur.dioxide:sulphates:alcohol + fixed.acidity:total.sulfur.dioxide:density:pH + fixed.acidity:total.sulfur.dioxide:density:alcohol + fixed.acidity:total.sulfur.dioxide:pH:sulphates + fixed.acidity:density:pH:sulphates + fixed.acidity:density:pH:alcohol + fixed.acidity:pH:sulphates:alcohol + volatile.acidity:citric.acid:residual.sugar:chlorides + volatile.acidity:citric.acid:residual.sugar:free.sulfur.dioxide + volatile.acidity:citric.acid:residual.sugar:density + volatile.acidity:citric.acid:residual.sugar:pH + volatile.acidity:citric.acid:residual.sugar:alcohol + volatile.acidity:citric.acid:chlorides:density + volatile.acidity:citric.acid:chlorides:pH + volatile.acidity:citric.acid:chlorides:alcohol + volatile.acidity:citric.acid:total.sulfur.dioxide:density + volatile.acidity:citric.acid:total.sulfur.dioxide:pH + volatile.acidity:citric.acid:density:pH + volatile.acidity:citric.acid:density:sulphates + volatile.acidity:citric.acid:pH:sulphates + volatile.acidity:citric.acid:pH:alcohol + volatile.acidity:citric.acid:sulphates:alcohol + volatile.acidity:residual.sugar:chlorides:density + volatile.acidity:residual.sugar:free.sulfur.dioxide:density + volatile.acidity:residual.sugar:free.sulfur.dioxide:pH + volatile.acidity:residual.sugar:free.sulfur.dioxide:sulphates + volatile.acidity:residual.sugar:free.sulfur.dioxide:alcohol + volatile.acidity:residual.sugar:total.sulfur.dioxide:density + volatile.acidity:residual.sugar:total.sulfur.dioxide:sulphates + volatile.acidity:residual.sugar:total.sulfur.dioxide:alcohol + volatile.acidity:residual.sugar:density:pH + volatile.acidity:residual.sugar:density:alcohol + volatile.acidity:residual.sugar:pH:sulphates + volatile.acidity:chlorides:free.sulfur.dioxide:density + volatile.acidity:chlorides:free.sulfur.dioxide:sulphates + volatile.acidity:chlorides:total.sulfur.dioxide:density + volatile.acidity:chlorides:total.sulfur.dioxide:alcohol + volatile.acidity:chlorides:density:pH + volatile.acidity:chlorides:pH:alcohol + volatile.acidity:free.sulfur.dioxide:density:alcohol + volatile.acidity:free.sulfur.dioxide:pH:sulphates + volatile.acidity:total.sulfur.dioxide:density:sulphates + volatile.acidity:total.sulfur.dioxide:density:alcohol + volatile.acidity:total.sulfur.dioxide:pH:sulphates + volatile.acidity:total.sulfur.dioxide:sulphates:alcohol + volatile.acidity:density:pH:alcohol + citric.acid:residual.sugar:chlorides:free.sulfur.dioxide + citric.acid:residual.sugar:chlorides:total.sulfur.dioxide + citric.acid:residual.sugar:chlorides:density + citric.acid:residual.sugar:chlorides:alcohol + citric.acid:residual.sugar:free.sulfur.dioxide:pH + citric.acid:residual.sugar:free.sulfur.dioxide:sulphates + citric.acid:residual.sugar:total.sulfur.dioxide:pH + citric.acid:residual.sugar:total.sulfur.dioxide:sulphates + citric.acid:residual.sugar:density:sulphates + citric.acid:residual.sugar:sulphates:alcohol + citric.acid:chlorides:free.sulfur.dioxide:pH + citric.acid:chlorides:total.sulfur.dioxide:density + citric.acid:chlorides:total.sulfur.dioxide:sulphates + citric.acid:chlorides:total.sulfur.dioxide:alcohol + citric.acid:chlorides:density:pH + citric.acid:chlorides:density:sulphates + citric.acid:chlorides:density:alcohol + citric.acid:chlorides:pH:sulphates + citric.acid:chlorides:pH:alcohol + citric.acid:free.sulfur.dioxide:total.sulfur.dioxide:pH + citric.acid:free.sulfur.dioxide:density:pH + citric.acid:free.sulfur.dioxide:density:sulphates + citric.acid:free.sulfur.dioxide:pH:sulphates + citric.acid:free.sulfur.dioxide:pH:alcohol + citric.acid:free.sulfur.dioxide:sulphates:alcohol + citric.acid:total.sulfur.dioxide:density:pH + citric.acid:total.sulfur.dioxide:density:sulphates + citric.acid:total.sulfur.dioxide:pH:alcohol + citric.acid:density:pH:sulphates + citric.acid:density:sulphates:alcohol + residual.sugar:chlorides:free.sulfur.dioxide:density + residual.sugar:chlorides:free.sulfur.dioxide:pH + residual.sugar:chlorides:free.sulfur.dioxide:alcohol + residual.sugar:chlorides:total.sulfur.dioxide:sulphates + residual.sugar:chlorides:density:alcohol + residual.sugar:chlorides:pH:sulphates + residual.sugar:free.sulfur.dioxide:total.sulfur.dioxide:pH + residual.sugar:free.sulfur.dioxide:total.sulfur.dioxide:sulphates + residual.sugar:free.sulfur.dioxide:density:sulphates + residual.sugar:free.sulfur.dioxide:density:alcohol + residual.sugar:free.sulfur.dioxide:pH:sulphates + residual.sugar:free.sulfur.dioxide:pH:alcohol + residual.sugar:free.sulfur.dioxide:sulphates:alcohol + residual.sugar:total.sulfur.dioxide:density:pH + residual.sugar:total.sulfur.dioxide:density:sulphates + residual.sugar:total.sulfur.dioxide:pH:sulphates + residual.sugar:total.sulfur.dioxide:pH:alcohol + residual.sugar:total.sulfur.dioxide:sulphates:alcohol + residual.sugar:density:pH:sulphates + residual.sugar:density:pH:alcohol + residual.sugar:density:sulphates:alcohol + residual.sugar:pH:sulphates:alcohol + chlorides:free.sulfur.dioxide:total.sulfur.dioxide:density + chlorides:free.sulfur.dioxide:density:pH + chlorides:free.sulfur.dioxide:density:sulphates + chlorides:free.sulfur.dioxide:density:alcohol + chlorides:free.sulfur.dioxide:pH:sulphates + chlorides:free.sulfur.dioxide:pH:alcohol + chlorides:free.sulfur.dioxide:sulphates:alcohol + chlorides:total.sulfur.dioxide:density:pH + chlorides:total.sulfur.dioxide:pH:sulphates + chlorides:total.sulfur.dioxide:pH:alcohol + chlorides:density:sulphates:alcohol + total.sulfur.dioxide:density:pH:alcohol + volatile.acidity:chlorides:sulphates:alcohol + fixed.acidity:residual.sugar:chlorides:pH + volatile.acidity:citric.acid:free.sulfur.dioxide:alcohol + volatile.acidity:pH:sulphates:alcohol + residual.sugar:chlorides:density:pH, data = WineQuality) multigraphsummary(model4) print(summary(model4)) cat("Model 3 AIC =", AIC(model4, k=2), "\n") cat("Model 1 AIC =", AIC(model1, k=2), "\n") cat("Model 2 AIC =", AIC(model2, k=2), "\n") cat("Model 3 AIC =", AIC(model3, k=2), "\n") cat("Model 4 AIC =", AIC(model4, k=2), "\n") # stepwise search for best regression model library(MASS) fit <- lm(quality~.^4, data=WineQuality) n = length(WineQuality$quality) step <- stepAIC(fit, direction="both", k=2) step$anova # display results