1
$\begingroup$

I am using linear mixed-effects (LME) models to investigate the longitudinal effects of maternal factors on infant adiposity indices. Infant adiposity was measured at 3-time points (birth, 3 months and 6 months) using different adiposity indices (e.g. fat mass (g), body fat percentage, fat mass index), and these indices are highly correlated (at all time points r> 0.9). But LME models result in different significant predictors for each outcome variable (e.g. maternal prepregnancy BMI is a significant predictor for body fat percentage, but not infant fat mass index, in 0-6 months old infants).

Predictors were added to the model one at a time and compared using ANOVA to decide whether to keep or not in the final model. Case-wise deletion was used to handle missing data, so the sample sizes for all variables are the same.

Is this result possible? If so, Could anyone provide me with an explanation to justify these results? Many thanks in advance!

E.g.

str(bb)
'data.frame':   478 obs. of  30 variables:
 $ infant_id          : Factor w/ 322 levels "P001","P002",..: 1 4 5 6 7 8 9 10 12 13 ...
 $ ethnicity          : Factor w/ 2 levels "Caucasian","Other": 1 1 1 1 1 1 1 1 1 1 ...
 $ smoking_antenatal  : Factor w/ 2 levels "0-3 days","4-7 days": 1 1 1 1 1 1 1 1 1 1 ...
 $ previous_births    : int  2 0 2 0 0 0 1 0 1 1 ...
 $ mode_delivery      : Factor w/ 2 levels "Vaginal","Caesarean": 1 2 1 2 1 1 1 1 1 2 ...
 $ antenatal_multivits: Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 1 2 ...
 $ antenatal_iron     : Factor w/ 2 levels "No","Yes": 2 2 1 2 2 1 1 1 2 1 ...
 $ antenatal_folicacid: Factor w/ 2 levels "No","Yes": 1 2 1 2 1 1 1 1 1 1 ...
 $ gdm_status         : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ prenatal_bmi       : num  22.8 26.2 23.4 24.7 39.9 30 24.8 35.4 27.7 26.2 ...
 $ net_wt_gain        : num  14.6 8.5 14.9 15.8 -2.4 7.3 17.4 -5.7 11.5 12.5 ...
 $ maternal_age       : int  27 30 33 32 27 30 35 34 21 31 ...
 $ gestational_age    : num  40.1 40.4 39 38.6 40.6 39.7 39.1 39.7 39.6 38.1 ...
 $ infant_sex         : Factor w/ 2 levels "Female","Male": 1 2 1 2 1 1 2 2 1 1 ...
 $ time_point         : Factor w/ 3 levels "Birth","3 months",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ weight_pp          : num  3601 3122 3217 3111 3330 ...
 $ pfm                : num  11.6 8.7 13.6 13.4 8.6 11.4 13 9.3 5.6 7.8 ...
 $ pffm               : num  88.4 91.3 86.4 86.6 91.4 88.6 87 90.7 94.4 92.2 ...
 $ fat_mass           : num  419 271 437 416 286 ...
 $ fatfree_mass       : num  3182 2851 2780 2695 3044 ...
 $ length             : num  50.3 49.5 47.5 49 49 ...
 $ infant_age         : int  2 2 1 2 1 1 2 2 1 1 ...
 $ fmi                : num  1.66 1.1 1.93 1.73 1.19 ...
 $ ffmi               : num  12.6 11.6 12.3 11.2 12.7 ...


> cor(bb$fat_mass,bb$pfm)
[1] 0.9574342

Model for fat mass

Models:
mod1: fat_mass ~ 1 + time_point + (1 | infant_id)
mod2: fat_mass ~ 1 + time_point + prenatal_bmi + (1 | infant_id)
     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1  5 6812.1 6833.0 -3401.1   6802.1                         
mod2  6 6812.5 6837.6 -3400.3   6800.5 1.5803      1     0.2087 

Model for percent body fat

Models:
mod1: pfm ~ 1 + time_point + (1 | infant_id)
mod2: pfm ~ 1 + time_point + prenatal_bmi + (1 | infant_id)
     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
mod1  5 2718.0 2738.9 -1354.0   2708.0                           
mod2  6 2715.7 2740.7 -1351.8   2703.7 4.3732      1    0.03651 *`



[![plot of fat mass vs body fat percentatge][1]][1]


  [1]: https://i.sstatic.net/qeJOI.png

#summary outputs for models of fat mass

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: fat_mass ~ 1 + time_point + (1 | infant_id)
   Data: bb

REML criterion at convergence: 6777.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0232 -0.5632  0.0261  0.4737  3.7745 

Random effects:
 Groups    Name        Variance Std.Dev.
 infant_id (Intercept) 23626    153.7   
 Residual              69320    263.3   
Number of obs: 478, groups:  infant_id, 240

Fixed effects:
                   Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)          353.84      19.87  448.13   17.81   <2e-16 ***
time_point3 months  1065.58      28.50  341.12   37.39   <2e-16 ***
time_point6 months  1511.61      33.54  360.24   45.07   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) tm_p3m
tm_pnt3mnth -0.526       
tm_pnt6mnth -0.445  0.353

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: fat_mass ~ 1 + time_point + prenatal_bmi + (1 | infant_id)
   Data: bb

REML criterion at convergence: 6772.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0035 -0.5518  0.0205  0.4594  3.7930 

Random effects:
 Groups    Name        Variance Std.Dev.
 infant_id (Intercept) 23400    153.0   
 Residual              69376    263.4   
Number of obs: 478, groups:  infant_id, 240

Fixed effects:
                   Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)         273.683     66.948  275.646   4.088 5.71e-05 ***
time_point3 months 1065.912     28.506  340.868  37.393  < 2e-16 ***
time_point6 months 1513.064     33.568  359.490  45.074  < 2e-16 ***
prenatal_bmi          2.975      2.373  259.498   1.254    0.211    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) tm_p3m tm_p6m
tm_pnt3mnth -0.165              
tm_pnt6mnth -0.168  0.353       
prenatal_bm -0.955  0.009  0.037

#summary outputs for models of body fat percentage (pfm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: pfm ~ 1 + time_point + (1 | infant_id)
   Data: bb

REML criterion at convergence: 2709.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3233 -0.6333 -0.0290  0.6385  2.6823 

Random effects:
 Groups    Name        Variance Std.Dev.
 infant_id (Intercept)  4.259   2.064   
 Residual              13.391   3.659   
Number of obs: 478, groups:  infant_id, 240

Fixed effects:
                   Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)         10.4433     0.2738 445.5427   38.13   <2e-16 ***
time_point3 months  13.3559     0.3955 322.3390   33.77   <2e-16 ***
time_point6 months  15.0301     0.4652 344.2752   32.31   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) tm_p3m
tm_pnt3mnth -0.531       
tm_pnt6mnth -0.450  0.352

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: pfm ~ 1 + time_point + prenatal_bmi + (1 | infant_id)
   Data: bb

REML criterion at convergence: 2709.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3046 -0.6541 -0.0076  0.6429  2.7259 

Random effects:
 Groups    Name        Variance Std.Dev.
 infant_id (Intercept)  4.076   2.019   
 Residual              13.407   3.662   
Number of obs: 478, groups:  infant_id, 240

Fixed effects:
                    Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)          8.62101    0.91247 251.71169   9.448   <2e-16 ***
time_point3 months  13.36357    0.39537 322.48135  33.800   <2e-16 ***
time_point6 months  15.06227    0.46524 344.41584  32.375   <2e-16 ***
prenatal_bmi         0.06762    0.03231 234.57422   2.093   0.0375 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) tm_p3m tm_p6m
tm_pnt3mnth -0.168              
tm_pnt6mnth -0.171  0.351       
prenatal_bm -0.954  0.009  0.038
Improve this question
$\endgroup$
12
  • 2
    $\begingroup$ Can you provide the p-values and coefficient estimates for each model? Some variation between models is expected with any non-perfect correlation, and if the p-values we're close to the significance threshold, it's likely that some happen to cross it, and some don't. Also, are the sample sizes equal in all models (e.g. maybe you have some missing values for one outcome, resulting in larger SEs)? $\endgroup$
    – juod
    Commented Aug 24, 2020 at 2:19
  • $\begingroup$ Thank you Juod. I am adding each predictor variable one at a time, and compare them using ANOVA to decide whether to keep or not in the model. mod1: fat_mass ~ 1 + time_point + (1 | infant_id) mod2: fat_mass ~ 1 + time_point + prenatal_bmi + (1 | infant_id); Pr(>Chisq)=0.2087 mod3: pfm ~ 1 + time_point + (1 | infant_id) mod4: pfm ~ 1 + time_point + prenatal_bmi + (1 | infant_id); Pr(>Chisq)=0.03651 I am keeping prenatal BMI only in the model for body fat percentage (pfm) $\endgroup$
    – Prabha
    Commented Aug 25, 2020 at 4:48
  • $\begingroup$ @juod, Missing values were handled with case-wise deletion, so the sample sizes are similar for models for all outcome variables. $\endgroup$
    – Prabha
    Commented Aug 25, 2020 at 4:58
  • $\begingroup$ your approach seems reasonable to me. I suggest adding these details into the original post, as well as any lme model outputs that you could share - maybe someone will notice something. My approach would be to try to determine what could have caused such differences: i.e. look at the actual coefficients of each term, compare the random effect variances, plot the distributions of the outcomes, and see if anything stands out. $\endgroup$
    – juod
    Commented Aug 25, 2020 at 21:50
  • 1
    $\begingroup$ You should always try to give the data itself. Maybe just a snippet of the data. To demonstrate the problem. Note that showing code is of use only to those who is familiar with that language/program. $\endgroup$
    – ttnphns
    Commented Aug 27, 2020 at 1:27

1 Answer 1

Reset to default
1
$\begingroup$

Your approach:

Predictors were added to the model one at a time and compared using ANOVA to decide whether to keep or not in the final model

is probably the source of your problem. This seems to be an attempt to use automated model selection, which is generally a bad idea. In particular, you seem to be using a forward stepwise approach, which might be the worst of all automated methods. In addition to the usual problems with placing too much emphasis on p-values, your modeling is not taking into account your selection of predictors based on their associations with outcome, so the p-values you get are highly unreliable.

I suspect that the following is happening. Your antenatal maternal predictors are likely to have some high correlations among themselves. So for any particular measure of infant adiposity, one of a set of correlated maternal predictors will happen to have the strongest relationship just by chance, depending on the vagaries of your particular data sample. Once that predictor has been added to the model, it will then dominate other members of that set of correlated predictors as you proceed, and prevent them from being incorporated later. That's a particular problem if your ANOVA is using Type I sums of squares, sometimes the default, which gives primacy to the first predictor specified in the model.

There are much better ways to proceed with this type of data. See this page among others for references. In particular, Frank Harrell's course notes and book provide much useful insight for this type of biomedical analysis. Chapter 7 of the notes discusses the pros and cons of different ways of analyzing repeated measures like you have; mixed models, although often useful, are not the only way.

With over 400 observations and approximately 30 predictors, you should be able to include all of your predictors in a single model without much risk of overfitting. That is a much more reliable way of proceeding than trying to build up stepwise with individual predictors, as you seem to have been doing thus far. Among other things, that minimizes the risk of omitted-variable bias, which occurs when you leave out from your model a predictor that is associated with outcome. Also, instead of simply deleting cases with missing values you should consider multiple imputation to avoid the bias that such deletion can cause.

In response to later comments:

You ask:

is it correct if I say that highly correlated outcome variables cannot have different significant predictors if we take the correct approach?

The answer is "No." This has much to do with arbitrary cutoffs of "statistical significance" based on p-values. This is discussed on many pages on this site, for example here and here and here. The underlying estimates of relationships are generally continuous, not all-or-none. For example, "statistical significance" is a function of the size of the data sample you have. With too small a sample you just might not be able to prove, based on p < 0.05, that a particular relationship is "statistically significant" even if it is, in practice, very important and would become evident with a larger data sample.

In a case like yours, having both outcomes and predictors correlated with each other, the specific relationships that turn out to be "statistically significant" based on a p < 0.05 cutoff can depend on the characteristics of the particular sample that you have. A different sample from the population might end up with different determinations of "significance"; you can see that by repeating analysis on multiple bootstrapped samples of the same data set.

You can also have a situation in which neither of two important correlated predictors passes the "significance" test, although the model would be much worse if you removed both of them, and if you included only one of them at a time either of them might be found to be "significant".

Focus on getting a model that describes your data well. You shouldn't make a claim of significance for a particular predictor if your p-value criterion isn't met, but don't over-interpret that lack of "significance." It might just say more about your data sample than about the underlying reality.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you @EdM. I am new to statistics and R, and learnt this approach from an online course. Just wondering if you could suggest me the best approach to analyse my data, it would be very helpful. Also, can I please ask another question, for fat mass and `percent body fat', intercept model (mod1 with time-point) is same and then in mod2, I add the maternal predictor. So, there can not be an impact of the correlation between maternal factors, on the outcome, until I add other maternal predictors later? I also have checked multicollinearity between the predictors (VIF factors<2 for all). $\endgroup$
    – Prabha
    Commented Aug 28, 2020 at 2:22
  • $\begingroup$ @Prabha I have added a bit to the end of the answer. You should be able to model all of your predictors at once and avoid the many problems that your one-at-a-time approach has led to. $\endgroup$
    – EdM
    Commented Aug 28, 2020 at 12:31
  • $\begingroup$ Thank you very much @EdM. So, as a concluding answer for my main question, is it correct if I say that highly correlated outcome variables cannot have different significant predictors if we take the correct approach? $\endgroup$
    – Prabha
    Commented Aug 29, 2020 at 2:29
  • $\begingroup$ I tried adding all my predictors in one model, I got this error "boundary (singular) fit: see ?isSingular" which is an indication of overfitting? $\endgroup$
    – Prabha
    Commented Aug 29, 2020 at 2:51
  • $\begingroup$ @Prabha I have addressed your next-to-last comment in the answer. A singular fit can happen if some of your predictors are exactly linearly related to each other. Look carefully at predictors that are expressed as in terms of each other (e.g., body mass=fat mass+fat-free mass, or proportional percentage values) and remove one or more to avoid linear dependence among the predictors. You might also have a singular fit if you don't have enough observations under different combination of predictors to estimate a random effect, so check the linear dependence without the random effect first. $\endgroup$
    – EdM
    Commented Aug 29, 2020 at 15:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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