tests/test.clm.model.matrix.R

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))

Try the ordinal package in your browser

Any scripts or data that you put into this service are public.

ordinal documentation built on Sept. 11, 2024, 7:44 p.m.