1
$\begingroup$

In Python, I'm trying to validate the LS means from a mixed model that I ran with R's lme4 after using the lsmeans library. I'm using Python's mixedlm() from statsmodels.

I've successfully obtained the fixed effects means (by "hand") by extracting the parameters and alternatively with .predict() and a pandas dataframe. But I'm hung up on the standard errors for my predicted values. I'm not quite sure how to obtain these, though R's lsmeans library spits them out. In truth, I'm not sure how R is calculating them.

Does anyone know how to calculate the LS mean standard errors?

Improve this question
$\endgroup$
3
  • $\begingroup$ standard errors from lsmeans are based the pooled standard deviation, obtained from a weighted average of the variances: sqrt(sum((n - 1) * s^2) / sum(n - 1)) $\endgroup$ – Robert Long Aug 25 '20 at 7:27
  • 2
    $\begingroup$ Your quesion is likely to be closed because it is asking about how to do something in specific software, which is off topic. However, at it's heard this is a statistical question so it would be a good idea to rephrase it to be a general one of how to calculate standard errors for lsmeans $\endgroup$ – Robert Long Aug 25 '20 at 7:29
  • $\begingroup$ Robert, thanks for the response; I edited the question accordingly. In my mixed model I have fixed effects for a binary variable, a categorical time variable with 5 levels, and two continuous variables. So do I need a weighted average of the variances of all the variables in the model? Or is this done at each factor level? Do I still need weighted variances for the binary and time variables? $\endgroup$ – SometimesZero Aug 25 '20 at 11:32
1
$\begingroup$

Each LS mean is a linear combination of the fixed-effect regression coefficients $b$, so that the LS mean is (in vector notation) $a'b = \sum_i a_ib_i$. Meanwhile, you can obtain the covariance matrix $V$ of $b$ from the object, via vcov(). The estimated SE of $a'b$ is equal to $\sqrt{a'Va}$.

In R, this fits together as follows: We can obtain the LS means, SEs, etc. from a fitted model via

library(emmeans)
emm <- emmeans(model, ...)
emm

And you can view the underlying sausage from what these are made via

  • $b$ is in emm@bhat
  • $V$ is in emm@V
  • Each row of emm@linfct is the row vector $a'$ corresponding to each LS mean

So the LS means are obtained via emm@linfct %*% emm@bhat, and their estimated covariance matrix is C <- emm@linfct %*% emm@V %*% t(emm@linfct), so that the SEs are sqrt(diag(C)

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you, Russ. Very helpful. I've successfully obtained the covariance matrix, but I'm getting tripped up on the notation (and probably some basic statistics!) that's preventing me from deriving the SE. What do a' and a represent in the SE formula? $\endgroup$ – SometimesZero Aug 25 '20 at 18:38
  • $\begingroup$ I'll expand my answer. $\endgroup$ – Russ Lenth Aug 25 '20 at 21:33
  • $\begingroup$ Wow, your expanded answer is perfect. Can't thank you enough for taking the time to walk me through the calculation. Thank you, again. $\endgroup$ – SometimesZero Aug 25 '20 at 23:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.