For developers: Extending emmeans

The recover_data method

It is usually an easy matter to write a recover_data method. Look at the one for lm objects:

emmeans:::recover_data.lm

## function (object, ...) 
## {
##     fcall = object$call
##     recover_data(fcall, delete.response(terms(object)), object$na.action, 
##         ...)
## }
## <bytecode: 0x000000001f978fc8>
## <environment: namespace:emmeans>

Note that all it does is obtain the call component and call the method for class call, with additional arguments for its terms component and na.action. It happens that we can access these attributes in exactly the same way as for lm objects; so:

recover_data.lqs = emmeans:::recover_data.lm

Let’s test it:

rec.fake = recover_data(fake.lts)
head(rec.fake)

##    A  B
## 1 a1 b1
## 2 a1 b1
## 3 a1 b1
## 4 a1 b1
## 5 a1 b1
## 6 a2 b1

Our recovered data excludes the response variable y (owing to the delete.response call), and this is fine.

The emm_basis method

The emm_basis method has four required arguments:

args(emmeans:::emm_basis.lm)

## function (object, trms, xlev, grid, ...) 
## NULL

These are, respectively, the model object, its terms component (at least for the right-hand side of the model), a list of levels of the factors, and the grid of predictor combinations that specify the reference grid.

The function must obtain six things and return them in a named list. They are the matrix X of linear functions for each point in the reference grid, the regression coefficients bhat; the variance-covariance matrix V; a matrix nbasis for non-estimable functions; a function dffun(k,dfargs) for computing degrees of freedom for the linear function sum(k*bhat); and a list dfargs of arguments to pass to dffun.

To write your own emm_basis function, examining some of the existing methods can help; but the best resource is the predict method for the object in question, looking carefully to see what it does to predict values for a new set of predictors (e.g., newdata in predict.lm). Following this advice, let’s take a look at it:

MASS:::predict.lqs

## function (object, newdata, na.action = na.pass, ...) 
## {
##     if (missing(newdata)) 
##         return(fitted(object))
##     Terms <- delete.response(terms(object))
##     m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels)
##     if (!is.null(cl <- attr(Terms, "dataClasses"))) 
##         .checkMFClasses(cl, m)
##     X <- model.matrix(Terms, m, contrasts = object$contrasts)
##     drop(X %*% object$coefficients)
## }
## <bytecode: 0x0000000025225b08>
## <environment: namespace:MASS>

emm_basis.lqs = function(object, trms, xlev, grid, ...) { 
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
    X = model.matrix(trms, m, contrasts.arg = object$contrasts) 
    bhat = coef(object) 
    Xmat = model.matrix(trms, data=object$model)                      # 5
    V = rev(object$scale)[1]^2 * solve(t(Xmat) %*% Xmat)
    nbasis = matrix(NA) 
    dfargs = list(df = nrow(Xmat) - ncol(Xmat))
    dffun = function(k, dfargs) dfargs$df
    list(X = X, bhat = bhat, nbasis = nbasis, V = V,                  #10
         dffun = dffun, dfargs = dfargs)
}

emmeans(fake.lts, ~ B | A)

## A = a1:
##  B  emmean    SE df lower.CL upper.CL
##  b1   11.9 0.228 24     11.4     12.3
##  b2   23.1 0.228 24     22.6     23.6
##  b3   17.8 0.228 24     17.3     18.2
## 
## A = a2:
##  B  emmean    SE df lower.CL upper.CL
##  b1   13.9 0.228 24     13.4     14.4
##  b2   24.1 0.228 24     23.6     24.5
##  b3   20.5 0.228 24     20.0     21.0
## 
## Confidence level used: 0.95

Dissecting emm_basis.lqs

Let’s go through the listing of this method, line-by-line:

• Lines 2–3: Construct the linear functions, X. This is a pretty standard two-step process: First obtain a model frame, m, for the grid of predictors, then pass it as data to model.matrix to create the associated design matrix. As promised, this code is essentially identical to what you find in predict.lqs.

• Line 4: Obtain the coefficients, bhat. Most model objects have a coef method.

• Lines 5–6: Obtain the covariance matrix, V, of bhat. In many models, this can be obtained using the object’s vcov method. But not in this case. Instead, I cobbled one together using the inverse of the X’X matrix as in ordinary regression, and the variance estimate found in the last element of the scale element of the object. This probably under-estimates the variances and distorts the covariances, because robust estimators have some efficiency loss.

• Line 7: Compute the basis for non-estimable functions. This applies only when there is a possibility of rank deficiency in the model. But lqs methods don’t allow rank deficiencies, so it we have fitted such a model, we can be sure that all linear functions are estimable; we signal that by setting nbasis equal to a 1 x 1 matrix of NA. If rank deficiency were possible, the estimability package (which is required by emmeans) provides a nonest.basis function that makes this fairly painless—I would have coded nbasis = estimability::nonest.basis(Xmat).

  There some subtleties you need to know regarding estimability. Suppose the model is rank-deficient, so that the design matrix X has p columns but rank r < p. In that case, bhat should be of length p (not r), and there should be p - r elements equal to NA, corresponding to columns of X that were excluded from the fit. Also, X should have all p columns. In other words, do not alter or throw-out columns of X or their corresponding elements of bhat—even those with NA coefficients—as they are essential for assessing estimability. V should be r x r, however—the covariance matrix for the non-excluded predictors.

• Lines 8–9: Obtain dffun and dfargs. This is a little awkward because it is designed to allow support for mixed models, where approximate methods may be used to obtain degrees of freedom. The function dffun is expected to have two arguments: k, the vector of coefficients of bhat, and dfargs, a list containing any additional arguments. In this case (and in many other models), the degrees of freedom are the same regardless of k. We put the required degrees of freedom in dfargs and write dffun so that it simply returns that value. (Note: If asymptotic tests and CIs are desired, return Inf degrees of freedom.)

• Line 10: Return these results in a named list.

Back to Contents

The recover_data method

Note that coding is a predictor transformation, not a response transformation (we could have that, too, as it’s already supported by the emmeans infrastructure). So, to handle the "decode" mode, we will need to actually decode the predictors used to construct he reference grid. That means we need to make recover_data a lot fancier! Here it is:

recover_data.rsm = function(object, data, mode = c("asis", "coded", "decoded"), ...) {
    mode = match.arg(mode)
    cod = rsm::codings(object)
    fcall = object$call
    if(is.null(data))                                                 # 5
        data = emmeans::recover_data(fcall, 
                   delete.response(terms(object)), object$na.action, ...)
    if (!is.null(cod) && (mode == "decoded")) {
        pred = cpred = attr(data, "predictors")
        trms = attr(data, "terms")                                    #10
        data = rsm::decode.data(rsm::as.coded.data(data, formulas = cod))
        for (form in cod) {
            vn = all.vars(form)
            if (!is.na(idx <- grep(vn[1], pred))) { 
                pred[idx] = vn[2]                                     #15
                cpred = setdiff(cpred, vn[1])
            }
        }
        attr(data, "predictors") = pred
        new.trms = update(trms, reformulate(c("1", cpred)))           #20
        attr(new.trms, "orig") = trms
        attr(data, "terms") = new.trms
        attr(data, "misc") = cod
    }
    data
}

Lines 2–7 ensure that mode is legal, retrieves the codings from the object, and obtain the results we would get from recover_data had it been an lm object. If mode is not "decoded", or if no codings were used, that’s all we need. Otherwise, we need to return the decoded data. However, it isn’t quite that simple, because the model equation is still defined on the coded scale. Rather than to try to translate the model coefficients and covariance matrix to the decoded scale, we elected to remember what we will need to do later to put things back on the coded scale. In lines 9–10, we retrieve the attributes of the recovered data that provide the predictor names and terms object on the coded scale. In line 11, we replace the recovered data with the decoded data.

By the way, the codings comprise a list of formulas with the coded name on the left and the original variable name on the right. It is possible that only some of the predictors are coded (for example, blocking factors will not be). In the for loop in lines 12–18, the coded predictor names are replaced with their decoded names. For technical reasons to be discussed later, we also remove these coded predictor names from a copy, cpred, of the list of all predictors in the coded model. In line 19, the "predictors" attribute of data is replaced with the modified version.

Now, there is a nasty technicality. The ref_grid function in emmeans has a few lines of code after recover_data is called that determine if any terms in the model convert covariates to factors or vice versa; and this code uses the model formula. That formula involves variables on the coded scale, and those variables are no longer present in the data, so an error will occur if it tries to access them. Luckily, if we simply take those terms out of the formula, it won’t hurt because those coded predictors would not have been converted in that way. So in line 20, we update trms with a simpler model with the coded variables excluded (the intercept is explicitly included to ensure there will be a right-hand side even is cpred is empty). We save that as the terms attribute, and the original terms as a new "orig" attribute to be retrieved later. The data object, modified or not, is returned. If data have been decoded, ref_grid will construct its grid using decoded variables.

In line 23, we save the codings as the "misc" attribute, to be accessed later by emm_basis().

The emm_basis method

Now comes the emm_basis method that will be called after the grid is defined. It is listed below:

emm_basis.rsm = function(object, trms, xlev, grid, 
                         mode = c("asis", "coded", "decoded"), misc, ...) {
    mode = match.arg(mode)
    cod = misc
    if(!is.null(cod) && mode == "decoded") {                          # 5
        grid = rsm::coded.data(grid, formulas = cod)
        trms = attr(trms, "orig")
    }
    
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)     #10
    X = model.matrix(trms, m, contrasts.arg = object$contrasts)
    bhat = as.numeric(object$coefficients) 
    V = emmeans::.my.vcov(object, ...)
    
    if (sum(is.na(bhat)) > 0)                                         #15
        nbasis = estimability::nonest.basis(object$qr)
    else
        nbasis = estimability::all.estble
    dfargs = list(df = object$df.residual)
    dffun = function(k, dfargs) dfargs$df                             #20

    list(X = X, bhat = bhat, nbasis = nbasis, V = V, 
         dffun = dffun, dfargs = dfargs, misc = list())
}

This is much simpler. The coding formulas are obtained from misc (line 4) so that we don’t have to re-obtain them from the object. All we have to do is determine if decoding was done (line 5); and, if so, convert the grid back to the coded scale (line 6) and recover the original terms attribute (line 7). The rest is borrowed directly from the emm_basis.lm method in emmeans. Note that line 13 uses one of the exported functions we described in the preceding section. Lines 15–18 use functions from the estimability package to handle the possibility that the model is rank-deficient.

A demonstration

Here’s a demonstration of this rsm support. The standard example for rsm fits a second-order model CR.rs2 to a dataset organized in two blocks and with two coded predictors.

library("rsm")
example("rsm")   ### (output is not shown) ###

First, let’s look at some results on the coded scale—which are the same as for an ordinary lm object.

emmeans(CR.rs2, ~ x1 * x2, mode = "coded", 
        at = list(x1 = c(-1, 0, 1), x2 = c(-2, 2)))

##  x1 x2 emmean    SE df lower.CL upper.CL
##  -1 -2   75.0 0.298  7     74.3     75.7
##   0 -2   77.0 0.240  7     76.4     77.5
##   1 -2   76.4 0.298  7     75.6     77.1
##  -1  2   76.8 0.298  7     76.1     77.5
##   0  2   79.3 0.240  7     78.7     79.9
##   1  2   79.2 0.298  7     78.5     79.9
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95

Now, the coded variables x1 and x2 are derived from these coding formulas for predictors Time and Temp:

codings(CR.rs1)

## $x1
## x1 ~ (Time - 85)/5
## 
## $x2
## x2 ~ (Temp - 175)/5

Thus, for example, a coded value of x1 = 1 corresponds to a time of 85 + 1 x 5 = 90. Here are some results working with decoded predictors. Note that the at list must now be given in terms of Time and Temp:

emmeans(CR.rs2, ~ Time * Temp, mode = "decoded", 
        at = list(Time = c(80, 85, 90), Temp = c(165, 185)))

##  Time Temp emmean    SE df lower.CL upper.CL
##    80  165   75.0 0.298  7     74.3     75.7
##    85  165   77.0 0.240  7     76.4     77.5
##    90  165   76.4 0.298  7     75.6     77.1
##    80  185   76.8 0.298  7     76.1     77.5
##    85  185   79.3 0.240  7     78.7     79.9
##    90  185   79.2 0.298  7     78.5     79.9
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95

Since the supplied settings are the same on the decoded scale as were used on the coded scale, the EMMs are identical to those in the previous output.