How do you deal with “nested” variables in a regression model? in R

Question

A conceptual solution for this scenario has been posted in:

How do you deal with "nested" variables in a regression model?

Problem is I am having trouble using this solution in R - glm() or lm().

I am using the model:

y ~ x1 + x1:x2 

Unfortunately if I encode the unmeaninful/missing data as NA the default na.action removes the rows with NAs and leaves x1 with only one level - making the model equivalent to just:

y ~ x2

If I use argument to glm:

na.action = na.pass

I get an error:

Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, : NA/NaN/Inf in 'x'

If instead I encode the missing variable as 0, the nested model:

y ~ x1 + x1:x2

Gives the exact same output as the non nested model:

y ~ x1 + x2

Here is my short script exploring this:

y = rnorm(100, 100, 10)
x1 = sample(c(0, 1), 100, rep=T)
x2 = sample(100:200, 100, rep=T)
x2[x1 == 0] = NA # when x1 is 0 x2 is NA
df = data.frame(y, x1, x2)

nest_NA = glm(y ~ x1 + x1:x2, data=df) # NAs removed 
add_NA = glm(y ~ x1 + x2,, data=df)

nest_NA = glm(y ~ x1 + x1:x2, data=df, na.action = na.pass) # NAs allowed = error
add_NA = glm(y ~ x1 + x2,, data=df, na.action = na.pass)

x2[x1 == 0] = 0 # when x1 is 0 x2 is also 0
df = data.frame(y, x1, x2)

nest_zero = glm(y ~ x1 + x1:x2, data=df)
add_zero = glm(y ~ x1 + x2, data=df)

Am I missing something here?

Answer 1

The issue here is that if you multiply NA by 0 you get NA instead of 0. Under the default options in R, when you run your glm it forms a model matrix by multiplying the elements, leading to many NA values through this multiplication, and then it removes all the rows with NA values.

What we need to do is to form the design matrix so that in the x1:x2 term we get a zero term when we have 0:NA (instead of this term being set to NA and then removed). As you correctly point out in your question, this can be done by changing the data-frame df to replace the NA values by zeros. However, you may prefer not to alter the raw data, and instead form the design matrix for the model directly.

#Form design matrix for your GLM
options(na.action = 'na.pass');
MAT    <- model.matrix( ~ 1 + x1 + x1:x2, data = df);
MAT[is.na(MAT[, 3]), 3] <- 0;

#Fit your GLM
MODEL <- glm(y ~ MAT, data = df);
summary(MODEL);

Call:
glm(formula = y ~ MAT, data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-25.497   -7.363    1.093    6.944   22.913  

Coefficients: (1 not defined because of singularities)
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    102.32290    1.42498  71.807   <2e-16 ***
MAT(Intercept)        NA         NA      NA       NA    
MATx1           -0.45977    7.35367  -0.063    0.950    
MATx1:x2        -0.02592    0.04761  -0.544    0.587    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 103.5589)

    Null deviance: 10540  on 99  degrees of freedom
Residual deviance: 10045  on 97  degrees of freedom
AIC: 752.76

Number of Fisher Scoring iterations: 2