How to do a marginal model analysis in R?

Question

Problem. I'd like to do a marginal model analysis in R--I think it is sometimes called population averaged model, marginal multilevel model or marginal linear regression model. However, I can't find any information on stackoverflow, Google, or Youtube about how to do this specifically in R.

Background. I'm talking about marginal models as described by The Analysis Factor here and here and as described on these PowerPoint slides. There's one person on CrossValidated who mentions this analysis in SPSS and R but he doesn't show his actual code, and his question hasn't been answered. Not sure if should be done in nlme package or not.

SPSS Code. I've described the nature of this data elsewhere on CrossValidated, but basically, we're interested in predicting participants emotion (measured twice in 2 different conditions) via personality (measured once). Here's the code I used in SPSS.

MIXED emotion BY condition WITH centeredPersonality
    /FIXED=condition centeredPersonality condition*centeredPersonality
    /METHOD = REML
    /REPEATED= condition | SUBJECT (ID) COVTYPE(UN)
    /PRINT=SOLUTION.

Question. How to do this in R?

Answer 1

I think geeglm for the geepack package can do that. My understanding is that generalized estimating equations are the same thing as marginal models. geeglm has syntax similar to glm, and if you use a gaussian family, you'll get a result similar to a standard marginal model. I'm sure there are other ways, but this should work.

edit: Here is an example you might use, regressing emotion on two variables, condition and personality, and their interaction. Condition is being treated as a factor, and errors are clustered by id. The default family for geeglm is gaussian/Normal, so we don't need to specify that.

> library(geepack)
> dat <- data.frame(id = c(1, 1, 2, 2, 3, 3, 4, 4), 
+                   condition = factor(c(1, 2, 1, 2, 1, 2, 1, 2)), 
+                   personality = c(2.5, 2.5, 4.0, 4.0, 3.3, 3.3, 4.2, 4.2),
+                   emotion = c(5.0, 4.9, 2.6, 2.3, 4.3, 2.9, 1.0, 1.0))
>   
> my_mod <- geeglm(emotion ~ condition*personality, data = dat, id = id)
> summary(my_mod)

Call:
geeglm(formula = emotion ~ condition * personality, data = dat, 
    id = id)

 Coefficients:
                       Estimate Std.err  Wald Pr(>|W|)    
(Intercept)              10.815   1.296 69.68  < 2e-16 ***
condition2               -0.902   1.284  0.49     0.48    
personality              -2.169   0.385 31.77  1.7e-08 ***
condition2:personality    0.129   0.322  0.16     0.69    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Estimated Scale Parameters:
            Estimate Std.err
(Intercept)    0.223  0.0427

Correlation: Structure = independenceNumber of clusters:   4   Maximum cluster size: 2