R Under development (unstable) (2023-12-02 r85657 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit 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))) [1] "y" "y.levels" "X" "offset" [5] "terms" "contrasts" "xlevels" "weights" [9] "S" "S.terms" "S.off" "S.contrasts" [13] "S.xlevels" "NOM" "nom.terms" "nom.contrasts" [17] "nom.xlevels" > (test <- mapply(function(x, y) isTRUE(all.equal(x, y)), + XX[keep], XX2[keep])) y y.levels X offset terms TRUE TRUE TRUE TRUE TRUE contrasts xlevels weights S S.terms TRUE TRUE TRUE TRUE TRUE S.off S.contrasts S.xlevels NOM nom.terms TRUE TRUE TRUE TRUE TRUE nom.contrasts nom.xlevels TRUE TRUE > 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) Warning message: (1) Hessian is numerically singular: parameters are not uniquely determined In addition: Absolute convergence criterion was met, but relative criterion was not met > 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))) [1] "y" "y.levels" "X" "offset" [5] "terms" "contrasts" "xlevels" "na.action" [9] "weights" "S" "S.terms" "S.off" [13] "S.contrasts" "S.xlevels" "NOM" "nom.terms" [17] "nom.contrasts" "nom.xlevels" > (test <- mapply(function(x, y) isTRUE(all.equal(x, y)), + XX[keep], XX2[keep])) y y.levels X offset terms TRUE TRUE TRUE TRUE TRUE contrasts xlevels na.action weights S TRUE TRUE TRUE TRUE TRUE S.terms S.off S.contrasts S.xlevels NOM TRUE TRUE TRUE TRUE TRUE nom.terms nom.contrasts nom.xlevels TRUE TRUE TRUE > 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)])) nlambda link gfun dfun pfun par clm.hess clm.grad TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE clm.nll wts fitted has.scale sigma k Soff S TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE n.psi o2 o1 B2 B1 TRUE TRUE TRUE TRUE TRUE > 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)) nlambda link gfun dfun pfun par clm.hess clm.grad TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE clm.nll wts fitted has.scale sigma k Soff S TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE n.psi o2 o1 B2 B1 TRUE TRUE TRUE TRUE TRUE > stopifnot(all(test)) ## same > (test <- mapply(function(x, y) isTRUE(all.equal(x, y)), + l3, l4[names(l3)])) nlambda link gfun dfun pfun par clm.hess clm.grad TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE clm.nll wts fitted has.scale sigma k Soff S TRUE FALSE FALSE TRUE FALSE TRUE FALSE TRUE n.psi o2 o1 B2 B1 TRUE FALSE FALSE FALSE FALSE > 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 Warning messages: 1: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored 2: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored > fm1 <- clm(rating ~ temp, scale=~contact, data=wine, + contrasts=contr) ## OK Warning messages: 1: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored 2: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored > ## These should give warnings: > fm1 <- clm(rating ~ temp, contrasts=c(contact="contr.sum"), data=wine) Warning messages: 1: variable 'contact' is absent: its contrasts will be ignored 2: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored > fm1 <- clm(rating ~ temp, contrasts=contr, data=wine) Warning messages: 1: variable 'contact' is absent: its contrasts will be ignored 2: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored > fm1 <- clm(rating ~ 1, scale=~contact, contrasts=c(temp="contr.sum"), + data=wine) Warning messages: 1: variable 'temp' is absent: its contrasts will be ignored 2: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored 3: In model.matrix.default(terms, data = fullmf, contrasts.arg = getContrasts(terms, : non-list contrasts argument ignored > fm1 <- clm(rating ~ 1, scale=~contact, contrasts=list(temp="contr.sum"), + data=wine) Warning message: variable 'temp' is absent: its contrasts will be ignored > > fm0 <- clm(rating ~ temp + contact, scale=~contact, data=wine) > ordinal:::checkContrasts(fm0$S.terms, fm0$contrasts) Warning message: variable 'temp' is absent: its contrasts will be ignored > 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)) (Intercept) temp1 contact1 1 1 1 1 2 1 1 1 3 1 1 -1 4 1 1 -1 5 1 -1 1 6 1 -1 1 (Intercept) contact1 1 1 1 2 1 1 3 1 -1 4 1 -1 5 1 1 6 1 1 > for(x in MM2) print(head(x)) (Intercept) temp1 contact1 1 1 1 1 2 1 1 1 3 1 1 -1 4 1 1 -1 5 1 -1 1 6 1 -1 1 (Intercept) contact1 1 1 1 2 1 1 3 1 -1 4 1 -1 5 1 1 6 1 1 > 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)) (Intercept) tempwarm contactyes 1 1 0 0 2 1 0 0 3 1 0 1 4 1 0 1 5 1 1 0 6 1 1 0 (Intercept) contactyes 1 1 0 2 1 0 3 1 1 4 1 1 5 1 0 6 1 0 > > > proc.time() user system elapsed 1.70 0.23 1.92