R: Calculate the expected rank of random coefficients that...
expectedRank
R Documentation
Calculate the expected rank of random coefficients that account for
uncertainty.
Description
expectedRank calculates the expected rank and the percentile expected
rank of any random term in a merMod object. A simple ranking of the estimated
random effects (as produced by ranef) is not satisfactory
because it ignores any amount of uncertainty.
Usage
expectedRank(merMod, groupFctr = NULL, term = NULL)
Arguments
merMod
An object of class merMod
groupFctr
The name of the grouping factor over which the random
coefficient of interest varies. This is the variable to the right of the
pipe, |, in the [g]lmer formula. This parameter is optional if only
a single grouping factor is included in the model, but required if there
are two or more.
term
The name of the random coefficient of interest. This is the
variable to the left of the pipe, |, in the [g]lmer formula. Partial
matching is attempted on the intercept term so the following character
strings will all return rankings based on the intercept (provided that
they do not match the name of another random coefficient for that factor):
c("(Intercept)", "Int", "intercep", ...).
Details
Inspired by Lingsma et al. (2010, see also Laird and Louis 1989),
expectedRank sums the probability that each level of the grouping factor is
greater than every other level of the grouping factor, similar to a
two-sample t-test.
where φ is the standard normal distribution function, θ
is the estimated random effect and var(θ) is the posterior
variance of the estimated random effect. We add one to the sum so that the
minimum rank is one instead of zero so that in the case where there is no
overlap between the variances of the random effects (or if the variances are
zero), the expected rank equals the actual rank. The ranks are ordered such
that the winners have ranks that are greater than the losers.
The formula for the percentile expected rank is:
100 * (ExpectedRank_i - 0.5) / N_grps
where N_grps is the number of grouping factor levels. The percentile
expected rank can be interpreted as the fraction of levels that score at or
below the given level.
NOTE: expectedRank will only work under conditions that lme4::ranef
will work. One current example of when this is not the case is for
models when there are multiple terms specified per factor (e.g. uncorrelated random
coefficients for the same term, e.g.
lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data = sleepstudy))
Value
A data.frame with the following five columns:
Column 1
The original grouping factor
Column 2
The estimated random effects (from
lme4::ranef(, condVar=TRUE)); name taken from term.
Column 3
The posterior variance of the estimate random effect
(from lme4::ranef(, condVar=TRUE)); named "term"_var.
ER
The expected rank.
pctER
The percentile expected rank.
References
Laird NM and Louis TA. Empirical Bayes Ranking Methods. Journal of
Education Statistics. 1989;14(1)29-46. Available at
http://www.jstor.org/stable/1164724.
Lingsma HF, Steyerberg EW, Eijkemans MJC, et al. Comparing and
ranking hospitals based on outcome: results from The Netherlands Stroke Survey.
QJM: An International Journal of Medicine. 2010;103(2):99-108.
doi:10.1093/qjmed/hcp169
Examples
#For a one-level random intercept model
require(lme4)
m1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
(m1.er <- expectedRank(m1))
#For a one-level random intercept model with multiple random terms
require(lme4)
m2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#ranked by the random slope on Days
(m2.er1 <- expectedRank(m2, term="Days"))
#ranked by the random intercept
(m2.er2 <- expectedRank(m2, term="int"))
## Not run:
#For a two-level model with random intercepts
require(lme4)
m3 <- lmer(y ~ service * dept + (1|s) + (1|d), InstEval)
#Ranked by the random intercept on 's'
(m3.er1 <- expectedRank(m3, groupFctr="s", term="Intercept"))
## End(Not run)