txt> print("R Descriptive Statistics from Day-Level-1 Data")
[1] "R Descriptive Statistics from Day-Level-1 Data"
txt> univSummary(data = Example5a, vars = c("lglucAM", "negmood", 
txt+ "WPgluc", "WPnm", "PMgluc", "PMnm"))

 Variable    N  Mean    SD Variance    Min   Max Median
 lglucAM  4347 4.978 0.337    0.113  3.881 6.066  4.971
 negmood  4347 0.402 0.948    0.899 -2.492 3.956  0.394
 WPgluc   4347 0.000 0.166    0.027 -0.566 0.637 -0.001
 WPnm     4347 0.000 0.726    0.527 -2.770 2.601 -0.008
 PMgluc   4347 4.978 0.293    0.086  4.131 5.696  4.949
 PMnm     4347 0.402 0.610    0.372 -1.093 1.982  0.448

txt> corrSummary(data = Example5a, vars = c("lglucAM", "negmood", 
txt+ "WPgluc", "WPnm", "PMgluc", "PMnm"), p.separate = TRUE)

Correlations:
 Variable lglucAM negmood WPgluc   WPnm PMgluc  PMnm
 lglucAM    1.000   0.143  0.492  0.059  0.871 0.152
 negmood    0.143   1.000  0.092  0.766  0.113 0.643
 WPgluc     0.492   0.092  1.000  0.120  0.000 0.000
 WPnm       0.059   0.766  0.120  1.000 -0.000 0.000
 PMgluc     0.871   0.113  0.000 -0.000  1.000 0.175
 PMnm       0.152   0.643  0.000  0.000  0.175 1.000

P-values:
 Variable lglucAM negmood WPgluc  WPnm PMgluc  PMnm
 lglucAM            <.001  <.001 <.001  <.001 <.001
 negmood    <.001          <.001 <.001  <.001 <.001
 WPgluc     <.001   <.001        <.001  1.000 1.000
 WPnm       <.001   <.001  <.001        1.000 1.000
 PMgluc     <.001   <.001  1.000 1.000        <.001
 PMnm       <.001   <.001  1.000 1.000  <.001      

txt> print("Model 0: Empty Means, Random Intercept for TV Negative Mood Predictor")
[1] "Model 0: Empty Means, Random Intercept for TV Negative Mood Predictor"
txt> NegMood = lme(data = Example5a, method = "REML", fixed = negmood ~ 
txt+ 1, random = ~1 | PersonID)
txt> lmeSummary(NegMood, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting negmood:
      AIC       BIC    logLik -2*logLik 
10323.935 10343.066 -5158.967 10317.935 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 PersonID  Intercept    0.347   0.589
 Residual               0.554   0.744
Number of obs: 4347 , groups: PersonID , 207 

Fixed Effects:
 Term         Est    SE       DF t.value p.value
 Intercept  0.402 0.042 4140.000   9.453   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.386 
  Unadjusted ICC: 0.386 

LRT for random effect(s):
 Random Effect Dropped  Npar    LL  -2LL   AIC   BIC      LRT DF p.value p.mix
 <none>                    3 -5159 10318 10324 10343                          
 Intercept | PersonID      2 -5941 11881 11885 11898 1563.156  1   <.001 <.001

txt> print("Model 1a: Empty Means, Random Intercept for TV AM Glucose Outcome")
[1] "Model 1a: Empty Means, Random Intercept for TV AM Glucose Outcome"
txt> Model1a = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1, random = ~1 | PersonID)
txt> lmeSummary(Model1a, ICC1 = TRUE, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting lglucAM:
      AIC       BIC    logLik -2*logLik 
-2219.239 -2200.108  1112.620 -2225.239 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 PersonID  Intercept    0.085   0.291
 Residual               0.029   0.170
Number of obs: 4347 , groups: PersonID , 207 

Fixed Effects:
 Term         Est    SE       DF t.value p.value
 Intercept  4.978 0.020 4140.000 243.896   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.747 
  Unadjusted ICC: 0.747 

LRT for random effect(s):
 Random Effect Dropped  Npar    LL  -2LL   AIC   BIC      LRT DF p.value p.mix
 <none>                    3  1113 -2225 -2219 -2200                          
 Intercept | PersonID      2 -1438  2875  2879  2892 5100.459  1   <.001 <.001

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                 FixedEffect Variance    SD Lower Upper
 Intercept | PersonID       4.978    0.085 0.291 4.407 5.549

txt> print("Model 1b: Add AR1 Residual Correlation Across Days")
[1] "Model 1b: Add AR1 Residual Correlation Across Days"
txt> Model1b = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1, random = ~1 | PersonID, correlation = corAR1(form = ~StudyDay | 
txt+ PersonID))
txt> lmeSummary(Model1b, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting lglucAM:
      AIC       BIC    logLik -2*logLik 
-3123.530 -3098.021  1565.765 -3131.530 

Model for the Variance:
 Groups    Name      Variance Std.Dev  Corr
 PersonID  Intercept    0.082   0.286      
 Residual               0.032   0.178      
 Residual  AR1 Phi                    0.480
Number of obs: 4347 , groups: PersonID , 207 

Fixed Effects:
 Term         Est    SE       DF t.value p.value
 Intercept  4.978 0.020 4140.000 244.224   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.722 
  Unadjusted ICC: 0.722 

LRT for each random effect and residual correlation, tested separately (full model vs. that piece dropped):
 Random Effect Dropped       Npar   LL  -2LL   AIC   BIC     LRT DF p.value p.mix
 <none>                         4 1566 -3132 -3124 -3098                         
 Intercept | PersonID           3 1290 -2579 -2573 -2554 552.427  1   <.001 <.001
 Residual correlation (AR1)     3 1113 -2225 -2219 -2200 906.290  1   <.001      

txt> V = getVarCov(Model1b, type = "marginal", individuals = 1)
txt> round(V[[1]][1:7, 1:7], 3)
      1     2     3     4     5     6     7
1 0.114 0.097 0.089 0.086 0.084 0.083 0.082
2 0.097 0.114 0.097 0.089 0.086 0.084 0.083
3 0.089 0.097 0.114 0.097 0.089 0.086 0.084
4 0.086 0.089 0.097 0.114 0.097 0.089 0.086
5 0.084 0.086 0.089 0.097 0.114 0.097 0.089
6 0.083 0.084 0.086 0.089 0.097 0.114 0.097
7 0.082 0.083 0.084 0.086 0.089 0.097 0.114
txt> round(cov2cor(V[[1]])[1:7, 1:7], 3)
      1     2     3     4     5     6     7
1 1.000 0.856 0.786 0.753 0.737 0.729 0.726
2 0.856 1.000 0.856 0.786 0.753 0.737 0.729
3 0.786 0.856 1.000 0.856 0.786 0.753 0.737
4 0.753 0.786 0.856 1.000 0.856 0.786 0.753
5 0.737 0.753 0.786 0.856 1.000 0.856 0.786
6 0.729 0.737 0.753 0.786 0.856 1.000 0.856
7 0.726 0.729 0.737 0.753 0.786 0.856 1.000
txt> MeansPlot = data.frame(StudyDay = 1:21)
txt> MeansPlot$Model1b = predict(Model1b, newdata = MeansPlot, level = 0)
txt> ggplot(MeansPlot, aes(x = StudyDay, y = Model1b)) + geom_line() + 
txt+ geom_point() + scale_x_continuous(breaks = 1:21) + ylim(4.9, 
txt+ 5.1) + labs(x = "Day in Study", y = "Predicted Log Glucose", 
txt+ title = "Model 1b: Same Mean for Every Day") + theme_classic(base_size = 14)
txt> print("Model 1c: Test Saturated Means for Day in Study as a Factor")
[1] "Model 1c: Test Saturated Means for Day in Study as a Factor"
txt> Model1c = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1 + StudyDayF, random = ~1 | PersonID, correlation = corAR1(form = ~StudyDay | 
txt+ PersonID))
txt> lmeSummary(Model1c, predMeans = TRUE)

Model Fit for Predicting lglucAM:
      AIC       BIC    logLik -2*logLik 
-2959.940 -2807.003  1503.970 -3007.940 

Model for the Variance:
 Groups    Name      Variance Std.Dev  Corr
 PersonID  Intercept    0.082   0.286      
 Residual               0.032   0.178      
 Residual  AR1 Phi                    0.481
Number of obs: 4347 , groups: PersonID , 207 

Overall F-test for All Fixed Slopes:
 Model   NumDF    DenDF F.value p.value
 Model1c    20 2135.160   0.930   0.548

Fixed Effects:
 Term                    Est    SE       DF t.value p.value
 Intercept             4.984 0.023  356.988 212.725   <.001
 StudyDayF (ref = 1)                                       
 StudyDayF2           -0.007 0.013 3553.743  -0.584   0.559
 StudyDayF3           -0.003 0.015 4093.263  -0.199   0.842
 StudyDayF4           -0.020 0.016 3512.913  -1.235   0.217
 StudyDayF5           -0.005 0.017 2931.161  -0.284   0.776
 StudyDayF6            0.004 0.017 2571.112   0.222   0.824
 StudyDayF7            0.014 0.017 2367.311   0.830   0.406
 StudyDayF8           -0.020 0.017 2254.172  -1.155   0.248
 StudyDayF9           -0.021 0.017 2191.863  -1.197   0.232
 StudyDayF10          -0.008 0.017 2157.811  -0.452   0.651
 StudyDayF11          -0.004 0.017 2139.365  -0.232   0.817
 StudyDayF12          -0.010 0.017 2129.461  -0.567   0.571
 StudyDayF13          -0.002 0.017 2124.190  -0.114   0.909
 StudyDayF14           0.010 0.017 2121.407   0.569   0.569
 StudyDayF15          -0.004 0.017 2119.949  -0.242   0.809
 StudyDayF16          -0.019 0.017 2119.189  -1.098   0.272
 StudyDayF17          -0.010 0.017 2118.795  -0.553   0.580
 StudyDayF18          -0.002 0.017 2118.593  -0.135   0.892
 StudyDayF19          -0.001 0.017 2118.489  -0.057   0.955
 StudyDayF20          -0.001 0.017 2118.435  -0.040   0.968
 StudyDayF21          -0.020 0.017 2118.408  -1.124   0.261

Model-Predicted Means for StudyDayF:
 Level  Mean    SE      DF
 1     4.984 0.023 356.988
 2     4.977 0.023 356.988
 3     4.981 0.023 356.988
 4     4.964 0.023 356.988
 5     4.979 0.023 356.988
 6     4.988 0.023 356.988
 7     4.998 0.023 356.988
 8     4.964 0.023 356.988
 9     4.963 0.023 356.988
 10    4.976 0.023 356.988
 11    4.980 0.023 356.988
 12    4.974 0.023 356.988
 13    4.982 0.023 356.988
 14    4.994 0.023 356.988
 15    4.980 0.023 356.988
 16    4.965 0.023 356.988
 17    4.974 0.023 356.988
 18    4.982 0.023 356.988
 19    4.983 0.023 356.988
 20    4.983 0.023 356.988
 21    4.964 0.023 356.988

txt> MeansPlot = data.frame(StudyDay = 1:21, StudyDayF = factor(1:21))
txt> MeansPlot$Model1c = predict(Model1c, newdata = MeansPlot, level = 0)
txt> ggplot(MeansPlot, aes(x = StudyDay, y = Model1c)) + geom_line() + 
txt+ geom_point() + scale_x_continuous(breaks = 1:21) + ylim(4, 
txt+ 6) + labs(x = "Day in Study", y = "Predicted Log Glucose", 
txt+ title = "Model 1c: Separate Mean for Each Day") + theme_classic(base_size = 14)
txt> print("Model 1d: Test Saturated Means for Day of Week as a Factor")
[1] "Model 1d: Test Saturated Means for Day of Week as a Factor"
txt> Model1d = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1 + DayName, random = ~1 | PersonID, correlation = corAR1(form = ~StudyDay | 
txt+ PersonID))
txt> lmeSummary(Model1d)

Model Fit for Predicting lglucAM:
      AIC       BIC    logLik -2*logLik 
-3069.165 -3005.409  1544.582 -3089.165 

Model for the Variance:
 Groups    Name      Variance Std.Dev  Corr
 PersonID  Intercept    0.082   0.286      
 Residual               0.032   0.178      
 Residual  AR1 Phi                    0.481
Number of obs: 4347 , groups: PersonID , 207 

Overall F-test for All Fixed Slopes:
 Model   NumDF    DenDF F.value p.value
 Model1d     6 3116.585   1.126   0.344

Fixed Effects:
 Term                         Est    SE       DF t.value p.value
 Intercept                  4.972 0.021  235.944 235.733   <.001
 DayName (ref = Saturday)                                       
 DayNameSunday             -0.006 0.007 3413.692  -0.796   0.426
 DayNameMonday              0.004 0.009 4118.959   0.482   0.630
 DayNameTuesday             0.003 0.009 4034.304   0.292   0.770
 DayNameWednesday           0.007 0.009 3965.812   0.731   0.465
 DayNameThursday            0.013 0.009 4101.565   1.450   0.147
 DayNameFriday              0.015 0.008 4006.466   1.930   0.054

txt> emmSat(Model1d, pairwise ~ DayName)
emmeans:
 DayName   emmean    SE      df lower.CL upper.CL
 Saturday   4.972 0.021 235.944    4.931    5.014
 Sunday     4.966 0.021 238.148    4.925    5.008
 Monday     4.976 0.021 238.652    4.935    5.018
 Tuesday    4.975 0.021 238.744    4.933    5.016
 Wednesday  4.979 0.021 238.652    4.937    5.021
 Thursday   4.985 0.021 238.148    4.943    5.027
 Friday     4.987 0.021 235.944    4.946    5.029
Degrees-of-freedom method: appx-satterthwaite
Confidence level used: 0.95

contrasts:
 contrast             estimate    SE       df t.ratio p.value
 Saturday - Sunday       0.006 0.007 3413.692   0.796   0.985
 Saturday - Monday      -0.004 0.009 4118.959  -0.482   0.999
 Saturday - Tuesday     -0.003 0.009 4034.304  -0.292   1.000
 Saturday - Wednesday   -0.007 0.009 3965.812  -0.731   0.991
 Saturday - Thursday    -0.013 0.009 4101.565  -1.450   0.775
 Saturday - Friday      -0.015 0.008 4006.466  -1.930   0.460
 Sunday - Monday        -0.010 0.007 3438.167  -1.370   0.818
 Sunday - Tuesday       -0.008 0.009 4127.963  -0.962   0.962
 Sunday - Wednesday     -0.013 0.009 3995.922  -1.346   0.830
 Sunday - Thursday      -0.019 0.009 3927.488  -1.996   0.417
 Sunday - Friday        -0.021 0.009 4101.565  -2.337   0.227
 Monday - Tuesday        0.001 0.007 3459.349   0.206   1.000
 Monday - Wednesday     -0.003 0.009 4130.797  -0.299   1.000
 Monday - Thursday      -0.009 0.009 3995.922  -0.947   0.965
 Monday - Friday        -0.011 0.009 3965.812  -1.187   0.899
 Tuesday - Wednesday    -0.004 0.007 3459.349  -0.568   0.998
 Tuesday - Thursday     -0.010 0.009 4127.963  -1.178   0.903
 Tuesday - Friday       -0.013 0.009 4034.304  -1.359   0.823
 Wednesday - Thursday   -0.006 0.007 3438.167  -0.857   0.979
 Wednesday - Friday     -0.008 0.009 4118.959  -0.971   0.960
 Thursday - Friday      -0.002 0.007 3413.692  -0.310   1.000
Degrees-of-freedom method: appx-satterthwaite
P value adjustment: tukey method for comparing a family of 7 estimates

txt> MeansPlot = data.frame(DayName = factor(levels(Example5a$DayName), 
txt+ levels = levels(Example5a$DayName)))
txt> MeansPlot$Model1d = predict(Model1d, newdata = MeansPlot, level = 0)
txt> ggplot(MeansPlot, aes(x = DayName, y = Model1d, group = 1)) + 
txt+ geom_line() + geom_point() + ylim(4, 6) + labs(x = "Day of Week", 
txt+ y = "Predicted Log Glucose", title = "Model 1d: Separate Mean for Each Day of Week") + 
txt+ theme_classic(base_size = 14)
txt> print("Model 2a: Add WPnm and PMnm0 Fixed Effects of Negative Mood")
[1] "Model 2a: Add WPnm and PMnm0 Fixed Effects of Negative Mood"
txt> Model2a = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1 + WPnm + PMnm0, random = ~1 | PersonID, correlation = corAR1(form = ~StudyDay | 
txt+ PersonID))
txt> lmeSummary(Model2a)

Model Fit for Predicting lglucAM:
      AIC       BIC    logLik -2*logLik 
-3196.769 -3158.509  1604.384 -3208.769 

Model for the Variance:
 Groups    Name      Variance Std.Dev  Corr
 PersonID  Intercept    0.080   0.282      
 Residual               0.031   0.177      
 Residual  AR1 Phi                    0.486
Number of obs: 4347 , groups: PersonID , 207 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model2a     2 386.129  46.581   <.001

Fixed Effects:
 Term         Est    SE       DF t.value p.value
 Intercept  4.944 0.024  205.267 205.319   <.001
 WPnm       0.027 0.003 3487.102   9.313   <.001
 PMnm0      0.084 0.033  205.267   2.538   0.012

txt> PseudoR2(baseModel = Model1b, altModel1 = Model2a)

Pseudo-R2
 Term       Model   Variance Pseudo.R2
 Intercept  Model1b    0.082          
 Intercept  Model2a    0.080     0.028
 Residual   Model1b    0.032          
 Residual   Model2a    0.031     0.012

txt> print("Model 2b: Add Random Within-Person (WPnm) Effect of Negative Mood")
[1] "Model 2b: Add Random Within-Person (WPnm) Effect of Negative Mood"
txt> Model2b = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1 + WPnm + PMnm0, random = ~1 + WPnm | PersonID, correlation = corAR1(form = ~StudyDay | 
txt+ PersonID))
txt> lmeSummary(Model2b, LRT = TRUE, digits = 4)

Model Fit for Predicting lglucAM:
       AIC        BIC     logLik  -2*logLik 
-3196.1189 -3145.1065  1606.0594 -3212.1189 

Model for the Variance:
 Groups    Name      Variance Std.Dev    Corr
 PersonID  Intercept   0.0798  0.2824        
           WPnm        0.0003  0.0164 -0.1873
 Residual              0.0310  0.1761        
 Residual  AR1 Phi                     0.4882
Number of obs: 4347 , groups: PersonID , 207 

Overall F-test for All Fixed Slopes:
 Model   NumDF    DenDF F.value p.value
 Model2b     2 193.4749 40.8822  <.0001

Fixed Effects:
 Term          Est     SE       DF  t.value p.value
 Intercept  4.9438 0.0241 204.5845 205.4476  <.0001
 WPnm       0.0273 0.0031 185.3197   8.6674  <.0001
 PMnm0      0.0842 0.0329 202.3973   2.5597  0.0112

LRT for each random effect and residual correlation, tested separately (full model vs. that piece dropped):
 Random Effect Dropped       Npar   LL  -2LL   AIC   BIC      LRT DF p.value  p.mix
 <none>                         8 1606 -3212 -3196 -3145                           
 WPnm | PersonID                6 1604 -3209 -3197 -3159   3.3502  2  0.1873 0.1272
 Residual correlation (AR1)     7 1140 -2280 -2266 -2221 932.3811  1  <.0001       

txt> print("Model 3: Add Moderation of Both Mood Effects by Sex (0=Men, 1=Women)")
[1] "Model 3: Add Moderation of Both Mood Effects by Sex (0=Men, 1=Women)"
txt> Model3 = lme(data = Example5a, method = "REML", fixed = lglucAM ~ 
txt+ 1 + WPnm + PMnm0 + sexMW + WPnm:sexMW + PMnm0:sexMW, random = ~1 | 
txt+ PersonID, correlation = corAR1(form = ~StudyDay | PersonID))
txt> lmeSummary(Model3)

Model Fit for Predicting lglucAM:
      AIC       BIC    logLik -2*logLik 
-3203.822 -3146.439  1610.911 -3221.822 

Model for the Variance:
 Groups    Name      Variance Std.Dev  Corr
 PersonID  Intercept    0.079   0.280      
 Residual               0.031   0.176      
 Residual  AR1 Phi                    0.487
Number of obs: 4347 , groups: PersonID , 207 

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model3     5 324.485  24.701   <.001

Fixed Effects:
 Term            Est    SE       DF t.value p.value
 Intercept     4.886 0.038  202.926 128.100   <.001
 WPnm          0.045 0.005 3505.954   9.796   <.001
 PMnm0         0.157 0.050  202.926   3.145   0.002
 sexMW         0.096 0.049  202.926   1.962   0.051
 WPnm:sexMW   -0.030 0.006 3491.648  -4.995   <.001
 PMnm0:sexMW  -0.125 0.066  202.926  -1.890   0.060

txt> con1DTable(Model3, L = list(`Intercept: Women` = c(Intercept = 1, 
txt+ sexMW = 1), `Within-Person Effect: Women` = c(WPnm = 1, `WPnm:sexMW` = 1), 
txt+ `Between-Person Effect: Women` = c(PMnm0 = 1, `PMnm0:sexMW` = 1)))

Requested Fixed Effects:
 Name                           Est    SE       DF t.value p.value
 Intercept: Women             4.982 0.031  202.926 162.207   <.001
 Within-Person Effect: Women  0.015 0.004 3470.448   4.034   <.001
 Between-Person Effect: Women 0.032 0.044  202.926   0.730   0.466

txt> conMDSummary(Model3, terms = c("sexMW", "WPnm:sexMW", "PMnm0:sexMW"))
Joint Test for:
sexMW, WPnm:sexMW, PMnm0:sexMW
 Model  NumDF   DenDF F.value p.value
 Model3     3 294.955   9.906   <.001

txt> PseudoR2(baseModel = Model1b, altModel1 = Model2a, altModel2 = Model3)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model1b    0.082                    
 Intercept  Model2a    0.080     0.028          
 Intercept  Model3     0.079     0.041     0.013
 Residual   Model1b    0.032                    
 Residual   Model2a    0.031     0.012          
 Residual   Model3     0.031     0.017     0.005

txt> TotalR2(altModel1 = Model2a, altModel2 = Model3)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model2a    0.027          
 Model3     0.045     0.019

txt> packageVersion("nlme")
[1] ‘3.1.168’
txt> citation("nlme")
To cite package ‘nlme’ in publications use:

  Pinheiro J, Bates D, R Core Team (2025). _nlme: Linear and Nonlinear Mixed Effects Models_. doi:10.32614/CRAN.package.nlme
  <https://doi.org/10.32614/CRAN.package.nlme>, R package version 3.1-168, <https://CRAN.R-project.org/package=nlme>.

  Pinheiro JC, Bates DM (2000). _Mixed-Effects Models in S and S-PLUS_. Springer, New York. doi:10.1007/b98882
  <https://doi.org/10.1007/b98882>.

To see these entries in BibTeX format, use 'print(<citation>, bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.
txt> packageVersion("emmeans")
[1] ‘2.0.2’
txt> citation("emmeans")
To cite package ‘emmeans’ in publications use:

  Lenth R, Piaskowski J (2026). _emmeans: Estimated Marginal Means, aka Least-Squares Means_.
  doi:10.32614/CRAN.package.emmeans <https://doi.org/10.32614/CRAN.package.emmeans>, R package version 2.0.2,
  <https://CRAN.R-project.org/package=emmeans>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
    author = {Russell V. Lenth and Julia Piaskowski},
    year = {2026},
    note = {R package version 2.0.2},
    url = {https://CRAN.R-project.org/package=emmeans},
    doi = {10.32614/CRAN.package.emmeans},
  }
txt> Persons5a = unique(Example5a[, c("PersonID", "sexMW", "PMnm0")])
txt> WPsd = sd(Example5a$WPnm)
txt> BPsd = sd(Persons5a$PMnm0)
txt> BPmean = mean(Persons5a$PMnm0)
txt> print("Check predicted outcomes for the sex-moderated mood effects")
[1] "Check predicted outcomes for the sex-moderated mood effects"
txt> GridPlot = expand.grid(WPnm = c(-WPsd, 0, WPsd), PMnm0 = BPmean + 
txt+ c(-BPsd, BPsd), sexMW = c(0, 1))
txt> GridPlot$lglucAM = predict(Model3, newdata = GridPlot, level = 0)
txt> print(GridPlot)
         WPnm      PMnm0 sexMW  lglucAM
1  -0.7261696 -0.2097031     0 4.820075
2   0.0000000 -0.2097031     0 4.852674
3   0.7261696 -0.2097031     0 4.885274
4  -0.7261696  1.0131096     0 5.012025
5   0.0000000  1.0131096     0 5.044624
6   0.7261696  1.0131096     0 5.077223
7  -0.7261696 -0.2097031     1 4.963913
8   0.0000000 -0.2097031     1 4.974973
9   0.7261696 -0.2097031     1 4.986033
10 -0.7261696  1.0131096     1 5.002799
11  0.0000000  1.0131096     1 5.013860
12  0.7261696  1.0131096     1 5.024920
txt> GridPlot$Group = factor(paste(ifelse(GridPlot$PMnm0 < BPmean, 
txt+ "Happy", "Grumpy"), ifelse(GridPlot$sexMW == 0, "Man", "Woman")), 
txt+ levels = c("Happy Man", "Grumpy Man", "Happy Woman", "Grumpy Woman"))
txt> legendLabels = c("Happy Man (PM = -1 SD)", "Grumpy Man (PM = +1 SD)", 
txt+ "Happy Woman (PM = -1 SD)", "Grumpy Woman (PM = +1 SD)")
txt> Figure1 = ggplot(GridPlot, aes(x = WPnm, y = lglucAM, color = Group, 
txt+ linetype = Group, shape = Group, fill = Group)) + geom_line(linewidth = 1.2) + 
txt+ geom_point(size = 3.5, stroke = 1) + scale_color_manual(values = c("blue", 
txt+ "blue", "red", "red"), labels = legendLabels) + scale_linetype_manual(values = c("dashed", 
txt+ "solid", "dashed", "solid"), labels = legendLabels) + scale_shape_manual(values = c(22, 
txt+ 22, 24, 24), labels = legendLabels) + scale_fill_manual(values = c("white", 
txt+ "blue", "white", "red"), labels = legendLabels) + scale_x_continuous(breaks = c(-WPsd, 
txt+ 0, WPsd), labels = c("1 SD Lower than Usual", "Usual", "1 SD Higher than Usual"), 
txt+ expand = expansion(mult = 0.12)) + scale_y_continuous(limits = c(4.6, 
txt+ 5.2), breaks = seq(4.6, 5.2, 0.1), labels = function(y) sprintf("%.2f", 
txt+ y)) + labs(x = "Within-Person Daily Negative Mood", y = "Log Glucose", 
txt+ title = expression(bold("Log Glucose by Sex and Person Mean Negative Mood")), 
txt+ subtitle = expression(bolditalic("(showing WP and BP mood effects)"))) + 
txt+ guides(color = guide_legend(nrow = 2, byrow = TRUE)) + theme_classic(base_size = 14) + 
txt+ theme(plot.title = element_text(face = "bold", hjust = 0.5, 
txt+ 
txt+ size = 16), plot.subtitle = element_text(face = "bold.italic", 
txt+ 
txt+ hjust = 0.5, size = 14, margin = margin(b = 2)), axis.title = element_text(face = "bold", 
txt+ 
txt+ size = 12), axis.text = element_text(size = 12, color = "black"), 
txt+ 
txt+ panel.grid.major.y = element_line(color = "black", linewidth = 0.4), 
txt+ 
txt+ panel.border = element_rect(color = "black", fill = NA, 
txt+ 
txt+ 
txt+ linewidth = 0.8), legend.position = "top", legend.title = element_blank(), 
txt+ 
txt+ legend.key.width = unit(1.6, "cm"), legend.text = element_text(size = 11), 
txt+ 
txt+ legend.box.margin = margin(0, 0, 0, 0), legend.margin = margin(0, 
txt+ 
txt+ 
txt+ 0, 0, 0), legend.box.spacing = unit(2, "pt"), legend.key.spacing.y = unit(0, 
txt+ 
txt+ 
txt+ "pt"), plot.margin = margin(5.5, 30, 5.5, 5.5))
txt> ggsave("Example5a_Sex_by_Mood_WP_BP.png", Figure1, width = 8.5, 
txt+ height = 5, dpi = 300)
txt> print(Figure1)
