What does predict.glm(, type="terms") actually do?

Question

I am confused with the way predict.glm function in R works. According to the help,

The "terms" option returns a matrix giving the fitted values of each term in the model formula on the linear predictor scale.

Thus, if my model has form f(y) = X*beta, then command

predict(model, X, type='terms')

is expected to produce the same matrix X, multiplied by beta element-wise. For example, if I train the following model

test.data = data.frame(y = c(0,0,0,1,1,1,1,1,1), x=c(1,2,3,1,2,2,3,3,3))
model = glm(y~(x==1)+(x==2), family = 'binomial', data = test.data)

the resulting coefficients are

beta <- model$coef

Design matrix is

X <- model.matrix(y~(x==1)+(x==2), data = test.data)

  (Intercept) x == 1TRUE x == 2TRUE
1           1          1          0
2           1          0          1
3           1          0          0
4           1          1          0
5           1          0          1
6           1          0          1
7           1          0          0
8           1          0          0
9           1          0          0

Then multiplied by coefficients it should look like

pred1 <- t(beta * t(X))

  (Intercept) x == 1TRUE x == 2TRUE
1    1.098612  -1.098612  0.0000000
2    1.098612   0.000000 -0.4054651
3    1.098612   0.000000  0.0000000
4    1.098612  -1.098612  0.0000000
5    1.098612   0.000000 -0.4054651
6    1.098612   0.000000 -0.4054651
7    1.098612   0.000000  0.0000000
8    1.098612   0.000000  0.0000000
9    1.098612   0.000000  0.0000000

However, actual matrix produced by predict.glm seems to be unrelated to this. The following code

pred2 <- predict(model, test.data, type = 'terms')

      x == 1     x == 2
1 -0.8544762  0.1351550
2  0.2441361 -0.2703101
3  0.2441361  0.1351550
4 -0.8544762  0.1351550
5  0.2441361 -0.2703101
6  0.2441361 -0.2703101
7  0.2441361  0.1351550
8  0.2441361  0.1351550
9  0.2441361  0.1351550
attr(,"constant")
[1] 0.7193212

How does one interpret such results?

Answer 1

I have already edited your question, to include "correct" way of getting (raw) model matrix, model coefficients, and your intended term-wise prediction. So your other question on how to get these are already solved. In the following, I shall help you understand predict.glm().

---

predict.glm() (actually, predict.lm()) has applied centring constraints for each model term when doing term-wise prediction.

Initially, you have a model matrix

X <- model.matrix(y~(x==1)+(x==2), data = test.data)

but it is centred, by dropping column means:

avx <- colMeans(X)
X1 <- sweep(X, 2L, avx)

> avx
(Intercept)  x == 1TRUE  x == 2TRUE 
  1.0000000   0.2222222   0.3333333 

> X1
  (Intercept) x == 1TRUE x == 2TRUE
1           0  0.7777778 -0.3333333
2           0 -0.2222222  0.6666667
3           0 -0.2222222 -0.3333333
4           0  0.7777778 -0.3333333
5           0 -0.2222222  0.6666667
6           0 -0.2222222  0.6666667
7           0 -0.2222222 -0.3333333
8           0 -0.2222222 -0.3333333
9           0 -0.2222222 -0.3333333

Then term-wise computation is done using this centred model matrix:

t(beta*t(X1))

  (Intercept) x == 1TRUE x == 2TRUE
1           0 -0.8544762  0.1351550
2           0  0.2441361 -0.2703101
3           0  0.2441361  0.1351550
4           0 -0.8544762  0.1351550
5           0  0.2441361 -0.2703101
6           0  0.2441361 -0.2703101
7           0  0.2441361  0.1351550
8           0  0.2441361  0.1351550
9           0  0.2441361  0.1351550

After centring, different terms are vertically shifted to have zero mean. As a result, intercept will be come 0. No worry, a new intercept is computed, by aggregating shifts of all model terms:

intercept <- as.numeric(crossprod(avx, beta))
# [1] 0.7193212

Now you should have seen what predict.glm(, type = "terms") gives you.

Answer 2

This code will compute uncentered term values.

## Extracted from stats::predict.lm (called by predict.Glm for type='terms')

ucmterms <-
  function (object, newdata, terms = NULL, na.action = na.pass, ...) {
    tt <- terms(object)
    if (missing(newdata) || is.null(newdata)) {
      mm <- X <- model.matrix(object)
      mmDone <- TRUE
      offset <- object$offset
    }
    else {
      Terms <- delete.response(tt)
      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.arg = object$contrasts)
      mmDone <- FALSE
    }
       
    beta <- object$coefficients
    nrp <- num.intercepts(object)   # in rms
    ## If > 1 intercept, remove all but one
    if(nrp > 1) beta <- cbind(beta[1], beta[-(1 : nrp)])
       
    if (!mmDone) {
      mm <- model.matrix(object)
      mmDone <- TRUE
    }
    aa <- attr(mm, "assign")
    ll <- attr(tt, "term.labels")
    hasintercept <- attr(tt, "intercept") > 0L
    if (hasintercept) ll <- c("(Intercept)", ll)
    aaa <- factor(aa, labels = ll)
    asgn <- split(order(aa), aaa)
    if(hasintercept) asgn$"(Intercept)" <- NULL

    nterms <- length(asgn)
    if (nterms > 0) {
      predictor <- matrix(ncol = nterms, nrow = NROW(X))
      dimnames(predictor) <- list(rownames(X), names(asgn))
      for(i in seq.int(1L, nterms, length.out = nterms)) {
        ii <- asgn[[i]]
        predictor[, i] <- X[, ii, drop = FALSE] %*% beta[ii]
      }
      if (!is.null(terms)) predictor <- predictor[, terms, drop = FALSE]
    } else predictor <- matrix(0, NROW(X), 0L)

    attr(predictor, "constant") <- 0

    if (missing(newdata) && ! is.null(na.act <- object$na.action))
      predictor <- napredict(na.act, predictor)

    predictor
}