1
$\begingroup$

I have a dataset where I collected cortisol samples 3 times a day, for 3 days at 2 timepoints. I am interested in looking at changes in cortisol right after awakening, at 0 min after awakening, 30 min and 60 min. However, I know that participants actually woke up earlier and took the cortisol samples much later than they were supposed to. So now I want to use the data I have to predict what the estimated cortisol values would have been if samples were taken precisely at 0, 30 and 60 min after awakening. Below is an generated example of what my data looks like:

    head(data)
    ID Trimester       Date Day sample         wakeup_time    sample_time       Cortisol time_since_wake
    1 1002        T0 01/01/2022   0      0 2022-01-01 05:40:00 2022-01-01 08:00:00 0.71343686             140
    2 1002        T0 01/01/2022   0     30 2022-01-01 05:40:00 2022-01-01 08:30:00 0.34559956             170
    3 1002        T0 01/01/2022   0     60 2022-01-01 05:40:00 2022-01-01 09:00:00 0.42300318             200
    4 1003        T0 01/01/2022   0      0 2022-01-01 03:00:00 2022-01-01 08:00:00 0.08194017             300
    5 1003        T0 01/01/2022   0     30 2022-01-01 03:00:00 2022-01-01 08:30:00 0.40100879             330
    6 1003        T0 01/01/2022   0     60 2022-01-01 03:00:00 2022-01-01 09:00:00 0.02570335             360

Using the wakeup time and sample time, I calculated the time difference. So the time_since_wake variable actually tells us when the sample was taken in relation to awakening. So in first row you can see that participant 1002 on day 0 took sample 0 at 140 min after awakening. I want to use this data and create a mixed effects prediction model, so that I can predict estimated values of cortisol samples at 0, 30 and 60 min.

For this purpose, I am implementing the below model and using the following code to predict:

model = lmer(AUCg ~ bs(time_since_wake, df = 3) + (1|ID/Trimester) + 
                     (1|ID/Trimester:Day), data = data)


> summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: AUCg ~ bs(time_since_wake, df = 3) + (1 | ID/Trimester) + 
           (1 |      ID/Trimester:Day)
   Data: data

REML criterion at convergence: 1.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9225 -0.5841  0.0369  0.5540  2.3695 

Random effects:
 Groups           Name        Variance Std.Dev.
 Trimester.Day.ID (Intercept) 0.00000  0.0000  
 Trimester.ID     (Intercept) 0.00000  0.0000  
 ID               (Intercept) 0.00000  0.0000  
 ID.1             (Intercept) 0.00000  0.0000  
 Residual                     0.05767  0.2401  
Number of obs: 36, groups:  Trimester:Day:ID, 12; Trimester:ID, 4; ID, 2

Fixed effects:
                             Estimate Std. Error t value
(Intercept)                   0.44780    0.19719   2.271
bs(time_since_wake, df = 3)1  0.19440    0.48305   0.402
bs(time_since_wake, df = 3)2  0.06442    0.32564   0.198
bs(time_since_wake, df = 3)3 -0.34822    0.30261  -1.151

Correlation of Fixed Effects:
            (Intr) b(__,d=3)1 b(__,d=3)2
bs(__,d=3)1 -0.898                      
bs(__,d=3)2 -0.117 -0.264               
bs(__,d=3)3 -0.765  0.809     -0.318    
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

predicted_data = expand.grid(ID = unique(data$ID), 
Trimester =unique(data$Trimester), 
    Day = unique(data$Day), time_since_wake = c(0, 30, 60))

    predictions = predict(model, newdata = predicted_data)

    predicted_data$cortisol_predicted = predictions

    head(predicted_data)
    ID Trimester Day time_since_wake cortisol_predicted
    1 1002        T0   0               0          0.3940981
    2 1003        T0   0               0          0.3940981
    3 1002        T1   0               0          0.3940981
    4 1003        T1   0               0          0.3940981
    5 1002        T0   1               0          0.3940981
    6 1003        T0   1               0          0.3940981

But as you can see in this predicted_data dataframe, the value for cortisol_predicted is repeated. I noticed that it generates unique predictions for estimated cortisol value at 0, 30 and 60 min time_since_wake, but it is the same across participants, trimesters and days. In my understanding of the model, I believed that including the random effects that I have specified in the model should allow to factor the day-level, trimester-level and individual-level variations. But why is this not reflected in the predictions? I expected that I would get a unique predicted cortisol value for each combination of ID, Trimester, Day and time_since_wake. What am I doing incorrectly? I would appreciate any and all help.

$\endgroup$
4
  • $\begingroup$ Can you add your model summary to the question text? Does it have non-zero variances for the random effects? $\endgroup$ Sep 21, 2023 at 13:01
  • $\begingroup$ @GeorgeSavva the question text has been edited to include the model summary. The random effects do not have non-zero variances with this randomly-generated example dataset. $\endgroup$
    – magg
    Sep 21, 2023 at 13:46
  • $\begingroup$ What are the variances in your real data? If the random effect variances are zero then there will not be different predictions at different levels of the predictors. You could also check your ranef()s from the model. $\endgroup$ Sep 21, 2023 at 15:38
  • $\begingroup$ @GeorgeSavva Thanks so much for that hint. Being fairly new to mixed effects modeling to create such prediction models, I missed out on this bit. I tried the same model with a different random dataset with greater variability in response, which yielded non-zero variances for random effects and then the prediction works as expected. Thanks! $\endgroup$
    – magg
    Sep 21, 2023 at 16:38

1 Answer 1

3
$\begingroup$

As alluded to in the comments, with the simulated data, the model converges to a singular fit:

boundary (singular) fit:

This can occur for several reasons but in your case it appears that you didn't simulate any variance in the random effects. Hence the estimate for the random effects would be close to zero (near or on the boundary of the parameter space).

If your real data does have variation in the grouping factors then hopefully the model will converge normally and you can then use the predict function.

$\endgroup$

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.