> Example5 = Example5[complete.cases(Example5[, c("cognition", 
+ "age", "grip", "sexMW", "demgroup")]), ]
> print("R Descriptive Statistics for Quantitative Variables")
[1] "R Descriptive Statistics for Quantitative Variables"
> psych::describe(x = Example5[, c("cognition", "age", "grip", 
+ "sexMW")])
          vars   n  mean    sd median trimmed   mad   min   max range  skew kurtosis   se
cognition    1 550 24.82 10.99  25.00   25.18 11.86  0.00 44.00 44.00 -0.26    -0.62 0.47
age          2 550 84.93  3.43  84.33   84.49  2.88 80.02 96.97 16.95  1.10     0.72 0.15
grip         3 550  9.11  2.98   9.00    9.11  2.97  0.00 19.00 19.00 -0.01    -0.17 0.13
sexMW        4 550  0.59  0.49   1.00    0.61  0.00  0.00  1.00  1.00 -0.35    -1.88 0.02
> Example5$age85 = Example5$age - 85
> Example5$grip9 = Example5$grip - 9
> print("R Descriptive Statistics for Categorical Variables")
[1] "R Descriptive Statistics for Categorical Variables"
> prop.table(table(x = Example5$sexMW, useNA = "ifany"))
x
     0      1 
0.4127 0.5873 
> prop.table(table(x = Example5$demgroup, useNA = "ifany"))
x
      1       2       3 
0.72545 0.19818 0.07636 
> Example5$demNF = NA
> Example5$demNC = NA
> Example5$demNF[which(Example5$demgroup == 1)] = 0
> Example5$demNC[which(Example5$demgroup == 1)] = 0
> Example5$demNF[which(Example5$demgroup == 2)] = 1
> Example5$demNC[which(Example5$demgroup == 2)] = 0
> Example5$demNF[which(Example5$demgroup == 3)] = 0
> Example5$demNC[which(Example5$demgroup == 3)] = 1
> print("Quadratic Age?")
[1] "Quadratic Age?"
> ModelQuadAge = lm(data = Example5, formula = cognition ~ 1 + 
+ age85 + I(age85^2))
> obj = LMsummary(ModelQuadAge)

Sums of Squares Table
             SS  DF      MS     F      p    R2
Model  1961.684   2 980.842 8.340 <0.001 0.030
Error 64334.854 547 117.614                   
Total 66296.538 549 120.759                   

Fixed Effects Table
              Est    SE      t      p    LCI    UCI
Intercept  24.571 0.601 40.912 <0.001 23.392 25.751
age85      -0.609 0.178 -3.433 <0.001 -0.958 -0.261
I(age85^2)  0.018 0.032  0.549  0.583 -0.045  0.080

> print("Quadratic Grip?")
[1] "Quadratic Grip?"
> ModelQuadGrip = lm(data = Example5, formula = cognition ~ 1 + 
+ grip9 + I(grip9^2))
> obj = LMsummary(ModelQuadGrip)

Sums of Squares Table
             SS  DF       MS      F      p    R2
Model  3897.774   2 1948.887 17.084 <0.001 0.059
Error 62398.764 547  114.075                    
Total 66296.538 549  120.759                    

Fixed Effects Table
              Est    SE      t      p    LCI    UCI
Intercept  24.579 0.566 43.420 <0.001 23.467 25.691
grip9       0.888 0.153  5.801 <0.001  0.587  1.188
I(grip9^2)  0.016 0.038  0.424  0.671 -0.058  0.090

> print("Demgroup")
[1] "Demgroup"
> ModelDem = lm(data = Example5, formula = cognition ~ 1 + demNF + 
+ demNC)
> obj = LMsummary(ModelDem)

Sums of Squares Table
             SS  DF       MS      F      p    R2
Model 11685.511   2 5842.755 58.523 <0.001 0.176
Error 54611.027 547   99.837                    
Total 66296.538 549  120.759                    

Fixed Effects Table
              Est    SE       t      p     LCI     UCI
Intercept  27.198 0.500  54.372 <0.001  26.215  28.181
demNF      -5.675 1.080  -5.255 <0.001  -7.796  -3.554
demNC     -16.388 1.621 -10.111 <0.001 -19.572 -13.205

> print("R Ask for missing model-implied group difference as beta3-beta2")
[1] "R Ask for missing model-implied group difference as beta3-beta2"
> glhtDem = multcomp::glht(model = ModelDem, linfct = rbind(`Future vs Current` = c(0, 
+ -1, 1)))
> obj = glhtSummary(glhtObject = glhtDem)

Linear Combinations Table
                      Est    SE      t      p     LCI    UCI
Future vs Current -10.713 1.815 -5.904 <0.001 -14.278 -7.149

> print("Pearson Correlation Matrix with p-values -- Must convert data frame to matrix")
[1] "Pearson Correlation Matrix with p-values -- Must convert data frame to matrix"
> Corrs = Hmisc::rcorr(x = as.matrix(Example5[, c("cognition", 
+ "age", "grip", "sexMW")]), type = "pearson")
> Corrs
          cognition   age  grip sexMW
cognition      1.00 -0.17  0.24 -0.24
age           -0.17  1.00 -0.18  0.05
grip           0.24 -0.18  1.00 -0.40
sexMW         -0.24  0.05 -0.40  1.00

n= 550 


P
          cognition age    grip   sexMW 
cognition           0.0000 0.0000 0.0000
age       0.0000           0.0000 0.2858
grip      0.0000    0.0000        0.0000
sexMW     0.0000    0.2858 0.0000       
> R2 = Corrs[["r"]]^2
> R2
          cognition      age    grip    sexMW
cognition   1.00000 0.029054 0.05848 0.055830
age         0.02905 1.000000 0.03391 0.002079
grip        0.05848 0.033906 1.00000 0.162605
sexMW       0.05583 0.002079 0.16260 1.000000
> print("Model 1: Main Effects Only Predicting Cognition (Equation 2.8)")
[1] "Model 1: Main Effects Only Predicting Cognition (Equation 2.8)"
> Model1 = lm(data = Example5, formula = cognition ~ 1 + age85 + 
+ grip9 + sexMW + demNF + demNC)
> obj = LMsummary(Model1, effectsizes = TRUE)

Sums of Squares Table
             SS  DF       MS      F      p    R2
Model 18385.979   5 3677.196 41.753 <0.001 0.277
Error 47910.559 544   88.071                    
Total 66296.538 549  120.759                    

Fixed Effects Table
              Est    SE       t      p     LCI     UCI
Intercept  29.264 0.699  41.895 <0.001  27.892  30.636
age85      -0.406 0.119  -3.412 <0.001  -0.639  -0.172
grip9       0.604 0.150   4.034 <0.001   0.310   0.898
sexMW      -3.657 0.891  -4.103 <0.001  -5.408  -1.906
demNF      -5.722 1.019  -5.615 <0.001  -7.724  -3.720
demNC     -16.480 1.523 -10.822 <0.001 -19.471 -13.489

Effect Sizes for Fixed Effects Table
          Est    SE      p      d     pr   sR2
age85  -0.406 0.119 <0.001 -0.293 -0.145 0.015
grip9   0.604 0.150 <0.001  0.346  0.170 0.022
sexMW  -3.657 0.891 <0.001 -0.352 -0.173 0.022
demNF  -5.722 1.019 <0.001 -0.481 -0.234 0.042
demNC -16.480 1.523 <0.001 -0.928 -0.421 0.156

> print("R Ask for missing model-implied group difference as beta5-beta4")
[1] "R Ask for missing model-implied group difference as beta5-beta4"
> SlopesModel1 = multcomp::glht(model = Model1, linfct = rbind(`Future vs Current` = c(0, 
+ 0, 0, 0, -1, 1)))
> obj = glhtSummary(glhtObject = SlopesModel1, effectsizes = TRUE)

Linear Combinations Table
                      Est    SE      t      p     LCI    UCI
Future vs Current -10.758 1.708 -6.299 <0.001 -14.113 -7.403

Effect Sizes for Linear Combinations Table
                      Est    SE      p      d     pr   sR2
Future vs Current -10.758 1.708 <0.001 -0.540 -0.261 0.053

> Model1NoDem = lm(data = Example5, formula = cognition ~ 1 + age85 + 
+ grip9 + sexMW)
> obj = R2compare(ReducedModel = Model1NoDem, FullModel = Model1, 
+ PredName = "Demgroup")

F-Test and R2 Change for Demgroup Slopes 
 R2reduced R2full R2diff DFnum DFden      F      p   pR2   sR2
     0.099  0.277  0.178     2   544 67.056 <0.001 0.198 0.178

> print("Model 2: Add 2 Interactions -- Sex by Age, Sex by Grip")
[1] "Model 2: Add 2 Interactions -- Sex by Age, Sex by Grip"
> Model2 = lm(data = Example5, formula = cognition ~ 1 + age85 + 
+ grip9 + sexMW + demNF + demNC + age85:sexMW + grip9:sexMW)
> obj = R2compare(ReducedModel = Model1, FullModel = Model2, PredName = "Sex*Age and Sex*Grip")

F-Test and R2 Change for Sex*Age and Sex*Grip Slopes 
 R2reduced R2full R2diff DFnum DFden     F     p   pR2   sR2
     0.277  0.279  0.002     2   542 0.814 0.444 0.003 0.002

> obj = LMsummary(Model2, effectsizes = TRUE)

Sums of Squares Table
             SS  DF       MS      F      p    R2
Model 18529.397   7 2647.057 30.035 <0.001 0.279
Error 47767.141 542   88.131                    
Total 66296.538 549  120.759                    

Fixed Effects Table
                Est    SE       t      p     LCI     UCI
Intercept    29.172 0.758  38.467 <0.001  27.682  30.662
age85        -0.232 0.189  -1.226  0.221  -0.604   0.140
grip9         0.697 0.239   2.918  0.004   0.228   1.166
sexMW        -3.627 0.911  -3.980 <0.001  -5.416  -1.837
demNF        -5.749 1.020  -5.635 <0.001  -7.753  -3.745
demNC       -16.495 1.523 -10.828 <0.001 -19.487 -13.502
age85:sexMW  -0.295 0.244  -1.208  0.227  -0.774   0.184
grip9:sexMW  -0.177 0.307  -0.576  0.565  -0.779   0.426

Effect Sizes for Fixed Effects Table
                Est    SE      p      d     pr   sR2
age85        -0.232 0.189  0.221 -0.105 -0.053 0.002
grip9         0.697 0.239  0.004  0.251  0.124 0.011
sexMW        -3.627 0.911 <0.001 -0.342 -0.169 0.021
demNF        -5.749 1.020 <0.001 -0.484 -0.235 0.042
demNC       -16.495 1.523 <0.001 -0.930 -0.422 0.156
age85:sexMW  -0.295 0.244  0.227 -0.104 -0.052 0.002
grip9:sexMW  -0.177 0.307  0.565 -0.050 -0.025 0.000

> print("Simple slopes of age by sex, simple slopes of sex by age")
[1] "Simple slopes of age by sex, simple slopes of sex by age"
> SlopesModel2 = multcomp::glht(model = Model2, linfct = rbind(`Age Slope for Men` = c(0, 
+ 1, 0, 0, 0, 0, 0, 0), `Age Slope for Women` = c(0, 1, 0, 
+ 0, 0, 0, 1, 0), `Grip Slope for Men` = c(0, 0, 1, 0, 0, 0, 
+ 0, 0), `Grip Slope for Women` = c(0, 0, 1, 0, 0, 0, 0, 1)))
> obj = glhtSummary(glhtObject = SlopesModel2, effectsizes = TRUE)

Linear Combinations Table
                        Est    SE      t      p    LCI    UCI
Age Slope for Men    -0.232 0.189 -1.226  0.221 -0.604  0.140
Age Slope for Women  -0.527 0.154 -3.428 <0.001 -0.829 -0.225
Grip Slope for Men    0.697 0.239  2.918  0.004  0.228  1.166
Grip Slope for Women  0.520 0.193  2.694  0.007  0.141  0.899

Effect Sizes for Linear Combinations Table
                        Est    SE      p      d     pr   sR2
Age Slope for Men    -0.232 0.189  0.221 -0.105 -0.053 0.002
Age Slope for Women  -0.527 0.154 <0.001 -0.294 -0.146 0.016
Grip Slope for Men    0.697 0.239  0.004  0.251  0.124 0.011
Grip Slope for Women  0.520 0.193  0.007  0.231  0.115 0.010

> print("Pred cognition outcomes holding demNF=none, and demNC=none")
[1] "Pred cognition outcomes holding demNF=none, and demNC=none"
> print("Provides predicted outcomes from min,max,by=increment of predictors")
[1] "Provides predicted outcomes from min,max,by=increment of predictors"
> PredModel2 = prediction::prediction(model = Model2, type = "response", 
+ at = list(demNF = 0, demNC = 0, grip9 = seq(-3, 3, by = 6), 
+ 
+ age85 = seq(-5, 5, by = 10), sexMW = 0:1))
> PlotModel2 = summary(PredModel2)
> print(PlotModel2, digits = 6)
 at(demNF) at(demNC) at(grip9) at(age85) at(sexMW) Prediction       SE       z            p   lower   upper
         0         0        -3        -5         0    28.2427 1.594901 17.7081  3.63113e-70 25.1167 31.3686
         0         0         3        -5         0    32.4240 1.145406 28.3078 2.76774e-176 30.1790 34.6689
         0         0        -3         5         0    25.9202 1.566544 16.5461  1.70845e-61 22.8498 28.9906
         0         0         3         5         0    30.1015 1.276913 23.5736 7.18750e-123 27.5988 32.6042
         0         0        -3        -5         1    26.6202 1.118191 23.8065 2.85875e-125 24.4286 28.8119
         0         0         3        -5         1    29.7412 1.115829 26.6539 1.61295e-156 27.5542 31.9282
         0         0        -3         5         1    21.3499 0.957061 22.3078 3.10557e-110 19.4741 23.2257
         0         0         3         5         1    24.4708 1.337982 18.2894  1.00594e-74 21.8484 27.0932
> png(file = "Sex by Grip=x GLM Plot.png")
> PlotModel2age80 = PlotModel2[which(PlotModel2$`at(age85)` == 
+ -5), ]
> plot(y = PlotModel2age80$Prediction, x = (PlotModel2age80$`at(grip9)` + 
+ 9), type = "n", ylim = c(15, 45), xlim = c(6, 12), xlab = "Grip Strength", 
+ ylab = "Predicted Cognition")
> lines(x = (PlotModel2age80$`at(grip9)` + 9)[1:2], y = PlotModel2age80$Prediction[1:2], 
+ type = "l", col = "blue1")
> lines(x = (PlotModel2age80$`at(grip9)` + 9)[3:4], y = PlotModel2age80$Prediction[3:4], 
+ type = "l", col = "red1")
> legend(x = 6, y = 45, legend = c("Men", "Women"), col = 1:2, 
+ lty = 1)
> dev.off()
null device 
          1 
> png(file = "Sex by Age=x GLM Plot.png")
> PlotModel2grip6 = PlotModel2[which(PlotModel2$`at(grip9)` == 
+ -3), ]
> plot(y = PlotModel2grip6$Prediction, x = (PlotModel2grip6$`at(grip9)` + 
+ 9), type = "n", ylim = c(15, 45), xlim = c(80, 90), xlab = "Years of Age", 
+ ylab = "Predicted Cognition")
> lines(x = (PlotModel2grip6$`at(age85)` + 85)[1:2], y = PlotModel2grip6$Prediction[1:2], 
+ type = "l", col = "blue1")
> lines(x = (PlotModel2grip6$`at(age85)` + 85)[3:4], y = PlotModel2grip6$Prediction[3:4], 
+ type = "l", col = "red1")
> legend(x = 80, y = 45, legend = c("Men", "Women"), col = 1:2, 
+ lty = 1)
> dev.off()
null device 
          1 
> print("Model 3: Remove 2 Sex Interactions, Add Age by Grip Interaction (Equation 2.9)")
[1] "Model 3: Remove 2 Sex Interactions, Add Age by Grip Interaction (Equation 2.9)"
> Model3 = lm(data = Example5, formula = cognition ~ 1 + age85 + 
+ grip9 + sexMW + demNF + demNC + age85:grip9)
> obj = LMsummary(Model3, effectsizes = TRUE)

Sums of Squares Table
             SS  DF       MS      F      p    R2
Model 19185.041   6 3197.507 36.854 <0.001 0.289
Error 47111.497 543   86.762                    
Total 66296.538 549  120.759                    

Fixed Effects Table
                Est    SE       t      p     LCI     UCI
Intercept    29.408 0.695  42.319 <0.001  28.043  30.773
age85        -0.334 0.120  -2.775  0.006  -0.570  -0.098
grip9         0.619 0.149   4.164 <0.001   0.327   0.912
sexMW        -3.456 0.887  -3.895 <0.001  -5.199  -1.713
demNF        -5.923 1.014  -5.843 <0.001  -7.914  -3.931
demNC       -16.300 1.513 -10.777 <0.001 -19.272 -13.329
age85:grip9   0.123 0.041   3.035  0.003   0.043   0.203

Effect Sizes for Fixed Effects Table
                Est    SE      p      d     pr   sR2
age85        -0.334 0.120  0.006 -0.238 -0.118 0.010
grip9         0.619 0.149 <0.001  0.357  0.176 0.023
sexMW        -3.456 0.887 <0.001 -0.334 -0.165 0.020
demNF        -5.923 1.014 <0.001 -0.501 -0.243 0.045
demNC       -16.300 1.513 <0.001 -0.925 -0.420 0.152
age85:grip9   0.123 0.041  0.003  0.260  0.129 0.012

> Model3NoDem = lm(data = Example5, formula = cognition ~ 1 + age85 + 
+ grip9 + sexMW + age85:grip9)
> obj = R2compare(ReducedModel = Model3NoDem, FullModel = Model3, 
+ PredName = "Demgroup")

F-Test and R2 Change for Demgroup Slopes 
 R2reduced R2full R2diff DFnum DFden      F      p   pR2   sR2
     0.112  0.289  0.177     2   543 67.701 <0.001 0.200 0.177

> print("R Ask for missing model-implied group difference as beta5-beta4")
[1] "R Ask for missing model-implied group difference as beta5-beta4"
> SlopesModel3 = multcomp::glht(model = Model3, linfct = rbind(`Future vs Current` = c(0, 
+ 0, 0, 0, -1, 1, 0)))
> obj = glhtSummary(glhtObject = SlopesModel3, effectsizes = TRUE)

Linear Combinations Table
                      Est    SE      t      p     LCI    UCI
Future vs Current -10.378 1.700 -6.105 <0.001 -13.717 -7.039

Effect Sizes for Linear Combinations Table
                      Est    SE      p      d     pr   sR2
Future vs Current -10.378 1.700 <0.001 -0.524 -0.253 0.049

> SlopesModel3 = multcomp::glht(model = Model3, linfct = rbind(`Age Slope at Grip = 6` = c(0, 
+ 1, 0, 0, 0, 0, -3), `Age Slope at Grip = 9` = c(0, 1, 0, 
+ 0, 0, 0, 0), `Age Slope at Grip = 12` = c(0, 1, 0, 0, 0, 
+ 0, 3), `Grip Slope at Age = 80` = c(0, 0, 1, 0, 0, 0, -5), 
+ `Grip Slope at Age = 85` = c(0, 0, 1, 0, 0, 0, 0), `Grip Slope at Age = 90` = c(0, 
+ 
+ 0, 1, 0, 0, 0, 5)))
> obj = glhtSummary(glhtObject = SlopesModel3, effectsizes = TRUE)

Linear Combinations Table
                          Est    SE      t      p    LCI    UCI
Age Slope at Grip = 6  -0.703 0.153 -4.584 <0.001 -1.004 -0.402
Age Slope at Grip = 9  -0.334 0.120 -2.775  0.006 -0.570 -0.098
Age Slope at Grip = 12  0.035 0.187  0.188  0.851 -0.333  0.403
Grip Slope at Age = 80  0.004 0.247  0.017  0.986 -0.482  0.490
Grip Slope at Age = 85  0.619 0.149  4.164 <0.001  0.327  0.912
Grip Slope at Age = 90  1.235 0.255  4.833 <0.001  0.733  1.736

Effect Sizes for Linear Combinations Table
                          Est    SE      p      d     pr   sR2
Age Slope at Grip = 6  -0.703 0.153 <0.001 -0.393 -0.193 0.027
Age Slope at Grip = 9  -0.334 0.120  0.006 -0.238 -0.118 0.010
Age Slope at Grip = 12  0.035 0.187  0.851  0.016  0.008 0.000
Grip Slope at Age = 80  0.004 0.247  0.986  0.002  0.001 0.000
Grip Slope at Age = 85  0.619 0.149 <0.001  0.357  0.176 0.023
Grip Slope at Age = 90  1.235 0.255 <0.001  0.415  0.203 0.031

> print("Simple slopes over range of moderator values using reghelper package instead")
[1] "Simple slopes over range of moderator values using reghelper package instead"
> reghelper::simple_slopes(model = Model3, levels = list(age85 = c(-5, 
+ 0, 5, "sstest"), grip9 = c(-3, 0, 3, "sstest")))
   age85  grip9 Test Estimate Std. Error t value  df   Pr(>|t|)
1 sstest     -3        -0.703      0.153  -4.584 543 0.00000567
2 sstest      0        -0.334      0.120  -2.775 543    0.00571
3 sstest      3         0.035      0.187   0.188 543    0.85132
4     -5 sstest         0.004      0.247   0.017 543    0.98605
5      0 sstest         0.619      0.149   4.164 543 0.00003631
6      5 sstest         1.235      0.255   4.833 543 0.00000175
> print("Pred cognition outcomes holding sexMW=0, demNF=none, and demNC=none")
[1] "Pred cognition outcomes holding sexMW=0, demNF=none, and demNC=none"
> print("Provides predicted outcomes from min,max,by=increment of predictors")
[1] "Provides predicted outcomes from min,max,by=increment of predictors"
> PredModel3 = prediction(model = Model3, type = "response", at = list(sexMW = 0, 
+ demNF = 0, demNC = 0, grip9 = seq(-3, 3, by = 3), age85 = seq(-5, 
+ 
+ 5, by = 5)))
> PlotModel3 = summary(PredModel3)
> print(PlotModel3, digits = 6)
 at(sexMW) at(demNF) at(demNC) at(grip9) at(age85) Prediction       SE       z            p   lower   upper
         0         0         0        -3        -5    31.0646 1.260473 24.6452 4.14143e-134 28.5941 33.5351
         0         0         0         0        -5    31.0776 0.916789 33.8983 7.04923e-252 29.2807 32.8745
         0         0         0         3        -5    31.0906 1.092408 28.4606 3.60205e-178 28.9495 33.2317
         0         0         0        -3         0    27.5495 0.930878 29.5952 1.72045e-192 25.7251 29.3740
         0         0         0         0         0    29.4078 0.694906 42.3191  0.00000e+00 28.0458 30.7698
         0         0         0         3         0    31.2661 0.705332 44.3281  0.00000e+00 29.8836 32.6485
         0         0         0        -3         5    24.0345 1.149080 20.9163  3.80801e-97 21.7823 26.2866
         0         0         0         0         5    27.7380 0.921723 30.0936 5.86735e-199 25.9315 29.5445
         0         0         0         3         5    31.4415 1.246179 25.2304 1.86086e-140 28.9991 33.8840
> png(file = "Age by Grip=x GLM Plot.png")
> PlotModel3 = PlotModel3[order(PlotModel3$`at(age85)`), ]
> plot(y = PlotModel3$Prediction, x = (PlotModel3$`at(grip9)` + 
+ 9), type = "n", ylim = c(15, 45), xlim = c(6, 12), xlab = "Grip Strength", 
+ ylab = "Predicted Cognition")
> lines(x = (PlotModel3$`at(grip9)` + 9)[1:3], y = PlotModel3$Prediction[1:3], 
+ type = "l", col = "blue1")
> lines(x = (PlotModel3$`at(grip9)` + 9)[4:6], y = PlotModel3$Prediction[4:6], 
+ type = "l", col = "red1")
> lines(x = (PlotModel3$`at(grip9)` + 9)[7:9], y = PlotModel3$Prediction[7:9], 
+ type = "l", col = "green1")
> legend(x = 6, y = 45, legend = c("Age 80", "Age 85", "Age 90"), 
+ col = 1:3, lty = 1)
> dev.off()
null device 
          1 
> png(file = "Grip by Age=x GLM Plot.png")
> PlotModel3 = PlotModel3[order(PlotModel3$`at(grip9)`), ]
> plot(y = PlotModel3$Prediction, x = (PlotModel3$`at(age85)` + 
+ 85), type = "n", ylim = c(15, 45), xlim = c(80, 90), xlab = "Years of Age", 
+ ylab = "Predicted Cognition")
> lines(x = (PlotModel3$`at(age85)` + 85)[1:3], y = PlotModel3$Prediction[1:3], 
+ type = "l", col = "blue1")
> lines(x = (PlotModel3$`at(age85)` + 85)[4:6], y = PlotModel3$Prediction[4:6], 
+ type = "l", col = "red1")
> lines(x = (PlotModel3$`at(age85)` + 85)[7:9], y = PlotModel3$Prediction[7:9], 
+ type = "l", col = "green1")
> legend(x = 80, y = 45, legend = c("Grip=6", "Grip=9", "Grip=12"), 
+ col = 1:3, lty = 1)
> dev.off()
null device 
          1 
> print("Values for regions of significance")
[1] "Values for regions of significance"
> vcov(Model3)
            (Intercept)      age85      grip9     sexMW     demNF      demNC age85:grip9
(Intercept)   0.4828946  0.0004536 -0.0307533 -0.450672 -0.182029 -0.2263470   0.0019165
age85         0.0004536  0.0144857  0.0033170  0.005024 -0.004131 -0.0011511   0.0009587
grip9        -0.0307533  0.0033170  0.0221243  0.053743 -0.013388 -0.0003011   0.0002029
sexMW        -0.4506716  0.0050240  0.0537427  0.787257 -0.071016  0.0237081   0.0026947
demNF        -0.1820286 -0.0041309 -0.0133877 -0.071016  1.027449  0.2129117  -0.0026791
demNC        -0.2263470 -0.0011511 -0.0003011  0.023708  0.212912  2.2877993   0.0023964
age85:grip9   0.0019165  0.0009587  0.0002029  0.002695 -0.002679  0.0023964   0.0016432
> interactions::johnson_neyman(model = Model3, pred = "age85", 
+ modx = "grip9", digits = 3, plot = TRUE)
JOHNSON-NEYMAN INTERVAL

When grip9 is OUTSIDE the interval [0.665, 9.521], the slope of age85 is p < .05.

Note: The range of observed values of grip9 is [-9.000, 10.000]

> interactions::johnson_neyman(model = Model3, pred = "grip9", 
+ modx = "age85", digits = 3, plot = TRUE)
JOHNSON-NEYMAN INTERVAL

When age85 is OUTSIDE the interval [-14.873, -2.281], the slope of grip9 is p < .05.

Note: The range of observed values of age85 is [-4.984, 11.967]

> print("Simple slope boundaries for age and grip given by regions of significance")
[1] "Simple slope boundaries for age and grip given by regions of significance"
> print(Model3[["coefficients"]], digits = 8)
 (Intercept)        age85        grip9        sexMW        demNF        demNC  age85:grip9 
 29.40780315  -0.33396058   0.61941863  -3.45563720  -5.92254309 -16.30040485   0.12301848 
> RegionSlopesModel3 = multcomp::glht(model = Model3, linfct = rbind(`Age Slope at Grip =  9.665` = c(0, 
+ 1, 0, 0, 0, 0, 0.665), `Age Slope at Grip = 18.521` = c(0, 
+ 1, 0, 0, 0, 0, 9.521), `Grip Slope at Age = 70.127` = c(0, 
+ 0, 1, 0, 0, 0, -14.873), `Grip Slope at Age = 82.719` = c(0, 
+ 0, 1, 0, 0, 0, -2.281)))
> obj = glhtSummary(glhtObject = RegionSlopesModel3, effectsizes = TRUE)

Linear Combinations Table
                              Est    SE      t     p    LCI    UCI
Age Slope at Grip =  9.665 -0.252 0.128 -1.964 0.050 -0.504  0.000
Age Slope at Grip = 18.521  0.837 0.426  1.964 0.050 -0.000  1.675
Grip Slope at Age = 70.127 -1.210 0.616 -1.964 0.050 -2.420 -0.000
Grip Slope at Age = 82.719  0.339 0.172  1.964 0.050  0.000  0.678

Effect Sizes for Linear Combinations Table
                              Est    SE     p      d     pr   sR2
Age Slope at Grip =  9.665 -0.252 0.128 0.050 -0.169 -0.084 0.005
Age Slope at Grip = 18.521  0.837 0.426 0.050  0.169  0.084 0.005
Grip Slope at Age = 70.127 -1.210 0.616 0.050 -0.169 -0.084 0.005
Grip Slope at Age = 82.719  0.339 0.172 0.050  0.169  0.084 0.005

> print("Model 4: Add Sex by Dementia Group Interaction (Equation 2.13)")
[1] "Model 4: Add Sex by Dementia Group Interaction (Equation 2.13)"
> print("Binary Predictors for Sex (0=Men) and Demgroup (0=None)")
[1] "Binary Predictors for Sex (0=Men) and Demgroup (0=None)"
> Model4 = lm(data = Example5, formula = cognition ~ 1 + age85 + 
+ grip9 + sexMW + demNF + demNC + age85:grip9 + sexMW:demNF + 
+ sexMW:demNC)
> obj = R2compare(ReducedModel = Model3, FullModel = Model4, PredName = "Sex*Demgroup")

F-Test and R2 Change for Sex*Demgroup Slopes 
 R2reduced R2full R2diff DFnum DFden     F     p   pR2   sR2
     0.289  0.298  0.009     2   541 3.492 0.031 0.013 0.009

> Model4NoAgeGrip = lm(data = Example5, formula = cognition ~ 1 + 
+ sexMW + sexMW:demNF + sexMW:demNC)
> obj = R2compare(ReducedModel = Model4NoAgeGrip, FullModel = Model4, 
+ PredName = "Age, Grip, and Age*Grip")

F-Test and R2 Change for Age, Grip, and Age*Grip Slopes 
 R2reduced R2full R2diff DFnum DFden      F      p   pR2   sR2
     0.193  0.298  0.106     5   541 16.329 <0.001 0.131 0.106

> obj = LMsummary(Model4, effectsizes = TRUE)

Sums of Squares Table
             SS  DF       MS      F      p    R2
Model 19785.461   8 2473.183 28.767 <0.001 0.298
Error 46511.077 541   85.972                    
Total 66296.538 549  120.759                    

Fixed Effects Table
                Est    SE      t      p     LCI    UCI
Intercept    29.070 0.748 38.838 <0.001  27.600 30.540
age85        -0.335 0.120 -2.793  0.005  -0.570 -0.099
grip9         0.618 0.148  4.173 <0.001   0.327  0.909
sexMW        -2.876 1.011 -2.844  0.005  -4.862 -0.889
demNF        -6.056 1.635 -3.704 <0.001  -9.268 -2.844
demNC       -11.971 2.245 -5.332 <0.001 -16.381 -7.561
age85:grip9   0.122 0.040  3.027  0.003   0.043  0.201
sexMW:demNF   0.164 2.070  0.079  0.937  -3.903  4.231
sexMW:demNC  -7.875 3.025 -2.604  0.009 -13.816 -1.934

Effect Sizes for Fixed Effects Table
                Est    SE      p      d     pr   sR2
age85        -0.335 0.120  0.005 -0.240 -0.119 0.010
grip9         0.618 0.148 <0.001  0.359  0.177 0.023
sexMW        -2.876 1.011  0.005 -0.245 -0.121 0.010
demNF        -6.056 1.635 <0.001 -0.318 -0.157 0.018
demNC       -11.971 2.245 <0.001 -0.459 -0.223 0.037
age85:grip9   0.122 0.040  0.003  0.260  0.129 0.012
sexMW:demNF   0.164 2.070  0.937  0.007  0.003 0.000
sexMW:demNC  -7.875 3.025  0.009 -0.224 -0.111 0.009

> print("Pred cognition outcomes --adjusted cell means-- holding age=85 and grip=9")
[1] "Pred cognition outcomes --adjusted cell means-- holding age=85 and grip=9"
> print("Will need to ignore impossible combinations of demNF and demNC for min:max")
[1] "Will need to ignore impossible combinations of demNF and demNC for min:max"
> PredModel4 = summary(prediction(model = Model4, type = "response", 
+ at = list(age85 = 0, grip9 = 0, sexMW = 0:1, demNF = 0:1, 
+ 
+ demNC = 0:1)))
> print(PredModel4, digits = 6)
 at(age85) at(grip9) at(sexMW) at(demNF) at(demNC) Prediction       SE         z           p    lower    upper
         0         0         0         0         0  29.070146 0.748499 38.837912 0.00000e+00 27.60311 30.53718
         0         0         1         0         0  26.194552 0.638834 41.003693 0.00000e+00 24.94246 27.44664
         0         0         0         1         0  23.014245 1.492760 15.417242 1.25343e-53 20.08849 25.94000
         0         0         1         1         0  20.302921 1.118633 18.149756 1.28974e-73 18.11044 22.49540
         0         0         0         0         1  17.099416 2.140218  7.989567 1.35414e-15 12.90467 21.29417
         0         0         1         0         1   6.348722 1.947880  3.259298 1.11688e-03  2.53095 10.16650
         0         0         0         1         1  11.043514 2.696429  4.095607 4.21063e-05  5.75861 16.32842
         0         0         1         1         1   0.457091 2.317874  0.197203 8.43669e-01 -4.08586  5.00004
> print("Create data frame for plotting and remove last 2 unneeded rows")
[1] "Create data frame for plotting and remove last 2 unneeded rows"
> PredModel4 = data.frame(PredModel4)
> PredModel4$sum = PredModel4$at.demNF. + PredModel4$at.demNC.
> PredModel4 = subset(x = PredModel4, PredModel4$sum < 2)
> PredModel4$demgroup = NA
> PredModel4$demgroup[which(PredModel4$at.demNF. == 0 & PredModel4$at.demNC. == 
+ 0)] = 1
> PredModel4$demgroup[which(PredModel4$at.demNF. == 1 & PredModel4$at.demNC. == 
+ 0)] = 2
> PredModel4$demgroup[which(PredModel4$at.demNF. == 0 & PredModel4$at.demNC. == 
+ 1)] = 3
> png(file = "Sex by Demgroup=x GLM Plot.png")
> plot(y = PredModel4$Prediction, x = PredModel4$demgroup, type = "n", 
+ ylim = c(0, 35), xlim = c(1, 3), xlab = "Dementia Group (1=None, 2=Future, 3=Current)", 
+ ylab = "Predicted Cognition")
> PredModel4 = PredModel4[order(PredModel4$at.sexMW.), ]
> lines(x = PredModel4$demgroup[1:3], y = PredModel4$Prediction[1:3], 
+ type = "l", col = "blue1")
> lines(x = PredModel4$demgroup[4:6], y = PredModel4$Prediction[4:6], 
+ type = "l", col = "red1")
> legend(x = 1, y = 35, legend = c("Sex=Men", "Sex=Women"), col = 1:2, 
+ lty = 1)
> dev.off()
RStudioGD 
        2 
> png(file = "Demgroup by Sex=x GLM Plot.png")
> plot(y = PredModel4$Prediction, x = PredModel4$at.sexMW., type = "n", 
+ ylim = c(0, 35), xlim = c(0, 1), xlab = "Sex (0=Men, 1=Women)", 
+ ylab = "Predicted Cognition")
> PredModel4 = PredModel4[order(PredModel4$demgroup), ]
> lines(x = PredModel4$at.sexMW.[1:2], y = PredModel4$Prediction[1:2], 
+ type = "l", col = "blue1")
> lines(x = PredModel4$at.sexMW.[3:4], y = PredModel4$Prediction[3:4], 
+ type = "l", col = "red1")
> lines(x = PredModel4$at.sexMW.[5:6], y = PredModel4$Prediction[5:6], 
+ type = "l", col = "green1")
> legend(x = 0, y = 36, legend = c("Demgroup=None", "Demgroup=Future", 
+ "Demgroup=Current"), col = 1:3, lty = 1)
> dev.off()
RStudioGD 
        2 
> print("GLHT pred cognition outcomes --adjusted cell means-- holding age=85 and grip=9")
[1] "GLHT pred cognition outcomes --adjusted cell means-- holding age=85 and grip=9"
> Pred2Model4 = multcomp::glht(model = Model4, linfct = rbind(`Yhat for Men   None` = c(1, 
+ 0, 0, 0, 0, 0, 0, 0, 0), `Yhat for Women None` = c(1, 0, 
+ 0, 1, 0, 0, 0, 0, 0), `Yhat for Men Future` = c(1, 0, 0, 
+ 0, 1, 0, 0, 0, 0), `Yhat for Women Future` = c(1, 0, 0, 1, 
+ 1, 0, 0, 1, 0), `Yhat for Men Current` = c(1, 0, 0, 0, 0, 
+ 1, 0, 0, 0), `Yhat for Women Current` = c(1, 0, 0, 1, 0, 
+ 1, 0, 0, 1)))
> obj = glhtSummary(glhtObject = Pred2Model4)

Linear Combinations Table
                          Est    SE      t      p    LCI    UCI
Yhat for Men   None    29.070 0.748 38.838 <0.001 27.600 30.540
Yhat for Women None    26.195 0.639 41.004 <0.001 24.940 27.449
Yhat for Men Future    23.014 1.493 15.417 <0.001 20.082 25.947
Yhat for Women Future  20.303 1.119 18.150 <0.001 18.106 22.500
Yhat for Men Current   17.099 2.140  7.990 <0.001 12.895 21.304
Yhat for Women Current  6.349 1.948  3.259  0.001  2.522 10.175

> print("DF=1 simple slopes for sex per demgroup, demgroup per sex, and interactions")
[1] "DF=1 simple slopes for sex per demgroup, demgroup per sex, and interactions"
> SlopesModel4 = multcomp::glht(model = Model4, linfct = rbind(`Sex Diff for No Dementia` = c(0, 
+ 0, 0, 1, 0, 0, 0, 0, 0), `Sex Diff for Future Dementia` = c(0, 
+ 0, 0, 1, 0, 0, 0, 1, 0), `Sex Diff for Current Dementia` = c(0, 
+ 0, 0, 1, 0, 0, 0, 0, 1), `None-Future Diff for Men` = c(0, 
+ 0, 0, 0, 1, 0, 0, 0, 0), `None-Future Diff for Women` = c(0, 
+ 0, 0, 0, 1, 0, 0, 1, 0), `None-Current Diff for Men` = c(0, 
+ 0, 0, 0, 0, 1, 0, 0, 0), `None-Current Diff for Women` = c(0, 
+ 0, 0, 0, 0, 1, 0, 0, 1), `Future-Current Diff for Men` = c(0, 
+ 0, 0, 0, -1, 1, 0, 0, 0), `Future-Current Diff for Women` = c(0, 
+ 0, 0, 0, -1, 1, 0, -1, 1), `A: Sex effect differ btw None and Future?` = c(0, 
+ 0, 0, 0, 0, 0, 0, 1, 0), `A: None-Future effect differ btw Men and Women?` = c(0, 
+ 0, 0, 0, 0, 0, 0, 1, 0), `B: Sex effect differ btw None and Current?` = c(0, 
+ 0, 0, 0, 0, 0, 0, 0, 1), `B: None-Current effect differ btw Men and Women?` = c(0, 
+ 0, 0, 0, 0, 0, 0, 0, 1), `C: Sex effect differ btw Future and Current?` = c(0, 
+ 0, 0, 0, 0, 0, 0, -1, 1), `C: Future-Current effect differ btw Men and Women?` = c(0, 
+ 0, 0, 0, 0, 0, 0, -1, 1)))
> obj = glhtSummary(glhtObject = SlopesModel4, effectsizes = TRUE)

Linear Combinations Table
                                                       Est    SE      t      p     LCI     UCI
Sex Diff for No Dementia                            -2.876 1.011 -2.844  0.005  -4.862  -0.889
Sex Diff for Future Dementia                        -2.711 1.874 -1.447  0.149  -6.393   0.970
Sex Diff for Current Dementia                      -10.751 2.899 -3.708 <0.001 -16.446  -5.055
None-Future Diff for Men                            -6.056 1.635 -3.704 <0.001  -9.268  -2.844
None-Future Diff for Women                          -5.892 1.278 -4.611 <0.001  -8.402  -3.382
None-Current Diff for Men                          -11.971 2.245 -5.332 <0.001 -16.381  -7.561
None-Current Diff for Women                        -19.846 2.029 -9.783 <0.001 -23.831 -15.861
Future-Current Diff for Men                         -5.915 2.587 -2.287  0.023 -10.996  -0.834
Future-Current Diff for Women                      -13.954 2.239 -6.233 <0.001 -18.352  -9.556
A: Sex effect differ btw None and Future?            0.164 2.070  0.079  0.937  -3.903   4.231
A: None-Future effect differ btw Men and Women?      0.164 2.070  0.079  0.937  -3.903   4.231
B: Sex effect differ btw None and Current?          -7.875 3.025 -2.604  0.009 -13.816  -1.934
B: None-Current effect differ btw Men and Women?    -7.875 3.025 -2.604  0.009 -13.816  -1.934
C: Sex effect differ btw Future and Current?        -8.039 3.415 -2.354  0.019 -14.748  -1.331
C: Future-Current effect differ btw Men and Women?  -8.039 3.415 -2.354  0.019 -14.748  -1.331

Effect Sizes for Linear Combinations Table
                                                       Est    SE      p      d     pr   sR2
Sex Diff for No Dementia                            -2.876 1.011  0.005 -0.245 -0.121 0.010
Sex Diff for Future Dementia                        -2.711 1.874  0.149 -0.124 -0.062 0.003
Sex Diff for Current Dementia                      -10.751 2.899 <0.001 -0.319 -0.157 0.018
None-Future Diff for Men                            -6.056 1.635 <0.001 -0.318 -0.157 0.018
None-Future Diff for Women                          -5.892 1.278 <0.001 -0.396 -0.194 0.028
None-Current Diff for Men                          -11.971 2.245 <0.001 -0.459 -0.223 0.037
None-Current Diff for Women                        -19.846 2.029 <0.001 -0.841 -0.388 0.124
Future-Current Diff for Men                         -5.915 2.587  0.023 -0.197 -0.098 0.007
Future-Current Diff for Women                      -13.954 2.239 <0.001 -0.536 -0.259 0.050
A: Sex effect differ btw None and Future?            0.164 2.070  0.937  0.007  0.003 0.000
A: None-Future effect differ btw Men and Women?      0.164 2.070  0.937  0.007  0.003 0.000
B: Sex effect differ btw None and Current?          -7.875 3.025  0.009 -0.224 -0.111 0.009
B: None-Current effect differ btw Men and Women?    -7.875 3.025  0.009 -0.224 -0.111 0.009
C: Sex effect differ btw Future and Current?        -8.039 3.415  0.019 -0.202 -0.101 0.007
C: Future-Current effect differ btw Men and Women?  -8.039 3.415  0.019 -0.202 -0.101 0.007

> print("Omnibus DF=2 F-test for Dementia Simple Main Effect for Men")
[1] "Omnibus DF=2 F-test for Dementia Simple Main Effect for Men"
> print("Comma separates each fixed effect to be tested jointly, already in model here")
[1] "Comma separates each fixed effect to be tested jointly, already in model here"
> DemforM = multcomp::glht(model = Model4, linfct = rbind(c(0, 
+ 0, 0, 0, 1, 0, 0, 0, 0), c(0, 0, 0, 0, 0, 1, 0, 0, 0)))
> print(summary(DemforM, test = Ftest()), digits = 4)

	 General Linear Hypotheses

Linear Hypotheses:
       Estimate
1 == 0   -6.056
2 == 0  -11.971

Global Test:
      F DF1 DF2        Pr(>F)
1 18.69   2 541 0.00000001419
> print("Omnibus DF=2 F-test for Dementia Simple Main Effect for Women")
[1] "Omnibus DF=2 F-test for Dementia Simple Main Effect for Women"
> print("Comma separates each fixed effect to be tested jointly, formed as a linear combination")
[1] "Comma separates each fixed effect to be tested jointly, formed as a linear combination"
> DemforW = multcomp::glht(model = Model4, linfct = rbind(c(0, 
+ 0, 0, 0, 1, 0, 0, 1, 0), c(0, 0, 0, 0, 0, 1, 0, 0, 1)))
> print(summary(DemforW, test = Ftest()), digits = 4)

	 General Linear Hypotheses

Linear Hypotheses:
       Estimate
1 == 0   -5.892
2 == 0  -19.846

Global Test:
      F DF1 DF2    Pr(>F)
1 53.16   2 541 8.378e-22
> print("Regions of significance for age and grip strength slopes")
[1] "Regions of significance for age and grip strength slopes"
> print("Values for regions of significance")
[1] "Values for regions of significance"
> print(Model4[["coefficients"]], digits = 8)
 (Intercept)        age85        grip9        sexMW        demNF        demNC  age85:grip9  sexMW:demNF  sexMW:demNC 
 29.07014634  -0.33479877   0.61789286  -2.87559408  -6.05590147 -11.97073055   0.12215159   0.16426999  -7.87509954 
> vcov(Model4)
            (Intercept)      age85      grip9     sexMW     demNF      demNC age85:grip9 sexMW:demNF sexMW:demNC
(Intercept)    0.560251  0.0016040 -0.0309178 -0.587371 -0.502778 -0.5097673  0.00195260  0.52140496   0.5131013
age85          0.001604  0.0143730  0.0032782  0.002989 -0.009677 -0.0029569  0.00095028  0.00899367   0.0031988
grip9         -0.030918  0.0032782  0.0219275  0.054020 -0.010633 -0.0006414  0.00020119 -0.00422662   0.0006726
sexMW         -0.587371  0.0029888  0.0540198  1.022601  0.485319  0.5143296  0.00257807 -0.89866720  -0.8821734
demNF         -0.502778 -0.0096769 -0.0106333  0.485319  2.673637  0.5110573 -0.00269673 -2.66391613  -0.5149121
demNC         -0.509767 -0.0029569 -0.0006414  0.514330  0.511057  5.0398171  0.00182362 -0.51427819  -5.0362444
age85:grip9    0.001953  0.0009503  0.0002012  0.002578 -0.002697  0.0018236  0.00162835  0.00007386   0.0010011
sexMW:demNF    0.521405  0.0089937 -0.0042266 -0.898667 -2.663916 -0.5142782  0.00007386  4.28686770   0.8856691
sexMW:demNC    0.513101  0.0031988  0.0006726 -0.882173 -0.514912 -5.0362444  0.00100108  0.88566905   9.1478209
> interactions::johnson_neyman(model = Model4, pred = "age85", 
+ modx = "grip9", digits = 3, plot = FALSE)
JOHNSON-NEYMAN INTERVAL

When grip9 is OUTSIDE the interval [0.680, 9.638], the slope of age85 is p < .05.

Note: The range of observed values of grip9 is [-9.000, 10.000]

> interactions::johnson_neyman(model = Model4, pred = "grip9", 
+ modx = "age85", digits = 3, plot = FALSE)
JOHNSON-NEYMAN INTERVAL

When age85 is OUTSIDE the interval [-15.003, -2.293], the slope of grip9 is p < .05.

Note: The range of observed values of age85 is [-4.984, 11.967]

