Predict Method for CLM fits

Obtains predictions from a cumulative link (mixed) model.

Usage

# S3 method for class 'clm2'
predict(object, newdata, ...)

<!-- %% \method{predict}{clmm}(object, newdata, ...) -->

Arguments

object: a fitted object of class inheriting from clm2 including clmm2 objects.
newdata: optionally, a data frame in which to look for variables with which to predict. Observe that the response variable should also be present.
...: further arguments passed to or from other methods.

Details

This method does not duplicate the behavior of predict.polr in package MASS which produces a matrix instead of a vector of predictions. The behavior of predict.polr can be mimiced as shown in the examples.

If newdata is not supplied, the fitted values are obtained. For clmm2 fits this means predictions that are controlled for the observed value of the random effects. If the predictions for a random effect of zero, i.e. an average 'subject', are wanted, the same data used to fit the model should be supplied in the newdata argument. For clm2 fits those two sets of predictions are identical.

Value

A vector of predicted probabilities.

Author

Rune Haubo B Christensen

Examples

options(contrasts = c("contr.treatment", "contr.poly"))

## More manageable data set for less voluminous printing:
(tab26 <- with(soup, table("Product" = PROD, "Response" = SURENESS)))
#>        Response
#> Product   1   2   3   4   5   6
#>    Ref  132 161  65  41 121 219
#>    Test  96  99  50  57 156 650
dimnames(tab26)[[2]] <- c("Sure", "Not Sure", "Guess", "Guess", "Not Sure", "Sure")
dat26 <- expand.grid(sureness = as.factor(1:6), prod = c("Ref", "Test"))
dat26$wghts <- c(t(tab26))
dat26
#>    sureness prod wghts
#> 1         1  Ref   132
#> 2         2  Ref   161
#> 3         3  Ref    65
#> 4         4  Ref    41
#> 5         5  Ref   121
#> 6         6  Ref   219
#> 7         1 Test    96
#> 8         2 Test    99
#> 9         3 Test    50
#> 10        4 Test    57
#> 11        5 Test   156
#> 12        6 Test   650

m1 <- clm2(sureness ~ prod, scale = ~prod, data = dat26,
          weights = wghts, link = "logistic")
predict(m1)
#>  [1] 0.18373313 0.20510835 0.08438208 0.06752717 0.16667487 0.29257440
#>  [7] 0.08288757 0.09840656 0.04839246 0.04385503 0.13834829 0.58811009

mN1 <-  clm2(sureness ~ 1, nominal = ~prod, data = dat26,
            weights = wghts)
predict(mN1)
#>  [1] 0.17861975 0.21786198 0.08795670 0.05548038 0.16373478 0.29634641
#>  [7] 0.08664260 0.08935018 0.04512635 0.05144404 0.14079422 0.58664260

predict(update(m1, scale = ~.-prod))
#>  [1] 0.19702319 0.19835927 0.07932257 0.06291764 0.15503077 0.30734657
#>  [7] 0.07246407 0.09986942 0.05111144 0.04674472 0.14759301 0.58221735


#################################
## Mimicing the behavior of predict.polr:
if(require(MASS)) {
    ## Fit model from polr example:
    fm1 <- clm2(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
    predict(fm1)

    set.seed(123)
    nlev <- 3
    y <- gl(nlev, 5)
    x <- as.numeric(y) + rnorm(15)
    fm.clm <- clm2(y ~ x)
    fm.polr <- polr(y ~ x)

    ## The equivalent of predict.polr(object, type = "probs"):
    (pmat.polr <- predict(fm.polr, type = "probs"))
    ndat <- expand.grid(y = gl(nlev,1), x = x)
    (pmat.clm <- matrix(predict(fm.clm, newdata = ndat), ncol=nlev,
                        byrow = TRUE))
    all.equal(c(pmat.clm), c(pmat.polr), tol = 1e-5) # TRUE

    ## The equivalent of predict.polr(object, type = "class"):
    (class.polr <- predict(fm.polr))
    (class.clm <- factor(apply(pmat.clm, 1, which.max)))
    all.equal(class.clm, class.polr) ## TRUE
}
#> [1] TRUE