(This article was originally published at The Data Game - Ilan Man, and syndicated at StatsBlogs.)
The moment you've all been waiting for has arrived. No, not the day Ryan Gosling actually eats his cereal, but rather the results of the highly anticipated study which attempted to predict the likelihood of a college student having smoked marijuana given a set of predictors. To recap, the predictors are:GPA score, number of classes skipped per week, ever owning a fake ID, ever cheating on a significant other, ever having smoked a cigarette and ever cheating on an exam.
Drum roll please...
....
Cigarettes take the day!
Indeed, by a wide margin, it appears that someone who has ever smoked a cigarette before is much more likely to smoke marijuana than someone who has partaken in one of the other vices mentioned above. And this is controlling for interaction among variables, that is, people who might commit more than 1 vice.
Without digging very deep into it, the odds of someone having ever gotten pollinated (yes, people say that), given someone smoked a cigarette is around 2 times more likely than someone who owns a fake ID and 3 times more likely than someone who just cheated on an exam.
Wow. Amazing. Or actually quite expected as evidenced by the overwhelming number of you who voted for smoking cigarettes as the key predictor.
However, let's not stop our analysis there. I also posed the question "which of the predictors is least significant?" This type of question is also important when performing an analysis, because you want to know which variables you can drop, since these variables might skew your more significant ones. Also, learning about what isn't important can make you double-take on any preconceived notions you might have. For example, some people (usually older, more conservative types) might argue that student's school grades are highly related to their having smoked marijuana. Well based on my analysis, these people are wrong. Not just wrong, but incorrect to the nth degree. Where n = 4.
It turns out that the GPA and cheated-on-exam variables are basically unrelated to someone's likelihood of smoking weed. There is no relationship between the two. This isn't a real shocker. I think we all know good and bad students who like to ride the magic dragon every now and again. Just make sure not to start too young, because apparently that can have detrimental affects (note that study was for teens who use it 4 times a week). Of course, that study only showed an association, not necessarily a strong correlation, let alone a causation.
(Side note - given how many people are interested in analyzing information from published studies, I'm going to do a post in the future on how to best interpret studies like the one I linked to above)
Nerd Alert
If you're interested in basic model selection, keep reading below. Otherwise, feel free to skip the next several paragraphs and head straight to the very bottom.
End Nerd Alert
Statistically speaking, this means that we can remove them (GPA and cheated) as predictors, and re-fit our model. Doing so still yields smoking cigarettes as the most likely predictor, which makes sense, since removing GPA and exam cheating should not have impacted our model, given that they were not related to the outcome variable. By the way, one of the reasons for removing them is that, generally speaking when performing statistical analysis people like to adhere to Occams razor and be parsimonious as long as it's not to the detriment of your fit.
Sometimes removing variables that are weakly related to the outcome is not necessarily the best bet. If you aren't sure whether or not to remove some variables, there are several diagnostic measures that you can perform. There is no right or best method, just rules of thumb. For this fit, we're going to look at the p-value of the Wald Statistic for each variable, the odds ratio, the pseudo R2, and the AIC/BIC values for our full model (i.e. the one with all the variables) and our reduced model (i.e. the resulting model when we remove GPA and exam-cheating).
p-value
A variable's p-value is usually the first thing people look at. It doesn't always tell the full story, especially if there are dependencies among variables, but it's generally a good place to start. Among the various interpretations, the p-value is the probability of observing a test statistic (i.e. the Wald in this case) at least as large as you would expect, if your null hypothesis is true. Our null hypothesis (basically our baseline assumption) is that there is NO relationship between the individual variable and the outcome. Our goal is to see if this is true or false. If false, then we assume that the alternative hypothesis, which says there IS a relationship, is true. Therefore, if the probability of our null hypothesis - our p-value - is small (usually we use a threshold of 5%, but this is quite arbitrary), we think there's such a small chance of it being true, that we reject it. We therefore conclude that the alternative hypothesis is true. If this doesn't make sense to you, don't worry. The concept of p-values is confusing to many people, including grad students. For our purposes, just know that p-values greater than 5% mean no relationship, and imply we should remove the variable from the model.
Let's check out our p-values:
Call: glm(formula = SmokedMJ ~ GPA + SkipClass + FakeID + ChtdSO + SmokeCig + ChtdExam, family = "binomial", data = students)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.2549 -0.9049 -0.8021 1.2100 1.6149
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.76141 0.53781 -1.416 0.156847
GPA -0.05746 0.15963 -0.360 0.718867 SkipClass 0.28251 0.08547 3.305 0.000949 *** FakeIDYes 0.68280 0.19300 3.538 0.000403 *** ChtdSOYes 0.58411 0.16701 3.497 0.000470 *** SmokeCigYes 1.36876 0.24946 5.487 4.09e-08 *** ChtdExamYes 0.17041 0.23966 0.711 0.477070
Looking at the right most column, we see that GPA and ChtdExamYes have larger than 5% p-values. You can ignore the first row, the Intercept. Okay, so the p-values tell us to drop these two variables. But let's continue with the diagnostics before doing so.
odds ratio
Briefly, the odds ratio tells you how significant a predictor is in influencing the outcome variable. The further away the odds ratio is from 1.0 (can be bigger or smaller), the more sensitive the outcome variable is to that predictor. The odds ratio is simply the exponentiated value of the beta parameter for that particular variable. In math notation it's:
Let's look at the odds ratios for each variable:
exp(fit.all$coef)
## output
(Intercept) GPA SkipClass FakeIDYes ChtdSOYes SmokeCigYes ChtdExamYes
0.46700 0.94415 1.32645 1.97942 1.79338 3.93046 1.18578
As you can see, smoking cigarettes and having a fake ID are the largest, while the GPA and cheating on the exam are closest to 1. This doesn't mean we can remove them, but it's a sign to consider it, which is consistent with the p-value approach.
pseudo R2
Traditional R2 is used in linear regression to test how much of the variability in the outcome can be explained by the predictors. In other words, how strong is their relationship? In logistic regression, however, we cannot rely on this measure as it's currently defined, and need to tweak it a bit. There are many different pseudo R2's, but the one we'll use is known as Nagelkerke's R2. It looks like this:
This has many of the same properties as the traditional R2, but modified for logistic regression. Our hope is that the reduced model does not have a significantly lower psuedo R2. Let's now test our full and reduced models:
require(rms)
full.model<-lrm(SmokedMJ~SkipClass+FakeID+ChtdSO+SmokeCig+ChtdExam+GPA,data=student)
reduced.model<-lrm(SmokedMJ~SkipClass+FakeID+ChtdSO+SmokeCig, data=student)
print(full.model) ## psuedo R^2 = 0.134
print(reduced.model) ## psuedo R^2 = 0.135
Great. The reduced model has a pseudo R2 very close to the full model (albeit neither is very high). Let's keep rolling!
AIC and BIC
Without going into the nitty gritty of these two measures (perhaps a topic for another post!) the idea is to select the model which has the lowest value for each. Again, let's see how our models did:
InfoCriteria<-as.data.frame(rbind(c(AIC(reduced.model),BIC(reduced.model)),c(AIC(full.model),BIC(full.model))))
row.names(InfoCriteria)=c("Reduced Model","Full Model")
names(InfoCriteria)=c("AIC","BIC")
print(InfoCriteria) ## output
Reduced Model Full Model
AIC 1136.964 1170.58
BIC 1133.618 1157.63
Unsurprisingly by now, we see that the reduced model has a lower BIC and AIC. More likely than not, this is the best model to use.
You should note that there are tons of different diagnostic measures to use when comparing models and you should never use an absolute rule such as "lower AIC" or "higher R2" = "better model". The analyst should always assess a variety of measures (hopefully measures which provide non-overlapping bits of information) and assess the variables in the model critically. Any prudent analysis is supplemented by a content/domain expert. These are the folks with knowledge of the stuff you're analyzing. Remember, these diagnostic measures are just used to help weed out other potential models.
And with that last marijuana reference, we are brought back full circle to the original intent of this post, which was to determine which predictors of those in our data set, are most likely to tell us if a student has smoked marijuana or not.
I hope you guys enjoyed this. Stay tuned for more (lighter) content!
Please comment on the article here: The Data Game - Ilan Man
