Hwk #2: Classification methods and Penalization
For this homework we will use NHANES data that exists in a package for R.
NHANES consists of survey data collected by the US National Center for Health Statistics (NCHS) which has conducted a series of health and nutrition surveys since the early 1960’s. Since 1999 approximately 5,000 individuals of all ages are interviewed in their homes every year and complete the health examination component of the survey. The health examination is conducted in a mobile examination center (MEC).
Note that there is the following warning on the NHANES website: “For NHANES datasets, the use of sampling weights and sample design variables is recommended for all analyses because the sample design is a clustered design and incorporates differential probabilities of selection. If you fail to account for the sampling parameters, you may obtain biased estimates and overstate significance levels.”
For this homework, please ignore this warning and just apply our analyses to the data as if they were randomly sampled! We will be using the data called NHANESraw.
For questions that ask for your comments, it suffices to answer with one or two sentences in each case.
Data Preparation
- Install the package
NHANESinto R, load theNHANESpackage, and then run the commanddata(NHANES)which will load the NHANES data. Type?NHANESand read about the dataset.
# install.packages("NHANES")
library(NHANES)
data(NHANES)
?NHANES- Make an object
nhanesthat is a subset versionNHANESrawthat does not include any missing data forDiabetes,BPSysAve,BPDiaAve, orAge.
library(tidyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnhanes <- NHANESraw %>% filter(!is.na(Diabetes) & !is.na(BPSysAve) & !is.na(BPDiaAve) & !is.na(Age))- (1 point) Further subset the data such the observations with
BPDiaAveequal to zero are removed.
nhanesnozeros <- nhanes %>% filter(BPDiaAve != 0)- (1 point) Make an object
nhanes09that is a subset ofnhanesto only the 2009_10 data. This will be your training dataset. Also make an objectnhanes11that is a subset ofnhanesto only the 2011_12 data. This will be your test dataset.
nhanes09 <- nhanesnozeros %>% filter(SurveyYr == "2009_10")
nhanes11 <- nhanesnozeros %>% filter(SurveyYr == "2011_12")
#rm(nhanesnozeros) # to save memory
#rm(nhanes) # to save memoryLogistic regression
- (2 point) Fit a logistic regression model (call it
glm1) using thenhanes09dataset. UseDiabetesas the outcome and averaged systolic blood pressure (BPSysAve) as a single predictor. Use the summary command to examine the fitted model. Generate the 95% confidence intervals for theBPSysAvecoefficient.
glm1 <- glm(Diabetes ~ BPSysAve, family="binomial", data=nhanes09)
summary(glm1)##
## Call:
## glm(formula = Diabetes ~ BPSysAve, family = "binomial", data = nhanes09)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5554 -0.4892 -0.4149 -0.3510 2.5957
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.645247 0.224502 -25.15 <2e-16 ***
## BPSysAve 0.028892 0.001761 16.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5245.1 on 7724 degrees of freedom
## Residual deviance: 4983.9 on 7723 degrees of freedom
## AIC: 4987.9
##
## Number of Fisher Scoring iterations: 5confint(glm1)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -6.08759789 -5.20726141
## BPSysAve 0.02544521 0.03235066- (1 point) Generate the estimate and 95% confidence interval for the odds-ratio associated with BPSysAve. Summarize the result.
exp(coefficients(glm1)["BPSysAve"])## BPSysAve
## 1.029313exp(confint(glm1)["BPSysAve",])## Waiting for profiling to be done...## 2.5 % 97.5 %
## 1.025772 1.032880# For each unit increase in Systolic Blood Pressure we see about a 2.93% increase in odds of diabetes, 95% CI 2.58% to 3.29%.- (1 point) Predict the probabilities of diabetes associated with each of the training observations of
BPSysAve. Make a vector of predictions for diabetes based on whether the predictions are above or below 0.5.
pred09 <- predict(glm1, type="response") ## Generate predicted probabilities
pred09binary <- ifelse(pred09 > 0.5, "Yes", "No") ## Threshold at 0.5- (1 point) Generate a confusion matrix that shows the number of false positives, false negatives, true positives, and true negatives in the training data. The rows should correspond to the true diabetes status and the columns should correspond to the predicted values.
table(true_disease=nhanes09$Diabetes, predictions=pred09binary)## predictions
## true_disease No Yes
## No 6885 16
## Yes 817 7- (1 point) Find the proportion of correctly classified observations in the training data.
mean(nhanes09$Diabetes == pred09binary)## [1] 0.8921683- (2 points) Now repeat questions 7 to 9 but for predicting the test dataset.
pred11 <- predict(glm1,type="response",newdata=nhanes11)
pred11binary <- ifelse(pred11 > 0.5,"Yes","No")
xt11 <- table(true_disease=nhanes11$Diabetes, predictions=pred11binary)
print(xt11)## predictions
## true_disease No Yes
## No 6195 13
## Yes 757 4print(paste("Test accuracy", mean(nhanes11$Diabetes == pred11binary)))## [1] "Test accuracy 0.889510690199455"- (1 point) Comment on the difference in results between the training and test prediction tables and classification accuracies.
# The results are actually very similar for the training and test datasets. This is not totally surprising given that nearly all predicted probabilities are below 0.5.- (1 point) Manually calculate the sensitivity and specificity estimates for the test dataset based on the 0.5 threshold.
sens <- xt11[2,2]/sum(xt11[2, ])
spec <- xt11[1,1]/sum(xt11[1, ])
print(paste("Sensitivity", sens))## [1] "Sensitivity 0.00525624178712221"print(paste("Specificity", spec))## [1] "Specificity 0.997905927835051"- (2 points) Generate an ROC curve using the test data. What is the AUC and its 95% confidence interval?
library(ROCit)
ROCit_obj <- rocit(score=pred11, class=nhanes11$Diabetes)
plot(ROCit_obj)
ciAUC(ROCit_obj)##
## estimated AUC : 0.718795191995069
## AUC estimation method : empirical
##
## CI of AUC
## confidence level = 95%
## lower = 0.697482998977986 upper = 0.740107385012152- What value can you use to threshold the predicted probability to achieve a sensitivity of at least 0.6 and a specificity of at least 0.3?
# The following cutoffs can all be used to achieve the desired sensitivity and specificity levels
thresholds <- ROCit_obj$Cutoff[(ROCit_obj$TPR > 0.6) & (ROCit_obj$FPR < 0.3)]
print(paste("Possible threshold:", thresholds[1]))## [1] "Possible threshold: 0.115702890399852"- (2 points) Comment on the results of the analyses for the different thresholds in terms of the tables, classification accuracies, and sensitivity and specificity. Under what circumstances might you prefer each of the thresholds?
# The accuracy is highest for the 0.5 threshold, is a little worse for the 0.2 threshold but quite a bit worse for the 0.1 threshold. However, for the 0.5 and 0.2 thresholds the specificity is really high but the sensitivity is poor (slightly better sensitivity for 0.2). Despite having considerably lower classification accuracy, there is more of a balance between sensitivity and specificity for the 0.1 threshold. This set of results reflects the fact that a large majority of the subjects do not have diabetes so in the absence of strong evidence for diabetes you are more likely to get it right by predicting no diabetes. However, if you want to screen for diabetes you need to get the sensitivity up and so may want to use the 0.1 threshold (albeit at the risk of having many false positives.)- (2 points) Fit a multiple predictor logistic regression (call it
glm2) withDiabetesas outcome and predictors:BPSysAve,BPDiaAve, andAge. Use thesummarycommand to examine the fitted model and determine the estimated coefficients, odds-ratios, and 95% confidence intervals thereof.
glm2 <- glm(Diabetes ~ BPSysAve + BPDiaAve + Age, family="binomial", data=nhanes09)
summary(glm2)##
## Call:
## glm(formula = Diabetes ~ BPSysAve + BPDiaAve + Age, family = "binomial",
## data = nhanes09)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1027 -0.4916 -0.2862 -0.1874 2.9814
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.076129 0.283412 -17.911 <2e-16 ***
## BPSysAve 0.004657 0.002281 2.042 0.0412 *
## BPDiaAve -0.002999 0.003181 -0.943 0.3458
## Age 0.051655 0.002423 21.318 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5245.1 on 7724 degrees of freedom
## Residual deviance: 4401.8 on 7721 degrees of freedom
## AIC: 4409.8
##
## Number of Fisher Scoring iterations: 6# The 95% confidence interval for the model coefficients
confint(glm2)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -5.634646626 -4.523306668
## BPSysAve 0.000171014 0.009114841
## BPDiaAve -0.009209046 0.003264903
## Age 0.046963207 0.056464977# The odds-ratios associated with each predictor
exp(coefficients(glm2))## (Intercept) BPSysAve BPDiaAve Age
## 0.006244034 1.004667430 0.997005198 1.053012909# The 95% confidence interval for the odds-ratios associated with each predictor
exp(confint(glm2))## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 0.003571939 0.01085308
## BPSysAve 1.000171029 1.00915651
## BPDiaAve 0.990833227 1.00327024
## Age 1.048083446 1.05808956- (2 points) Generate an ROC curve for the
glm2model using the test data. What is the AUC and its 95% confidence interval?
pred11_glm2 <- predict(glm2, newdata = nhanes11, type = "response")
ROCit_obj2 <- rocit(score=pred11_glm2, class=nhanes11$Diabetes)
plot(ROCit_obj2)
ciAUC(ROCit_obj2)##
## estimated AUC : 0.812555352256255
## AUC estimation method : empirical
##
## CI of AUC
## confidence level = 95%
## lower = 0.793592417743851 upper = 0.831518286768659- (1 point) What is the maximum sensitivity level you can achieve if we require the specificity to be at least 0.3?
max(ROCit_obj2$TPR[ROCit_obj2$FPR < 0.7])## [1] 0.9868594- (1 point) Would you prefer the single predictor or multiple predictor model if your objective was to maximize classification accuracy, and which threshold level would you choose? Comment on the reason for your choices.
# Probably would threshold based on 0.5 since that is the best threshold for classification accuracy for both models. You might prefer the multi predictor model since it has very slightly higher accuracy on the test set, but it is so marginal that you may opt for the simpler single predictor model, even with the tiny higher performance on the test set.Linear discriminant analysis
- (2 points) Fit a linear discriminant analysis (
lda1) withDiabetesas outcome and predictors ofBPSysAve,BPDiaAve, andAgein the training dataset. Examine the fit by typinglda1.
lda1 <- lda(Diabetes ~ BPSysAve + BPDiaAve + Age, data=nhanes09)
lda1## Call:
## lda(Diabetes ~ BPSysAve + BPDiaAve + Age, data = nhanes09)
##
## Prior probabilities of groups:
## No Yes
## 0.8933333 0.1066667
##
## Group means:
## BPSysAve BPDiaAve Age
## No 116.5382 65.68106 37.78844
## Yes 128.3180 66.86529 60.83374
##
## Coefficients of linear discriminants:
## LD1
## BPSysAve 0.01042485
## BPDiaAve -0.01806242
## Age 0.04469600- (2 points) Generate the confusion matrix for
lda1using the test set. Compute the classification accuracy, sensitivity, and specificity.
predlda1 <- predict(lda1, newdata=nhanes11)
lda1tab <- table(true_disease=nhanes11$Diabetes, predictions=predlda1$class)
print(lda1tab)## predictions
## true_disease No Yes
## No 6158 50
## Yes 738 23print(paste("Accuracy", mean(predlda1$class==nhanes11$Diabetes)))## [1] "Accuracy 0.886927823217104"senslda1 <- lda1tab[2,2]/sum(lda1tab[2, ])
speclda1 <- lda1tab[1,1]/sum(lda1tab[1, ])
print(paste("Sensitivity", senslda1))## [1] "Sensitivity 0.0302233902759527"print(paste("Specificity", speclda1))## [1] "Specificity 0.99194587628866"- How do these measures compare with that of the logistic regression model with these predictors and 0.5 threshold?
# The accuracy is pretty similar 0.8869 (lda) vs. 0.8895 (logistic) between the two approaches. The lda has higher sensitivity relative to logistic at the expense of some specificity.- (3 points) Redo question 21 but with prior probabilities set to 0.5 for diabetes.
predlda2 <- predict(lda1, prior=c(0.5,0.5), newdata=nhanes11)
lda2tab <- table(true_disease=nhanes11$Diabetes, predictions=predlda2$class)
print(lda2tab)## predictions
## true_disease No Yes
## No 4630 1578
## Yes 196 565print(paste("Accuracy", mean(predlda2$class==nhanes11$Diabetes)))## [1] "Accuracy 0.745444109628354"senslda2 <- lda2tab[2,2]/sum(lda2tab[2, ])
speclda2 <- lda2tab[1,1]/sum(lda2tab[1, ])
print(paste("Sensitivity", senslda2))## [1] "Sensitivity 0.742444152431012"print(paste("Specificity", speclda2))## [1] "Specificity 0.745811855670103"- (2 points) Comment on how LDA’s performance changed when we changed the prior probabilities.
# The second model assumedtest patients had an equal probability of being diabetic and non-diabetic a priori. This led to an increase in sensitivity but decreased the specificity. The overall accuracy decreased because the prior distribution did not reflect the proportions seen in the test data.Penalized regression
- Read in the dementia data “dementia2.csv” into a data frame called
dementia_dat. This dataset contains measurements obtained from MRI brain scans and whether or not the patient has dementia. We’ll try to build a prediction model for diagnosing dementia based on these derived measurements. How many observations are in this dataset? How many predictors are in this dataset?
dementia_dat <- read.csv("dementia2.csv")
print(paste("Number of observations", nrow(dementia_dat)))## [1] "Number of observations 660"print(paste("Number of predictors", ncol(dementia_dat) - 1))## [1] "Number of predictors 142"head(dementia_dat)## Dementia MacCohort_kr gm wm csf tiv ageAtScan
## 1 No CONTROL 595.483 521.213 284.383 1401.08 70
## 2 Yes AD 559.454 549.746 292.001 1401.20 72
## 3 Yes AD 687.109 551.160 298.904 1537.17 73
## 4 Yes AD 644.104 586.846 515.496 1746.45 78
## 5 No CONTROL 550.418 422.150 285.142 1257.71 80
## 6 Yes AD 445.386 451.453 356.003 1252.84 78
## X3rd_Ventricle X4th_Ventricle Right_Accumbens_Area Left_Accumbens_Area
## 1 131.49612 166.5999 114.73468 153.24140
## 2 106.60645 140.9746 106.33279 98.76677
## 3 91.13264 194.6998 92.30153 88.31427
## 4 95.40594 194.4836 93.00046 57.81628
## 5 143.28459 181.6107 124.05885 125.37046
## 6 105.41526 140.4759 72.20149 77.43875
## Right_Amygdala Left_Amygdala Brain_Stem Right_Caudate Left_Caudate
## 1 284.2214 312.0922 840.5216 644.7019 745.0837
## 2 253.8847 287.5789 667.2526 511.1797 465.4505
## 3 271.0102 236.9666 951.5518 719.4063 572.9821
## 4 182.7516 147.1373 1035.7474 539.8652 547.1518
## 5 326.3988 379.5273 992.3929 747.4604 777.3925
## 6 140.0999 166.6878 647.0480 583.3195 611.7902
## Right_Cerebellum_Exterior Left_Cerebellum_Exterior
## 1 13712.39 13881.63
## 2 11960.62 11063.44
## 3 12450.55 12168.28
## 4 12214.42 11763.49
## 5 13106.70 13271.89
## 6 11616.47 11339.14
## Right_Cerebellum_White_Matter Left_Cerebellum_White_Matter
## 1 3796.907 3429.957
## 2 3379.344 3056.046
## 3 4534.826 3862.235
## 4 3438.232 2955.012
## 5 3764.693 3343.781
## 6 4219.977 3694.138
## Right_Cerebral_White_Matter Left_Cerebral_White_Matter Cerebrospinal_Fluid
## 1 56278.86 56108.17 170.4707
## 2 58783.69 58325.69 130.3656
## 3 52133.48 51225.53 142.8433
## 4 46805.58 48817.08 113.6674
## 5 50415.45 50328.39 184.8811
## 6 50792.75 50039.58 168.4659
## Right_Hippocampus Left_Hippocampus Right_Inf_Lat_Vent Left_Inf_Lat_Vent
## 1 884.5345 812.9423 52.83013 17.386439
## 2 935.3056 894.4070 53.20210 19.496608
## 3 997.9761 735.0826 53.42905 12.989925
## 4 686.8895 564.1407 32.30037 10.335909
## 5 977.1922 982.9319 55.51559 21.801584
## 6 539.9982 535.7009 19.71797 9.992932
## Right_Lateral_Ventricle Left_Lateral_Ventricle Right_Pallidum Left_Pallidum
## 1 501.5410 615.7695 75.92260 97.81442
## 2 520.4170 537.7467 62.05112 73.44707
## 3 824.7920 785.8833 57.92768 58.18801
## 4 605.7268 797.7107 48.36185 35.96183
## 5 580.2155 701.4868 101.26055 126.58661
## 6 589.7373 790.7854 51.79003 77.73289
## Right_Putamen Left_Putamen Right_Thalamus_Proper Left_Thalamus_Proper
## 1 871.9013 1028.9238 1314.4115 1182.6963
## 2 633.0704 648.5695 1112.7609 1004.9455
## 3 625.7079 573.4950 828.2146 663.0622
## 4 391.4689 309.8539 952.4882 808.7958
## 5 874.3787 922.7681 1326.9221 1259.8541
## 6 547.8218 712.2014 986.8584 812.8650
## Right_Ventral_DC Left_Ventral_DC Right_vessel Left_vessel Optic_Chiasm
## 1 336.7650 387.7478 0.5590566 2.100953 3.117903
## 2 305.9562 340.2070 0.4801367 1.976466 3.151992
## 3 350.2180 334.7909 0.4324495 1.828934 5.370881
## 4 306.0463 290.2501 0.3840831 1.246457 2.522688
## 5 381.3044 435.3156 0.6120818 2.491986 3.759846
## 6 297.7132 303.0005 0.3662870 2.253577 2.658108
## Cerebellar_Vermal_Lobules_I_V Cerebellar_Vermal_Lobules_VI_VII
## 1 1235.934 579.0423
## 2 1004.452 422.5326
## 3 1049.574 484.1025
## 4 1079.287 458.2170
## 5 1227.715 530.5586
## 6 1151.713 460.2984
## Cerebellar_Vermal_Lobules_VIII_X Left_Basal_Forebrain Right_Basal_Forebrain
## 1 651.0357 76.86333 68.23146
## 2 571.5840 64.68308 57.92234
## 3 638.2878 64.79633 59.01797
## 4 612.3502 39.81159 51.16755
## 5 774.3089 87.34242 80.43258
## 6 539.0609 56.17912 43.48826
## Right_ACgG_anterior_cingulate_gyrus Left_ACgG_anterior_cingulate_gyrus
## 1 917.9102 1196.245
## 2 936.2062 1144.866
## 3 994.1086 1160.793
## 4 676.9367 912.061
## 5 789.1565 1047.792
## 6 907.3145 1187.140
## Right_AIns_anterior_insula Left_AIns_anterior_insula
## 1 1044.7071 1083.4238
## 2 1071.9465 982.9163
## 3 1015.6598 917.6926
## 4 794.7559 651.9497
## 5 1265.2635 1130.7485
## 6 950.8980 1095.9243
## Right_AOrG_anterior_orbital_gyrus Left_AOrG_anterior_orbital_gyrus
## 1 362.6913 324.0319
## 2 382.2687 300.9694
## 3 363.4014 310.4434
## 4 300.8172 271.8262
## 5 374.4690 320.1677
## 6 361.3199 345.9985
## Right_AnG_angular_gyrus Left_AnG_angular_gyrus Right_Calc_calcarine_cortex
## 1 2468.202 1819.668 822.5055
## 2 2259.860 1872.589 821.7743
## 3 3000.104 2514.393 907.7510
## 4 2135.954 1710.856 598.0245
## 5 2538.029 2009.256 753.2055
## 6 1424.652 1356.299 631.5251
## Left_Calc_calcarine_cortex Right_CO_central_operculum
## 1 733.3573 863.1964
## 2 824.9935 759.0901
## 3 961.9054 792.6900
## 4 590.0507 583.0569
## 5 784.9695 766.6014
## 6 772.2251 748.5800
## Left_CO_central_operculum Right_Cun_cuneus Left_Cun_cuneus
## 1 889.6048 1047.8334 1072.0937
## 2 775.4344 1035.6598 1104.8444
## 3 816.3310 1238.5967 1257.2403
## 4 583.5119 1130.2484 1134.2509
## 5 744.5943 1135.8323 1058.2455
## 6 728.4763 807.4513 724.3054
## Right_Ent_entorhinal_area Left_Ent_entorhinal_area Right_FO_frontal_operculum
## 1 470.8081 512.4500 438.7602
## 2 371.5236 456.4357 409.6604
## 3 393.1950 415.7966 350.9608
## 4 229.7719 218.7622 317.6053
## 5 504.1409 629.5387 452.3321
## 6 268.3865 362.3349 390.5451
## Left_FO_frontal_operculum Right_FRP_frontal_pole Left_FRP_frontal_pole
## 1 399.2624 420.4775 389.3749
## 2 378.8756 515.9630 394.0116
## 3 355.8374 495.0902 422.7534
## 4 247.9598 389.7580 267.4202
## 5 401.1970 509.3405 411.5239
## 6 415.1297 496.6748 406.0127
## Right_FuG_fusiform_gyrus Left_FuG_fusiform_gyrus Right_GRe_gyrus_rectus
## 1 2190.041 2196.570 385.9814
## 2 1965.592 1918.292 320.9145
## 3 2045.268 1811.830 404.0503
## 4 1571.594 1795.276 324.9839
## 5 2248.107 2305.690 467.7999
## 6 1437.789 1621.597 380.6918
## Left_GRe_gyrus_rectus Right_IOG_inferior_occipital_gyrus
## 1 462.4741 1284.2315
## 2 398.4675 1228.1011
## 3 448.3582 1299.8986
## 4 337.7166 1072.2421
## 5 526.6272 1370.1291
## 6 406.5202 630.6016
## Left_IOG_inferior_occipital_gyrus Right_ITG_inferior_temporal_gyrus
## 1 1460.3330 3045.044
## 2 1291.8200 2611.533
## 3 1406.7272 2566.334
## 4 1310.3683 2421.563
## 5 1282.0306 2942.766
## 6 863.5296 1846.332
## Left_ITG_inferior_temporal_gyrus Right_LiG_lingual_gyrus
## 1 2934.649 2052.104
## 2 2571.158 1845.108
## 3 2416.326 2023.537
## 4 2119.568 1478.046
## 5 3149.072 2052.361
## 6 1962.119 1477.693
## Left_LiG_lingual_gyrus Right_LOrG_lateral_orbital_gyrus
## 1 1890.704 390.2485
## 2 1661.376 413.4683
## 3 1958.926 452.7257
## 4 1511.658 320.5236
## 5 1976.760 417.0670
## 6 1517.256 404.9983
## Left_LOrG_lateral_orbital_gyrus Right_MCgG_middle_cingulate_gyrus
## 1 335.3368 1090.5962
## 2 317.9800 993.7555
## 3 402.5816 1362.7835
## 4 309.6721 925.8665
## 5 354.8905 995.1886
## 6 380.1913 1162.9486
## Left_MCgG_middle_cingulate_gyrus Right_MFC_medial_frontal_cortex
## 1 1075.660 442.9942
## 2 1017.042 364.7236
## 3 1348.294 460.3636
## 4 1005.792 308.9991
## 5 1089.013 432.2927
## 6 1266.525 386.9875
## Left_MFC_medial_frontal_cortex Right_MFG_middle_frontal_gyrus
## 1 436.7819 4308.928
## 2 356.3494 4123.326
## 3 419.6581 4710.869
## 4 312.2126 2811.137
## 5 382.9107 3815.229
## 6 366.8386 3706.915
## Left_MFG_middle_frontal_gyrus Right_MOG_middle_occipital_gyrus
## 1 3850.901 1061.5525
## 2 3557.332 1184.9642
## 3 4553.527 1346.0106
## 4 3585.992 971.3984
## 5 3575.704 1166.0505
## 6 3559.475 543.6746
## Left_MOG_middle_occipital_gyrus Right_MOrG_medial_orbital_gyrus
## 1 1144.9400 783.5301
## 2 1354.1675 745.6793
## 3 1767.1266 796.2907
## 4 1331.0365 659.3365
## 5 1265.8821 924.7532
## 6 749.2799 775.4500
## Left_MOrG_medial_orbital_gyrus Right_MPoG_postcentral_gyrus_medial_segment
## 1 818.9714 119.2732
## 2 748.7998 118.9847
## 3 739.9167 202.8169
## 4 586.3209 104.3297
## 5 815.1683 120.5339
## 6 698.9626 106.4237
## Left_MPoG_postcentral_gyrus_medial_segment
## 1 106.2494
## 2 127.7356
## 3 222.2967
## 4 184.2743
## 5 135.7243
## 6 141.6829
## Right_MPrG_precentral_gyrus_medial_segment
## 1 428.1524
## 2 497.7414
## 3 777.1414
## 4 430.6259
## 5 543.9846
## 6 477.5870
## Left_MPrG_precentral_gyrus_medial_segment
## 1 448.2152
## 2 505.4001
## 3 766.3531
## 4 563.4391
## 5 563.1077
## 6 515.3850
## Right_MSFG_superior_frontal_gyrus_medial_segment
## 1 1614.799
## 2 1355.547
## 3 1674.341
## 4 1297.392
## 5 1695.683
## 6 1540.470
## Left_MSFG_superior_frontal_gyrus_medial_segment
## 1 1507.976
## 2 1290.434
## 3 1446.438
## 4 1088.450
## 5 1429.245
## 6 1328.178
## Right_MTG_middle_temporal_gyrus Left_MTG_middle_temporal_gyrus
## 1 3776.873 3412.776
## 2 3237.662 3030.586
## 3 3531.994 2826.949
## 4 3292.460 2727.113
## 5 3378.029 3466.174
## 6 1723.473 2157.270
## Right_OCP_occipital_pole Left_OCP_occipital_pole
## 1 411.0311 382.8387
## 2 491.2327 509.0663
## 3 551.8073 603.8612
## 4 479.8948 546.2559
## 5 457.1943 409.3231
## 6 262.4983 360.4347
## Right_OFuG_occipital_fusiform_gyrus Left_OFuG_occipital_fusiform_gyrus
## 1 869.9975 932.2508
## 2 771.0413 792.2971
## 3 926.9804 921.5601
## 4 824.5367 966.4971
## 5 864.7238 853.2304
## 6 587.9675 671.1180
## Right_OpIFG_opercular_part_of_the_inferior_frontal_gyrus
## 1 795.9392
## 2 678.8026
## 3 632.0436
## 4 462.3425
## 5 707.1959
## 6 635.3059
## Left_OpIFG_opercular_part_of_the_inferior_frontal_gyrus
## 1 643.7469
## 2 580.7676
## 3 604.5363
## 4 515.1460
## 5 600.3105
## 6 586.4046
## Right_OrIFG_orbital_part_of_the_inferior_frontal_gyrus
## 1 321.5030
## 2 298.0486
## 3 303.7324
## 4 217.6333
## 5 307.3435
## 6 273.9788
## Left_OrIFG_orbital_part_of_the_inferior_frontal_gyrus
## 1 230.4304
## 2 235.6734
## 3 249.3638
## 4 158.6697
## 5 226.0364
## 6 222.7689
## Right_PCgG_posterior_cingulate_gyrus Left_PCgG_posterior_cingulate_gyrus
## 1 868.3144 930.4568
## 2 886.3022 983.2542
## 3 928.6333 998.1955
## 4 740.7224 912.8874
## 5 945.8035 1058.3773
## 6 719.5483 883.9857
## Right_PCu_precuneus Left_PCu_precuneus Right_PHG_parahippocampal_gyrus
## 1 2740.016 2498.985 556.8363
## 2 2699.754 2596.997 526.1274
## 3 3560.742 3537.677 545.0616
## 4 2157.984 2483.972 359.2331
## 5 2530.221 2530.690 663.0359
## 6 2246.718 2296.198 451.4971
## Left_PHG_parahippocampal_gyrus Right_PIns_posterior_insula
## 1 696.6645 528.7890
## 2 703.7195 482.5632
## 3 672.7597 487.9964
## 4 504.5110 396.0665
## 5 912.2598 533.1570
## 6 563.6182 475.0022
## Left_PIns_posterior_insula Right_PO_parietal_operculum
## 1 500.4888 536.7314
## 2 493.0755 414.2915
## 3 448.1238 493.8383
## 4 347.5335 305.6484
## 5 522.9435 401.3681
## 6 516.6780 368.8541
## Left_PO_parietal_operculum Right_PoG_postcentral_gyrus
## 1 594.1453 1518.092
## 2 464.4155 1487.181
## 3 545.3961 2146.631
## 4 408.8729 1215.698
## 5 403.9972 1629.955
## 6 406.5298 1485.266
## Left_PoG_postcentral_gyrus Right_POrG_posterior_orbital_gyrus
## 1 1826.470 555.5699
## 2 1788.792 589.3325
## 3 2888.647 595.2095
## 4 1735.881 444.6932
## 5 1801.243 649.7661
## 6 1604.268 496.0492
## Left_POrG_posterior_orbital_gyrus Right_PP_planum_polare
## 1 560.6758 350.5620
## 2 480.2874 305.6652
## 3 573.7225 258.9913
## 4 406.9689 221.4874
## 5 584.8397 355.7933
## 6 512.7895 256.0361
## Left_PP_planum_polare Right_PrG_precentral_gyrus Left_PrG_precentral_gyrus
## 1 400.4917 2114.178 1980.061
## 2 399.3511 2053.886 2253.542
## 3 321.9628 2811.211 2832.379
## 4 274.1484 1837.671 2124.787
## 5 432.0939 2531.505 2468.314
## 6 342.9697 1794.439 1916.948
## Right_PT_planum_temporale Left_PT_planum_temporale Right_SCA_subcallosal_area
## 1 424.0398 537.5284 216.8584
## 2 322.1405 374.5510 208.1503
## 3 334.0758 423.5789 238.4470
## 4 218.0493 303.2508 197.8472
## 5 363.0033 362.6841 245.5656
## 6 222.7591 309.8191 180.2060
## Left_SCA_subcallosal_area Right_SFG_superior_frontal_gyrus
## 1 272.2707 2559.076
## 2 243.0930 2337.149
## 3 239.7276 3654.610
## 4 197.0983 2363.672
## 5 268.9394 3087.998
## 6 178.3619 2205.228
## Left_SFG_superior_frontal_gyrus Right_SMC_supplementary_motor_cortex
## 1 2598.888 1099.5152
## 2 2264.801 931.7497
## 3 3774.769 1499.5520
## 4 2504.778 880.2152
## 5 3093.241 1284.5298
## 6 2294.438 989.5975
## Left_SMC_supplementary_motor_cortex Right_SMG_supramarginal_gyrus
## 1 1192.735 1937.321
## 2 1023.665 1467.564
## 3 1520.260 2194.323
## 4 1074.839 1304.714
## 5 1255.400 1628.800
## 6 1112.322 1369.122
## Left_SMG_supramarginal_gyrus Right_SOG_superior_occipital_gyrus
## 1 2061.067 592.3486
## 2 1495.502 648.7897
## 3 2499.219 891.3454
## 4 1513.370 711.4169
## 5 1730.106 728.4830
## 6 1536.322 480.5952
## Left_SOG_superior_occipital_gyrus Right_SPL_superior_parietal_lobule
## 1 462.5548 2159.874
## 2 518.9065 2031.026
## 3 707.6336 3381.726
## 4 657.7616 1594.188
## 5 605.1031 2371.317
## 6 386.7248 1851.989
## Left_SPL_superior_parietal_lobule Right_STG_superior_temporal_gyrus
## 1 2228.288 1553.7788
## 2 2004.108 1392.9486
## 3 3185.380 1381.0026
## 4 1803.571 1264.0813
## 5 2503.216 1539.4887
## 6 1895.775 847.4437
## Left_STG_superior_temporal_gyrus Right_TMP_temporal_pole
## 1 1451.0920 1886.398
## 2 1281.2290 1724.003
## 3 1224.0699 1535.828
## 4 1166.7331 1289.619
## 5 1439.4310 2327.944
## 6 980.4075 1214.448
## Left_TMP_temporal_pole
## 1 1750.815
## 2 1711.409
## 3 1424.926
## 4 988.262
## 5 2222.024
## 6 1564.378
## Right_TrIFG_triangular_part_of_the_inferior_frontal_gyrus
## 1 761.2683
## 2 655.4947
## 3 562.2262
## 4 448.5125
## 5 609.3857
## 6 557.9730
## Left_TrIFG_triangular_part_of_the_inferior_frontal_gyrus
## 1 539.7367
## 2 588.2388
## 3 581.0683
## 4 415.7482
## 5 547.6799
## 6 541.2018
## Right_TTG_transverse_temporal_gyrus Left_TTG_transverse_temporal_gyrus
## 1 320.1260 415.5840
## 2 256.5999 328.4706
## 3 240.7111 286.9600
## 4 168.2859 242.2923
## 5 256.3081 252.1471
## 6 198.1721 271.1793- Load the
glmnetandcaretpackages.
library(glmnet)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 3.0-2library(caret)## Loading required package: lattice## Loading required package: ggplot2- (1 point) Set the random seed to 4 and then split the data into 2 sets (400 train, and 260 test)
set.seed(4)
groups <- c(rep("train",400),rep("test",260))
groups <- sample(groups,length(groups))
datTrain <- dementia_dat[groups=="train", ]
datTest <- dementia_dat[groups=="test", ]- (4 points) Perform cross-validated lasso in the training data to select the optimal penalty parameter lambda. Use 5 folds and search over the range \(\lambda = 10^{3}\) to \(\lambda = 10^{-3}\). Set
Dementiaas outcome with all other variables as predictors exceptMacCohort_kr. Use thecaretpackage to do CV.
train_control <- trainControl(method="cv", number=5)
grid <- 10^seq(3,-3,length=100)
caret_grid <- data.frame("lambda" = grid, "alpha"= 1)
cv_model <- train(Dementia ~ . - MacCohort_kr, data=datTrain, trControl=train_control, method="glmnet", tuneGrid=caret_grid)- (1 point) What is the optimal value of lambda?
cv_model$bestTune$lambda## [1] 0.01072267- (1 point) Generate the confusion matrix for this final model.
lassoPred <- predict(cv_model, newdata = datTest)
table(true_disease=datTest$Dementia, predictions=lassoPred)## predictions
## true_disease No Yes
## No 111 15
## Yes 11 123- (1 point) How many non-zero coefficients are in the final model?
lasso_coefs <- coefficients(cv_model$finalModel, s = cv_model$bestTune$lambda)
sum(abs(lasso_coefs) > 0)## [1] 40