Regression models for data sciense, страница 11

In R, if our variable is a factor variable, it will create thedummy variables for us and pick one of the levels to be the reference level. Let’s go through anexample to see.Insect SpraysLet’s consider a model with factors. Consider the InsectSprays dataset in R. The data modelsthe number of dead insects from different pesticides. Since it’s not clear from the documentation,let’s assume (probably accurately) that these were annoying bad insects, like fleas, mosquitoes orcockroaches, and not good ones like butterflies or ladybugs. After getting over that mental hurdle,let’s plot the data.65Multivariable examples and tricksrequire(datasets);data(InsectSprays); require(stats); require(ggplot2)g = ggplot(data = InsectSprays, aes(y = count, x = spray, fill = spray))g = g + geom_violin(colour = "black", size = 2)g = g + xlab("Type of spray") + ylab("Insect count")gHere’s the plot.

There’s probably better ways to model this data, but let’s use a linear model just toillustrate factor variables.Insect spray datasetFirst, let’s set Spray A as the reference (the default, since it has the lowest alphanumeric factor level).> summary(lm(count ~ spray, data = InsectSprays))$coefEstimate Std. Error t value Pr(>|t|)(Intercept) 14.50001.132 12.8074 1.471e-19sprayB0.83331.601 0.5205 6.045e-01sprayC-12.41671.601 -7.7550 7.267e-11sprayD-9.58331.601 -5.9854 9.817e-08sprayE-11.00001.601 -6.8702 2.754e-09sprayF2.16671.601 1.3532 1.806e-0166Multivariable examples and tricksTherefore, 0.8333 is the estimated mean comparing Spray B to Spray A (as B - A), -12.4167 comparesSpray C to Spray A (as C - A) and so on.

The inferencial statistics: standard errors, t value and P-valueall correspond to those comparisons. The intercept, 14.5, is the mean for Spray A. So, its inferentialstatistics are testing whether or not the mean for Spray A is zero. As is often the case, this test isn’tterribly informative and often yields extremely small statistics (since we know the spray kills somebugs). The estimated mean for Spray B is its effect plus the intercept (14.5 + 0.8333); the estimatedmean for Spray C is 14.5 - 12.4167 (its effect plus the intercept) and so on for the rest of the factorlevels.Let’s hard code the factor levels so we can directly see what’s going on. Remember, we simply leaveout the dummy variable for the reference level.> summary(lm(count ~I(1 * (spray == 'B')) + I(1 * (sprayI(1 * (spray == 'D')) + I(1 * (sprayI(1 * (spray == 'F')), data = InsectSprays))$coefEstimate Std.

Error t value(Intercept)14.50001.132 12.8074I(1 * (spray == "B"))0.83331.601 0.5205I(1 * (spray == "C")) -12.41671.601 -7.7550I(1 * (spray == "D")) -9.58331.601 -5.9854I(1 * (spray == "E")) -11.00001.601 -6.8702I(1 * (spray == "F"))2.16671.601 1.3532== 'C')) +== 'E')) +Pr(>|t|)1.471e-196.045e-017.267e-119.817e-082.754e-091.806e-01Of course, it’s identical. You might further ask yourself, what would happen if I included a dummyvariable for Spray A? Would the world implode? No, it just realizes that one of the dummy variablesis redundant and drops it.> summary(lm(countI(1 * (spray ==I(1 * (spray ==I(1 * (spray ==(Intercept)I(1 * (sprayI(1 * (sprayI(1 * (sprayI(1 * (sprayI(1 * (spray==========~'B')) + I(1 * (spray == 'C')) +'D')) + I(1 * (spray == 'E')) +'F')) + I(1 * (spray == 'A')), data = InsectSprays))$coef"B"))"C"))"D"))"E"))"F"))Estimate Std.

Error t value Pr(>|t|)14.50001.132 12.8074 1.471e-190.83331.601 0.5205 6.045e-01-12.41671.601 -7.7550 7.267e-11-9.58331.601 -5.9854 9.817e-08-11.00001.601 -6.8702 2.754e-092.16671.601 1.3532 1.806e-01However, if we drop the intercept, then the Spray A term is no longer redundant. The each coefficientis the mean for that Spray.67Multivariable examples and tricks> summary(lm(count ~ spray - 1, data = InsectSprays))$coefEstimate Std.

Error t value Pr(>|t|)sprayAsprayBsprayCsprayDsprayEsprayF14.50015.3332.0834.9173.50016.6671.1321.1321.1321.1321.1321.13212.80713.5431.8404.3433.09114.7211.471e-191.002e-207.024e-024.953e-052.917e-031.573e-22So, for example, 14.5 is the mean for Spray A (as we already knew), 15.33 is the mean for Spray B (14.5+ 0.8333 from our previous model formulation), 2.083 is the mean for Spray C (14.5 - 12.4167 fromour previous model formluation) and so on. This is a nice trick if you want your model formulatedin the terms of the group means, rather than the group comparisons relative to the reference group.Also, if there’s no other covariates, the estimated coefficients for this mode are exactly the empiricalmeans of the groups.

We can use dplyr to check this really easily and grab the mean for each group.> library(dplyr)> summarise(group_by(InsectSprays, spray), mn = mean(count))Source: local data frame [6 x 2]123456spraymnA 14.500B 15.333C 2.083D 4.917E 3.500F 16.667Often your lowest alphanumeric level isn’t the level that you’re most interested in as a referencegroup.

There’s an easy fix for that with factor variables; use the relevel function. Here we give asimple example. We created a variable spray2 that has Spray C as the reference level.Multivariable examples and tricks68> spray2 <- relevel(InsectSprays$spray, "C")> summary(lm(count ~ spray2, data = InsectSprays))$coefEstimate Std. Error t value Pr(>|t|)2.0831.132 1.8401 7.024e-02(Intercept)spray2A12.4171.601 7.7550 7.267e-11spray2B13.2501.601 8.2755 8.510e-12spray2D2.8331.601 1.7696 8.141e-02spray2E1.4171.601 0.8848 3.795e-01spray2F14.5831.601 9.1083 2.794e-13Now the intercept is the mean for Spray C and all of the coefficients are interpreted with respect toSpray C. So, 12.417 is the comparison between Spray A and Spray C (as A - C) and so on.Summary of dummy variablesIf you haven’t seen this before, it might seem rather strange.

However, it’s essential to understandhow dummy variables are treated, as otherwise huge interpretation errors can be made. Here wegive a brief bullet summary of dummy variables to help solidify this information.• If we treat a variable as a factor, R includes an intercept and omits the alphabetically first levelof the factor.– The intercept is the estimated mean for the reference level.– The intercept t-test tests for whether or not the mean for the reference level is 0.– All other t-tests are for comparisons of the other levels versus the reference level.– Other group means are obtained the intercept plus their coefficient.• If we omit an intercept, then it includes terms for all levels of the factor.– Group means are now the coefficients.– Tests are tests of whether the groups are different than zero.• If we want comparisons between two levels, neither of which is the reference level, we couldrefit the model with one of them as the reference level.Other thoughts on this dataWe don’t suggest that this is in anyway a thorough analysis of this data.

For example, the data arecounts which are bounded from below by 0. This clearly violates the assumption of normality ofthe errors. Also there are counts near zero, so both the actual assumption and the intent of thisassumption are violated. Furthermore, the variance does not appear to be constant (look back atthe violin plots). Perhaps taking logs of the counts would help. But, since there are 0 counts, maybelog(Count + 1). Also, we’ll cover Poisson GLMs for fitting count data.Multivariable examples and tricks69Further analysis of the swiss datasetWatch this video before beginning.⁸³Then watch this video.⁸⁴Let’s create some dummy variables in the swiss dataset to illustrate them in a more multivariablecontext.

Just to remind ourselves of the dataset, here’s the first few rows.> spray2 <- relevel(InsectSprays$spray, "C")> library(datasets); data(swiss)> head(swiss)Fertility Agriculture Examination Education Catholic Infant.MortalityCourtelary80.217.015129.9622.2Delemont83.145.16984.8422.2Franches-Mnt92.539.75593.4020.2Moutier85.836.512733.7720.3Neuveville76.943.517155.1620.6Porrentruy76.135.39790.5726.6Let’s create a binary variable out of the variable Catholic to illustrate dummy variables inmultivariable models.

However, it should be noted that this isn’t patently absurd, since the variableis highly bimodal anyway. Let’s just split at majority Catholic or not:> spray2 <- relevel(InsectSprays$spray, "C")> library(dplyr)> swiss = mutate(swiss, CatholicBin = 1 * (Catholic > 50))Since we’re interested in Agriculture as a variable and Fertility as an outcome, let’s plot those twocolor coded by the binary Catholic variable:g = ggplot(swiss, aes(x = Agriculture, y = Fertility, colour = factor(CatholicBi\n)))g = g + geom_point(size = 6, colour = "black") + geom_point(size = 4)g = g + xlab("% in Agriculture") + ylab("Fertility")g⁸³https://youtu.be/Xjjbv42KCaM?list=PLpl-gQkQivXjqHAJd2t-J_One_fYE55tC⁸⁴https://youtu.be/HB4owlrqvDE?list=PLpl-gQkQivXjqHAJd2t-J_One_fYE55tC70Multivariable examples and tricksPlot of the Swiss dataset color coded by majority catholic.Our model is:Yi = β0 + Xi1 β1 + Xi2 β2 + ϵiwhere Yi is Fertility, Xi1 is ‘Agriculture and Xi2 is CatholicBin.

Let’s first fit the model withXi2 removed.> summary(lm(Fertility ~ Agriculture, data = swiss))$coef(Intercept)AgricultureEstimate Std. Error t value Pr(>|t|)60.30444.25126 14.185 3.216e-180.19420.076712.532 1.492e-02This model just assumes one line through the data (linear regression). Now let’s add our secondvariable. Notice that the model isYi = β0 + Xi1 β1 + ϵiwhen Xi2 = 0 and71Multivariable examples and tricksYi = (β0 + β2 ) + Xi1 β1 + ϵiwhen Xi2 = 1. Thus, the coefficient in front of the binary variable is the change in the interceptbetween non-Catholic and Catholic majority provinces. In other words, this model fits parallel linesfor the two levels of the factor variable. If the factor variable had 4 levels, it would fit 4 parallel lines,where the coefficients for the factors are the change in the intercepts to the reference level.## Parallel linessummary(lm(Fertility ~ Agriculture + factor(CatholicBin), data = swiss))$coefEstimate Std.