linearcontrasts

Emmeans for linear contrasts

library(lme4)

Loading required package: Matrix

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(emmeans)

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

setwd("D:/Canna_d/EM5_stuff/EM5_data")
study_time <- read.csv("Em5_1_data.csv", header = T)

# Make condsFile an ordered factor
study_time$condsFile <- factor(study_time$condsFile,
                            levels = c("noDelay", "longDelay", "totalResponse"),
                            ordered = TRUE)

# Model for study time
studyTimeGlmm <- glmer(mean_StudyTime ~ condsFile + age.ct + age.ct*condsFile + (1|participant), 
                         data = study_time, 
                         family = inverse.gaussian(link = "identity"), nAGQ = 0)

# Model for MU
MU_LMM <- lmer(mean_MU ~ condsFile + age.ct + age.ct*condsFile + (1|participant), 
                    data = study_time)

# Study time contrasts
st_contrasts <- emmeans(studyTimeGlmm, "condsFile")

NOTE: Results may be misleading due to involvement in interactions

st_result <- contrast(st_contrasts, "poly")
st_result

 contrast  estimate    SE  df z.ratio p.value
 linear         0.5 0.336 Inf   1.491  0.1359
 quadratic     -1.1 0.641 Inf  -1.718  0.0857

# MU contrasts
MU_contrasts <- emmeans(MU_LMM, "condsFile")

NOTE: Results may be misleading due to involvement in interactions

MU_result <- contrast(MU_contrasts, "poly")
MU_result

 contrast  estimate    SE  df t.ratio p.value
 linear       0.740 0.259 159   2.852  0.0049
 quadratic    0.175 0.446 158   0.392  0.6958

Degrees-of-freedom method: kenward-roger 

Just testing things. Condition study time on age, and also try it without the interaction just to see if that is the issue.

# Get quantiles for age
quantile(study_time$age)

  0%  25%  50%  75% 100% 
 4.5  5.1  5.7  6.3  6.9 

quantile(study_time$age.ct)

          0%          25%          50%          75%         100% 
-1.192592593 -0.592592593  0.007407407  0.607407407  1.207407407 

# Use -.6, 0, .6

# do the contrasts, but separately for 
st_contrasts_new <- emmeans(studyTimeGlmm, ~ condsFile | age.ct, 
                        at = list(age.ct = c(-.6, 0, .6)))  

contrast(st_contrasts_new, "poly")

age.ct = -0.6:
 contrast  estimate    SE  df z.ratio p.value
 linear      -0.795 0.397 Inf  -2.002  0.0452
 quadratic   -2.310 0.880 Inf  -2.626  0.0087

age.ct =  0.0:
 contrast  estimate    SE  df z.ratio p.value
 linear       0.500 0.336 Inf   1.491  0.1359
 quadratic   -1.101 0.641 Inf  -1.718  0.0857

age.ct =  0.6:
 contrast  estimate    SE  df z.ratio p.value
 linear       1.796 0.490 Inf   3.662  0.0003
 quadratic    0.108 0.798 Inf   0.135  0.8923

^ That makes sense. If you look at the age plot, the younger kids have an opposite trend, so we have a negative estimate. Kids in the middle age are in the middle, and then older kids have a positive estimate because they are showing the expected linear increase across delay conditions.

# OK no age model
st_noAge <- glmer(mean_StudyTime ~ condsFile + (1|participant), 
                         data = study_time, 
                         family = inverse.gaussian(link = "identity"), nAGQ = 0)

st_contrasts_noAge <- emmeans(st_noAge, "condsFile")
st_result_noAge <- contrast(st_contrasts_noAge, "poly")
st_result_noAge

 contrast  estimate    SE  df z.ratio p.value
 linear       0.383 0.336 Inf   1.139  0.2547
 quadratic   -1.061 0.656 Inf  -1.617  0.1058

^ That is poor model specification anyway (we know age has an effect) I just wanted to see what it looked like.

Possible linear contrasts method (old)

library(lme4)
setwd("D:/Canna_d/EM5_stuff/EM5_data")
study_time <- read.csv("Em5_1_data.csv", header = T)

# Make condsFile an ordered factor
study_time$condsFile <- factor(study_time$condsFile,
                            levels = c("noDelay", "longDelay", "totalResponse"),
                            ordered = TRUE)
# Polynomial contrasts
contrasts(study_time$condsFile) <- contr.poly(3)

# Model for study time
studyTimeGlmm <- glmer(mean_StudyTime ~ condsFile + age.ct + age.ct*condsFile + (1|participant), 
                         data = study_time, 
                         family = inverse.gaussian(link = "identity"), nAGQ = 0)

# Model for MU
MU_LMM <- lmer(mean_MU ~ condsFile + age.ct + age.ct*condsFile + (1|participant), 
                    data = study_time)

# Summaries
summary(studyTimeGlmm)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: inverse.gaussian  ( identity )
Formula: mean_StudyTime ~ condsFile + age.ct + age.ct * condsFile + (1 |  
    participant)
   Data: study_time

      AIC       BIC    logLik -2*log(L)  df.resid 
   1196.4    1224.4    -590.2    1180.4       235 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7997 -0.4809  0.0240  0.7019  2.8309 

Random effects:
 Groups      Name        Variance Std.Dev.
 participant (Intercept) 3.99117  1.9978  
 Residual                0.02886  0.1699  
Number of obs: 243, groups:  participant, 82

Fixed effects:
                   Estimate Std. Error t value Pr(>|z|)    
(Intercept)         6.47942    0.28489  22.743  < 2e-16 ***
condsFile.L         0.35386    0.23727   1.491   0.1359    
condsFile.Q        -0.44938    0.26151  -1.718   0.0857 .  
age.ct              0.06039    0.38651   0.156   0.8758    
condsFile.L:age.ct  1.52669    0.34654   4.406 1.06e-05 ***
condsFile.Q:age.ct  0.82251    0.36966   2.225   0.0261 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) cndF.L cndF.Q age.ct cF.L:.
condsFile.L  0.066                            
condsFile.Q -0.047  0.101                     
age.ct       0.011  0.153  0.092              
cndsFl.L:g.  0.173  0.210  0.238  0.041       
cndsFl.Q:g.  0.102  0.236 -0.098 -0.065  0.119

summary(MU_LMM)

Linear mixed model fit by REML ['lmerMod']
Formula: mean_MU ~ condsFile + age.ct + age.ct * condsFile + (1 | participant)
   Data: study_time

REML criterion at convergence: 972.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8939 -0.4581  0.0386  0.4729  3.3342 

Random effects:
 Groups      Name        Variance Std.Dev.
 participant (Intercept) 0.5515   0.7426  
 Residual                2.6985   1.6427  
Number of obs: 243, groups:  participant, 82

Fixed effects:
                   Estimate Std. Error t value
(Intercept)         1.91184    0.13362  14.308
condsFile.L         0.52311    0.18337   2.853
condsFile.Q         0.07132    0.18206   0.392
age.ct              0.75940    0.18022   4.214
condsFile.L:age.ct  0.79834    0.24749   3.226
condsFile.Q:age.ct  0.23387    0.24554   0.952

Correlation of Fixed Effects:
            (Intr) cndF.L cndF.Q age.ct cF.L:.
condsFile.L -0.014                            
condsFile.Q  0.008 -0.012                     
age.ct       0.000  0.000  0.000              
cndsFl.L:g.  0.000  0.000  0.000 -0.015       
cndsFl.Q:g.  0.000  0.000  0.000  0.009 -0.014

In this output, condsFile.L is the overall linear trend; condsFile.Q is a quadratic trend (not of interest here); age.ct is the fixed effect of age (centered); condsFile.L:age.ct is the interaction between age and the linear trend of condition; condsFile.Q:age.ct is the interaction between age and the quadratic trend of condition (not of interest).