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?


1 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);

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
  • 2
    $\begingroup$ The problem remains that when you replace the NAs with zeros the model y ~ x1 + x1:x2 gives the exact same output as y ~ x1 + x2. Therefore it is not a nested model. $\endgroup$ Jan 17, 2019 at 13:19
  • $\begingroup$ I think the problem is that they are the same thing but I am expecting something different. Thanks for your response. $\endgroup$ Jan 17, 2019 at 13:35
  • $\begingroup$ @Ben I'm also running into the same issue, and I agree with Adam's assessment that the model y ~ x1 + x1:x2 you are suggesting is equivalent to y ~ x1 + x2, i.e. (and please correct me if I'm wrong) the model treats the 0's actually as 0's, and not missing data. Do you have any suggestions on other approaches? $\endgroup$
    – platypus17
    Aug 31, 2021 at 5:10
  • $\begingroup$ @AdamWaring: The issue here is that the data has been formed in a particular way where x2 = NA if and only if x1 = 0. That means that the explanatory term is.na(x2) is equivalent to x1 in the model matrix and so it is redundant in the model; consequently, it is not included in the above answer. You are right that the model reduces to y ~ x1 + x2 in this case, but that is a peculiarity of the fact that occurrence of nesting is already indicated perfectly by the binary model term x1. $\endgroup$
    – Ben
    Aug 31, 2021 at 5:22
  • 1
    $\begingroup$ @platypus17: You don't need a way around it in this case since the missingness indicator is already x1. In more general cases, you would add the model term (x2 == 0) as the missingness indicator. $\endgroup$
    – Ben
    Aug 31, 2021 at 6:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.