Testing contrast of contrast (first/second difference) in outcome

Question

I want to conduct a contrast of a contrast (i.e., testing for an interaction effect through 1st/2nd differences), where the contrast refers to an outcome (predicted probability) **

It involves 3 steps:

(1) estimate predicted prob (several ways of doing this; I extracted from plot)

(2) calc if there is a difference in predict prob (1st difference) using “pairs”

(3) calc if there is a difference in the difference (2nd difference) using pairs.

I failed in steps (2) and (3). See my code in reprex below with fictious data. Is there a better way to do this?

** My recent S/O post showed HERE showed how to do for “differences in marginal effects of regresser on probability”. But this is a parallel question for “differences in probability of an outcome” instead.

suppressPackageStartupMessages({
  library(emmeans)})

# create ex. data set.  1 row per respondent (dataset shows 2 resp). 
cedata.1 <- data.frame( id    =  c(1,1,1,1,1,1,2,2,2,2,2,2),    
                        QES    = c(1,1,2,2,3,3,1,1,2,2,3,3),   # Choice set   
                        Alt    = c(1,2,1,2,1,2,1,2,1,2,1,2),   # Alt 1 or Alt 2 in  choice set 
                        Choice = c(0,1,1,0,1,0,0,1,0,1,0,1),   # Dep variable.  if  Chosen (1) or not (0)
                        LOC    = c(0,0,1,1,0,1,0,1,1,0,0,1),   # Indep variable per Choice set, binary categorical 
                        SIZE   = c(1,1,1,0,0,1,0,0,1,1,0,1),   # Indep variable per Choice set, binary categorical 
                        gender = c(1,1,1,1,1,1,0,0,0,0,0,0)    # Indep variable per indvidual, binary categorical 
)


# estimate model
glm.model <- glm(Choice ~  LOC*SIZE, data=cedata.1, family = binomial(link = "logit"))

# Plot interaction on response scale (i.e., predict prob) 
zzz <- emmip(glm.model, LOC ~ SIZE, type = "response")   

# (1) Estimate predicted prob   (extract values where yvar=predicted prob)
zzz$data
#>   LOC SIZE      yvar        SE  df tvar xvar
#> 1   0    0 0.3333333 0.2721655 Inf    0    0
#> 2   1    0 0.5000000 0.3535534 Inf    1    0
#> 3   0    1 0.6666667 0.2721655 Inf    0    1
#> 4   1    1 0.5000000 0.2500000 Inf    1    1

# (2) calc 1st diff. 
### I tried following, but got an error ->  pairs(zzz$data, simple = "SIZE", by= NULL?

# (3) calc 2nd diff  
### I tried following, but got an error ->  pairs(pairs(zzz$data, simple = "SIZE"), by = NULL)   

Answer 1

emm <- regrid(emmeans(glm.model, ~ SIZE | LOC))
pairs(pairs(emm), by = NULL)

emm contains the predicted probabilities for the four factor combinations. They are not marginal results because we didn't average over anything. pairs(emm) obtains the SIZE differences at each location, then disregard the by variable the next time around.

BTW you could also do

contrast(emm, interaction = "pairwise")