Statistics

Statistics 108, Final Project

Due : December 2nd, In Class

The final project is an opportunity to put all you have learned in STA 108 in action. Your task is to conduct a start-to-finish statistical analysis of a real data set. I will guide you through how to analysis a real data set, diabetes data. This data consist of 16 variables on 403 subjects who were interviewed in a study to understand the prevalence of obesity, diabetes, and other cardiovascular risk factors in central Virginia for African Americans. We will consider building regression models with glyhb as the response variable as Glycosolated Hemoglobin > 70 is often taken as a positive diagnostics of diabetes. The data is on smartsite under Resources/Homework/diabetesfull.txt. Please attach your R codes and plots.

Processing and exploration of the data.

1. Drop all the cases having missing value(s). (There are other better ways to deal with missing data and here we take the brutal easy way. You would end up with 366 cases.) You may use the code:

mydata = read.table("diabetesfull.txt", header=T)

index.na = apply(is.na(mydata), 1, any) ## identify cases having NA values.

newdata = mydata[index.na==FALSE,] ## drop cases with NA entries.

any(is.na(newdata)) ## this should return FALSE -- no NA in data.s

dim(newdata) ## this should return 366 16: 366 cases, 16 variables.

2. Which of the variables are quantitative variables? Which are qualitative variables? Draw histogram for each quantitative variable and comment on its distribution. Draw pie chart for each qualitative variable and comment on how its classes are distributed.

3. Draw scatterplot matrix and obtain the pairwise correlation matrix for all quantitative vari- ables in the data. Comment on their relationships.

Diagnostic and Remedy.

4. Regress glybh on all predictor variables (Model 1). Draw the diagnostic plots of the model and comment.

5. You want to check whether any transformation on the response variable is needed. You use the function ‘boxcox’ to help you make the decision. State the transformation you decide to use. In the following, we denote the transformed response variable to be glyhb∗. Regress glyhb∗ on all predictor variables (Model 2). Draw the diagnostic plots of this model and comment. Apply boxcox again on Model 2; what do you find?

Split data for validation later on.

6. Randomly split data into two equal halves: a training data set and a validation data set. You may use the code.

1

set.seed(10) ## set seed for random number generator

## so everyone gets the same split of the data.

N=nrow(newdata) ## number of cases in the data (366)

index=sample(1: N, size=N/2, replace=FALSE)

## randomly sample 183 cases to form the training data.

data.t=newdata[index,] ## get the training data set.

data.v=newdata[-index,] ## the remaining 183 cases form the validation set.

Selection of first-order effects. We now consider subsets selection from the pool of all first-order effects of the 15 predictors. glyhb* is used as the response variable for the following problems.

7. Fit a model with all first-order effects (Model 3). How many regression coefficients are there in this model? What is the MSE from this model?

Hint: you may use the pseudo-code:

model3 = lm(f(glyhb)~., data=data.t)

where f(glyhb) is the transformation you decided to use for the response variable (e.g., if you decide to use the ‘log’ transformation, you would put log(glyhb) there), and data.t denotes the training data.

8. Consider best subsets selection using the R function regsubsets() from the leaps library with Model 3 as the full model. Return the top 1 best subset of all subset sizes (i.e., number of X variables) up to 16 (because frame has 3 levels). Get SSEp,R

2 p, R

2 a,p, Cp, AICp,BICp for

each of these models, as well as the none-model (the model with only an intercept). Identify the best model according to each criterion. For the best model according to Cp criterion, what do you observe about its Cp value? Do you have a possible explanation for it?

Denote the best models according to AIC, BIC, and adjusted R2 be Model 3.1, Model 3.2, Model 3.3, respectively. (It is possible that some of the three models are the same.)

Hint: you may use the pseudo-code:

best = regsubsets(f(glyhb)~., data=data.t, nbest=1, nvmax=16)

Selection of first- and second- order effects. We now consider subsets selection from the pool of first-order effects as well as 2-way interaction effects of the 15 predictors.

9. Fit a model with all first-order and 2-way interaction effects (Model 4). How many regression coefficients are there in this model? What is the MSE from this model? Do you have any concern about the fitting of this model and why?

Hint: you may use the pseudo-code:

model4 = lm(f(glyhb)~.^2, data=data.t)

10. Apply the forward stepwise procedure using R function step() (or stepAIC()), starting from the none-model and using the AICp criterion. What is the model being selected? Denote this model by Model fs1. Compare its AIC value with that of Model3.1. What do you find?

2

11. Apply the forward stepwise procedure using R function step() (or stepAIC()), starting from the full model (Model 3) and using the AICp criterion. What is the model being se- lected? Denote this model by Model fs2. Compare its AIC value with that of Model fs1. What do you find?

12. Compare the BIC values of Model fs1 and Model fs2. What do you find? Do AIC and BIC choose the same model among these two models or not? Denote the model selected by AIC among the two models by Model 4.1 and that selected by BIC be Model 4.2. (It is possible that Model 4.1 and Model 4.2 are the same model.) (In real data analysis, we may run ‘step()’ function with more starting models and choose among the resulting models the best ones according to their AIC and BIC values.)

Notes: If you are curious, you can have a try of best subsets selection using the R pseudo-code

best2 = regsubsets(f(glyhb)~.^2, data=data.t, nbest=1, really.big=T)

However, be careful, you may have to force your way out of R session due to the slowness of this procedure! So save all that you want to save before you try this.

Model validation. We now consider validation of the models (Model 3.1, Model 3.2, Model 3.3, Model 4.1, Model4.2) you selected in the previous studies.

13. Internal validation. We use PRESS for this purpose. Calculate PRESS for each of these models. Comment.

14. External validation using the validation set. For each of these models (Model 3.1, Model 3.2, Model 3.3, Model 4.1, Model4.2), calculate the mean squared prediction error (MSPR), i.e., you use the model to predict the 183 observations in the validation set and calculate the averaged squared prediction error. How do these MSPRs compare with the respective PRSSE/n (here n is the sample size of the training data, i.e., 183). Which model has the smallest MSPR?

15. Based on both internal and external validation, which model you would choose as the final model? Fit the final model using the entire data set (training and validation combined) (Model 5). Write down the fitted regression function and report the R summary() and anova() output.

Model Diagnostics.

16. Draw the diagnostic plots for the final model, Model 5. Comment.

17. Which case has the highest leverage?

18. Calculate studentized deleted residuals. Test whether there are outliers with respect to the response variable using the t-test with Bonferroni correction. (Significance level α = 0.05)

19. Calculate Cook’s distance. Which case has the largest Cook’s distance?

20. Is multicollinearity a big issue here? i.e., what is the variance inflation factor for each of the predictor variables in the model?

3