tests/test.clm.model.matrix.R

In ordinal: Regression Models for Ordinal Data

library(ordinal)
## source("test.clm.model.matrix.R")

## library(devtools)
## r2path <- "/Users/rhbc/Documents/Rpackages/ordinal/pkg/ordinal"
## clean_dll(pkg = r2path)
## load_all(r2path)

## Check that get_clmDesign works in standard setting:
fm1 <- clm(rating ~ temp, scale=~contact, nominal=~contact, data=wine)
contr <- c(fm1$contrasts, fm1$S.contrasts, fm1$nom.contrasts)
XX <- ordinal:::get_clmDesign(fm1$model, terms(fm1, "all"), contrasts=contr)
XX2 <- update(fm1, method="design")
(keep <- intersect(names(XX), names(XX2)))
(test <- mapply(function(x, y) isTRUE(all.equal(x, y)),
                XX[keep], XX2[keep]))
stopifnot(all(test))

## Check that get_clmDesign works with singular fit and NAs:
cy <- with(wine, which(temp == "cold" & contact == "yes"))
wine2 <- subset(wine, subset=(!1:nrow(wine) %in% cy))
wine2[c(9, 15, 46), "rating"] <- NA
fm1 <- clm(rating ~ temp, scale=~contact, nominal=~contact,
           data=wine2)
contr <- c(fm1$contrasts, fm1$S.contrasts, fm1$nom.contrasts)
XX <- ordinal:::get_clmDesign(fm1$model, terms(fm1, "all"), contrasts=contr)
XX2 <- update(fm1, method="design")
(keep <- intersect(names(XX), names(XX2)))
(test <- mapply(function(x, y) isTRUE(all.equal(x, y)),
                XX[keep], XX2[keep]))
stopifnot(all(test))

## In this situation update and get_clmRho give the same results:
wine2 <- wine
fm1 <- clm(rating ~ temp + contact, data=wine2) ## OK
rho1 <- ordinal:::get_clmRho.clm(fm1)
l1 <- as.list(rho1)
l2 <- as.list(update(fm1, doFit=FALSE))
(test <- mapply(function(x, y) isTRUE(all.equal(x, y)),
                l1, l2[names(l1)]))
stopifnot(all(test))
## If we modify the data (or other subset, weights, formulae, etc.)
## used in the model call, the results from update no longer correspond
## to the elements of the fitted model object. get_clmRho gets it
## right on the other hand:
wine2[10:13, "rating"] <- NA
l3 <- as.list(ordinal:::get_clmRho.clm(fm1))
l4 <- as.list(update(fm1, doFit=FALSE))
(test <- mapply(function(x, y) isTRUE(all.equal(x, y)),
                l1, l3))
stopifnot(all(test)) ## same
(test <- mapply(function(x, y) isTRUE(all.equal(x, y)),
                l3, l4[names(l3)]))
stopifnot(sum(!test) == 8) ## not all the same anymore!
## In conclusion l1, l2, and l3 are identical. l4 is different.

#################################
## Test that checkContrasts give appropriate warnings:
contr <- c(temp="contr.sum", contact="contr.sum")
fm1 <- clm(rating ~ temp + contact, scale=~contact, data=wine) ## OK
fm1 <- clm(rating ~ temp + contact, scale=~contact, data=wine,
           contrasts=contr) ## OK
fm1 <- clm(rating ~ temp, scale=~contact, data=wine,
           contrasts=contr) ## OK
## These should give warnings:
fm1 <- clm(rating ~ temp, contrasts=c(contact="contr.sum"), data=wine)
fm1 <- clm(rating ~ temp, contrasts=contr, data=wine)
fm1 <- clm(rating ~ 1, scale=~contact, contrasts=c(temp="contr.sum"),
           data=wine)
fm1 <- clm(rating ~ 1, scale=~contact, contrasts=list(temp="contr.sum"),
           data=wine)

fm0 <- clm(rating ~ temp + contact, scale=~contact, data=wine)
ordinal:::checkContrasts(fm0$S.terms, fm0$contrasts)
ordinal:::checkContrasts(fm0$S.terms, fm0$S.contrasts)
ordinal:::checkContrasts(fm0$terms, fm0$contrasts)
ordinal:::checkContrasts(fm0$terms, fm0$S.contrasts)

#################################
## Check that clm and model.matrix respects contrast settings:
options("contrasts" = c("contr.treatment", "contr.poly"))
fm0 <- clm(rating ~ temp + contact, data=wine)
options("contrasts" = c("contr.sum", "contr.poly"))
fm1 <- clm(rating ~ temp + contact, data=wine)
stopifnot(all(model.matrix(fm0)$X[, 2] %in% c(0, 1)))
stopifnot(all(model.matrix(fm1)$X[, 2] %in% c(1, -1)))

#################################
## Check that model.matrix results do not depend on global contrast
## setting:
options("contrasts" = c("contr.sum", "contr.poly"))
fm0 <- clm(rating ~ temp + contact, scale=~contact, data=wine)
MM <- model.matrix(fm0)
options("contrasts" = c("contr.treatment", "contr.poly"))
MM2 <- model.matrix(fm0)
for(x in MM) print(head(x))
for(x in MM2) print(head(x))
stopifnot(all(mapply(all.equal, MM, MM2)))

#################################
## This gave a warning before getContrasts was implemented:
fm0 <- clm(rating ~ temp + contact, scale=~contact, data=wine)
MM <- model.matrix(fm0)
## > fm0 <- clm(rating ~ temp + contact, scale=~contact, data=wine)
## > MM <- model.matrix(fm0)
## Warning message:
## In model.matrix.default(res$S.terms, data = fullmf, contrasts.arg = getContrasts(res$S.terms,  :
##   variable 'temp' is absent, its contrast will be ignored
for(x in MM) print(head(x))