CODEWITHSUNDEEP
rregression
Miss Ewura
Miss Ewura

Reputation: 23

NAs in linear regression model summary

I am trying to run a linear regression model on the built-in ChickWeight data with predictors chick, time, diet, time diet interaction as fixed effects and weight as the outcome.

I am getting some NAs in my model summary. Am I doing something wrong?

model2<-lm(weight~ Chick + Time + Diet + Time*Diet, data = Data)
summary(model2)

Upvotes: 2

Views: 769

Answers (2)

dipetkov
dipetkov

Reputation: 3690

As @zephryl points out, each chick was fed the same diet throughout the experiment, so if we were to include Chick as a fixed effect, the Chick effect is confounded with the Diet effect. Intuitively, for example, for a chick that didn't gain much weight, we can't say whether it's because their diet is not very good or that chick was unfortunately a frail chick.

So we take into account the design of the experiment: we have multiple measurements from each chick, so each chick is a cluster of measurements.

We get very similar estimates as with @zephryl's mixed effects model so all is good.

fit.geeglm <- geepack::geeglm(
  weight ~ Time * Diet,
  id = ChickWeight$Chick,
  family = "gaussian",
  data = ChickWeight
)
broom::tidy(fit.geeglm)
#> # A tibble: 8 × 5
#>   term        estimate std.error  statistic  p.value
#>   <chr>          <dbl>     <dbl>      <dbl>    <dbl>
#> 1 (Intercept)   30.9       3.06  102.       0       
#> 2 Time           6.84      0.736  86.3      0       
#> 3 Diet2         -2.30      5.22    0.194    0.660   
#> 4 Diet3        -12.7       4.87    6.77     0.00925 
#> 5 Diet4         -0.139     4.85    0.000820 0.977   
#> 6 Time:Diet2     1.77      1.42    1.55     0.213   
#> 7 Time:Diet3     4.58      1.29   12.6      0.000387
#> 8 Time:Diet4     2.87      0.968   8.80     0.00301

Created on 2022-03-18 by the reprex package (v2.0.1)

Upvotes: 2

zephryl
zephryl

Reputation: 17069

If you look at the summary for your model, you’ll see:

Coefficients: (3 not defined because of singularities)

This means that some of your predictors are perfectly correlated — in this case, Chick and Diet — which is why some coefficients are shown as NA.

Even if this weren’t the case, I’m not sure why you would include a fixed effect for Chick. If you want to account for the fact that there are multiple observations per chick, you can instead include random intercepts for Chick. Using lme4:

library(lme4)

Data <- ChickWeight

model2_mixed <- lmer(
  weight ~ Time*Diet + (1 | Chick), 
  data = Data
)
summary(model2_mixed)

Output:

Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Time * Diet + (1 | Chick)
   Data: Data

REML criterion at convergence: 5466.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3158 -0.5900 -0.0693  0.5361  3.6024 

Random effects:
 Groups   Name        Variance Std.Dev.
 Chick    (Intercept) 545.7    23.36   
 Residual             643.3    25.36   
Number of obs: 578, groups:  Chick, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)  31.5143     6.1163   5.152
Time          6.7115     0.2584  25.976
Diet2        -2.8807    10.5479  -0.273
Diet3       -13.2640    10.5479  -1.258
Diet4        -0.4016    10.5565  -0.038
Time:Diet2    1.8977     0.4284   4.430
Time:Diet3    4.7114     0.4284  10.998
Time:Diet4    2.9506     0.4340   6.799

Correlation of Fixed Effects:
           (Intr) Time   Diet2  Diet3  Diet4  Tm:Dt2 Tm:Dt3
Time       -0.426                                          
Diet2      -0.580  0.247                                   
Diet3      -0.580  0.247  0.336                            
Diet4      -0.579  0.247  0.336  0.336                     
Time:Diet2  0.257 -0.603 -0.431 -0.149 -0.149              
Time:Diet3  0.257 -0.603 -0.149 -0.431 -0.149  0.364       
Time:Diet4  0.254 -0.595 -0.147 -0.147 -0.432  0.359  0.359

Upvotes: 2

Related Questions