# Reproducing output from emmeans R package

I'm trying to better understand estimated marginal means for relatively simple linear models. To do so, I'm using the very nice emmeans as a reference but also trying to reproduce the results from scratch. The data I'm simulating has a binary treatment assignment and a blocking factor and is simulated to have a within-block treatment effect.

library(emmeans)
library(dplyr)
n_samples = 1e3L
error_variance  = 1
set.seed(4355)
# blocking factors
block_names = c("Group A", "Group B", "Group C", "Group D")
block_probabilities = c(0.04, 0.27, 0.16, 0.53)
# treatments
treatment_names = c("Control", "Treatment")
treatment_probabilities = c(0.5, 0.5)
# build data matrix
dm = data.frame(Level = sample(x = treatment_names,
size = n_samples,
replace = T,
prob = treatment_probabilities),
Block = sample(x = block_names,
size = n_samples,
replace = T,
prob = block_probabilities))
# design matrix
Xlmbi = model.matrix( ~ Level * Block, dm)
# invent marginal parameter values and simulate response
error = rnorm(n_samples, 0, error_variance)
betas = data.frame(
value = c(4.23,
0.02,
0.03,
0.46,
-0.23,
0.01,
0.62,
-0.83),
betanames = colnames(X)
)
zebra = data.frame(dm, y = Xlmbi %*% betas$value + error)  For complicated reasons, I'm interested in fitting a comprehensive model (lmbi = lm(y ~ Level * Block, zebra)), but only interested in the contrast between treatment and control. I can replicate, by hand, the reference grid given by summary(ref_grid(lmbi)), but when trying to get marginalize over Block and it's Level interactions, I can get neither the mean, nor the standard errors right. From Ch. 6.7 of the 5th edition of Kutner, Nachtsheim, Neter & Li, I know that the predicted mean response using the matrix of each unique set of indicators $$X_h$$ is: $$E(\hat{Y}) = X_h \hat{\beta}$$ with estimated variance-covariance $$V(\hat{Y}) = mse * X'_h (X' X)^{-1} X_h = X'_h V(\hat{\beta}) X_h$$ which is how I get the reference grid: # least squares expectation beta_mean_lmbi = solve(t(Xlmbi) %*% Xlmbi) %*% t(Xlmbi) %*% zebra$y
# least squares variance-covariance
epsilon = zebra\$y - Xlmbi %*% beta_mean_lmbi
mse = c((t(epsilon) %*% epsilon) / (nrow(Xlmbi) - ncol(Xlmbi)))
beta_varcovar_lmbi = mse * solve(t(Xlmbi) %*% Xlmbi)
# each unique row of the design matrix
Xlmbi_h = as.matrix(distinct(data.frame(Xlmbi)))
# reproducing the reference grid
y_hat_estimate = c(Xlmbi_h %*% beta_mean_lmbi)
y_hat_std_error = sqrt(diag(Xlmbi_h %*% beta_varcovar_lmbi %*% t(Xlmbi_h)))


But then I can't marginalize correctly:

data.frame(y_hat_estimate, y_hat_std_error) %>%
cbind(Xlmbi_h) %>%
mutate(ses = y_hat_std_error ^ 2) %>%
group_by(LevelTreatment) %>%
summarize(
group_means = mean(y_hat_estimate),
standard_error = sqrt(sum(ses)))


I suspect that there's an elegant linear algebra solution, and that I'm just missing a weight matrix somewhere, but I don't know what the weights are or where they go.

This is documented. I suggest reading the vignette on "basics", where EMMs are described: EMMs defined

The reference grid consists of combinations of predictors. The predictions for the reference grid are each linear combinations of the regression coefficients. You can find out what these are by doing something like this:

rg <- ref_grid(model)
rg@linfct


Each row corresponds to one prediction.

EMMs are averages of these predictions, and their estimates are thus obtained by averaging together appropriate rows of the above. You may find out what these are via:

emm <- emmeans(rg, "treatment")
emm@linfct


Then the SEs are obtained using the formulas in the question:

vcov.emm <- emm@linfct %*% vcov(model) %*% t(emm@linfct)
sqrt(diag(vcov.emm))


I'll leave it to you to translate all this into your own context.