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) > ## library(devtools) > ## r2path <- "/Users/rhbc/Documents/Rpackages/ordinal/pkg/ordinal" > ## clean_dll(pkg = r2path) > ## load_all(r2path) > > ################################# > ## Appropriate evaluation of formulas: > > ## These fail and give appropriate error messages: > ## fm1 <- clm(rating ~ contact, scale=temp, data=wine) > ## fm1 <- clm(rating ~ contact, scale=~Temp, data=wine) > ## fm1 <- clm(rating ~ contact, scale="temp", data=wine) > ## sca <- "temp" > ## fm1 <- clm(rating ~ contact, scale=sca, data=wine) > ## sca <- as.formula(sca) > ## sca <- as.formula(temp) > ## sca <- with(wine, as.formula(temp)) > > ## These all work as intended with no warnings or errors: > fm1 <- clm(rating ~ contact, scale="~temp", data=wine) > fm1 <- clm(rating ~ contact, scale=~temp, data=wine) > sca <- "~temp" > fm1 <- clm(rating ~ contact, scale=sca, data=wine) > sca <- as.formula("~temp") > fm1 <- clm(rating ~ contact, scale=sca, data=wine) > fm1 <- clm(rating ~ contact, scale=as.formula(~temp), data=wine) > fm1 <- clm(rating ~ contact, scale=as.formula("~temp"), data=wine) > > ################################# > ## can evaluate if 'formula' is a character: > f <- "rating ~ contact + temp" > clm(f, data=wine) formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01 Coefficients: contactyes tempwarm 1.528 2.503 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > clm(as.formula(f), data=wine) formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01 Coefficients: contactyes tempwarm 1.528 2.503 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > > ################################# > > ### finding variables in the environment of the formula: > makeform <- function() { + f1 <- as.formula(rating ~ temp + contact) + rating <- wine$rating + temp <- wine$temp + contact <- wine$contact + f1 + } > ## 'makeform' makes are formula object in the environment of the > ## function makeform: > f1 <- makeform() > f1 # print rating ~ temp + contact > class(f1) [1] "formula" > ## If we give the data, we can evaluate the model: > fm1 <- clm(f1, data=wine) > ## We can also evaluate the model because the data are available in > ## the environment associated with the formula: > fm1 <- clm(f1) > ## For instance, the 'rating' variable is not found in the Global > ## environment; we have to evaluate the 'name' of 'rating' in the > ## appropriate environment: > (try(rating, silent=TRUE)) [1] "Error : object 'rating' not found\n" attr(,"class") [1] "try-error" attr(,"condition") > eval(as.name("rating"), envir=environment(f1)) [1] 2 3 3 4 4 4 5 5 1 2 1 3 2 3 5 4 2 3 3 2 5 5 4 4 3 2 3 2 3 2 5 3 2 3 4 3 3 3 [39] 3 3 3 2 3 2 2 4 5 4 1 1 2 2 2 3 2 3 2 2 2 3 3 3 3 4 1 2 3 2 3 2 4 4 Levels: 1 < 2 < 3 < 4 < 5 > ## If instead we generate the formula in the Global environment where > ## the variables are not found, we cannot evaluate the model: > f2 <- as.formula(rating ~ temp + contact) > (try(fm2 <- clm(f2), silent=TRUE)) [1] "Error in eval(predvars, data, env) : object 'rating' not found\n" attr(,"class") [1] "try-error" attr(,"condition") > environment(f2) <- environment(f1) > fm2 <- clm(f2) > > > ################################# > ## Use of formula-objects in location, scale and nominal: > ## Bug-report from Lluís Marco Almagro > ## 5 May 2010 17:58 > f <- formula(rating ~ temp) > fs <- formula( ~ contact) > m2 <- clm(f, scale = fs, data = wine) > summary(m2) formula: rating ~ temp scale: ~contact data: wine link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -91.92 195.83 8(0) 1.57e-09 5.6e+01 Coefficients: Estimate Std. Error z value Pr(>|z|) tempwarm 2.1690 0.5527 3.924 8.69e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 log-scale coefficients: Estimate Std. Error z value Pr(>|z|) contactyes -0.1222 0.2775 -0.44 0.66 Threshold coefficients: Estimate Std. Error z value 1|2 -1.8957 0.4756 -3.986 2|3 0.3668 0.3442 1.066 3|4 2.2483 0.5932 3.790 4|5 3.5393 0.8336 4.246 > > ################################# > ## Other ways to construct formulas: > set.seed(12345) > y <- factor(sample(1:4,20,replace=TRUE)) > x <- rnorm(20) > data <- data.frame(y=y,x=x) > rm(x, y) > fit <- clm(data$y ~ data$x) > fit formula: data$y ~ data$x link threshold nobs logLik AIC niter max.grad cond.H logit flexible 20 -25.37 58.74 5(0) 2.68e-09 8.4e+00 Coefficients: data$x 0.2527 Threshold coefficients: 1|2 2|3 3|4 -2.18234 0.06059 0.92746 > fit <- clm(data[,1] ~ data[,2]) > fit formula: data[, 1] ~ data[, 2] link threshold nobs logLik AIC niter max.grad cond.H logit flexible 20 -25.37 58.74 5(0) 2.68e-09 8.4e+00 Coefficients: data[, 2] 0.2527 Threshold coefficients: 1|2 2|3 3|4 -2.18234 0.06059 0.92746 > ## This previously failed, but now works: > fit <- clm(data$y ~ data$x, ~data$x) > fit formula: data$y ~ data$x scale: ~data$x link threshold nobs logLik AIC niter max.grad cond.H logit flexible 20 -25.32 60.64 7(0) 1.86e-09 9.7e+00 Coefficients: data$x 0.2477 log-scale coefficients: data$x -0.1136 Threshold coefficients: 1|2 2|3 3|4 -2.16950 0.07039 0.91381 > > ################################# > ## Evaluation within other functions: > ## date: January 18th 2012. > ## > ## The problem was raised by Stefan Herzog (stefan.herzog@unibas.ch) > ## January 12th 2012 in trying to make clm work with glmulti. > > fun.clm <- function(formula, data) + ### This only works because clm via eclm.model.frame is careful to + ### evaluate the 'formula' in the parent environment such it is not the + ### character "formula" that is attempted evaluated. + clm(formula, data = data) > > fun2.clm <- function(formula, data, weights, subset) { + ### This should be the safe way to ensure evaluation of clm in the + ### right environment. + mc <- match.call() + mc[[1]] <- as.name("clm") + eval.parent(mc) + } > > fun.clm(rating ~ temp + contact, data=wine) ## works formula: rating ~ temp + contact data: data link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 Coefficients: tempwarm contactyes 2.503 1.528 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > fun2.clm(rating ~ temp + contact, data=wine) ## works formula: rating ~ temp + contact data: wine link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 Coefficients: tempwarm contactyes 2.503 1.528 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > > form1 <- "rating ~ temp + contact" > fun.clm(form1, data=wine) ## works formula: rating ~ temp + contact data: data link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 Coefficients: tempwarm contactyes 2.503 1.528 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > fun2.clm(form1, data=wine) ## works formula: rating ~ temp + contact data: wine link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 Coefficients: tempwarm contactyes 2.503 1.528 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > > form2 <- formula(rating ~ temp + contact) > fun.clm(form2, data=wine) ## works formula: rating ~ temp + contact data: data link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 Coefficients: tempwarm contactyes 2.503 1.528 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > fun2.clm(form2, data=wine) ## works formula: rating ~ temp + contact data: wine link threshold nobs logLik AIC niter max.grad cond.H logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 Coefficients: tempwarm contactyes 2.503 1.528 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.344 1.251 3.467 5.006 > ## Notice that clm is not able to get the name of the data (wine) > ## correct when using fun.clm. > > ################################# > ## Evaluation of long formulas: no line breaking in getFullForm: > data(soup, package="ordinal") > > rhs <- paste(names(soup)[c(3, 5:12)], collapse=" + ") > Location <- as.formula(paste("SURENESS ~ ", rhs, sep=" ")) > Scale <- as.formula("~ PROD") > > fm5 <- clm(Location, scale=Scale, data=soup) > summary(fm5) formula: SURENESS ~ PRODID + DAY + SOUPTYPE + SOUPFREQ + COLD + EASY + GENDER + AGEGROUP + LOCATION scale: ~PROD data: soup link threshold nobs logLik AIC niter max.grad cond.H logit flexible 1847 -2651.65 5367.30 8(1) 4.65e-07 1.1e+04 Coefficients: Estimate Std. Error z value Pr(>|z|) PRODID2 1.09324 0.14790 7.392 1.45e-13 *** PRODID3 1.54233 0.21252 7.257 3.95e-13 *** PRODID4 0.95602 0.18180 5.258 1.45e-07 *** PRODID5 1.52404 0.20311 7.504 6.21e-14 *** PRODID6 1.83964 0.21750 8.458 < 2e-16 *** DAY2 -0.26999 0.10280 -2.626 0.00863 ** SOUPTYPECanned -0.31340 0.11745 -2.668 0.00762 ** SOUPTYPEDry-mix 0.22429 0.14295 1.569 0.11665 SOUPFREQ1-4/month -0.06459 0.10551 -0.612 0.54046 SOUPFREQ<1/month 0.03664 0.21310 0.172 0.86350 COLDYes 0.31373 0.14894 2.106 0.03517 * EASY2 0.16689 0.55063 0.303 0.76182 EASY3 0.34765 0.50432 0.689 0.49061 EASY4 0.30008 0.48283 0.621 0.53428 EASY5 0.23492 0.47686 0.493 0.62226 EASY6 -0.03205 0.46993 -0.068 0.94562 EASY7 0.02823 0.46963 0.060 0.95208 EASY8 -0.04600 0.47066 -0.098 0.92215 EASY9 -0.04084 0.49595 -0.082 0.93437 EASY10 1.03641 0.60080 1.725 0.08452 . GENDERFemale -0.05092 0.10738 -0.474 0.63534 AGEGROUP31-40 0.02098 0.14822 0.142 0.88745 AGEGROUP41-50 0.23712 0.15824 1.498 0.13401 AGEGROUP51-65 -0.13791 0.14393 -0.958 0.33798 LOCATIONRegion 2 -0.09641 0.13225 -0.729 0.46599 LOCATIONRegion 3 0.06062 0.11981 0.506 0.61290 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 log-scale coefficients: Estimate Std. Error z value Pr(>|z|) PRODTest 0.14545 0.06621 2.197 0.028 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Threshold coefficients: Estimate Std. Error z value 1|2 -1.64700 0.49982 -3.295 2|3 -0.58840 0.49681 -1.184 3|4 -0.23496 0.49669 -0.473 4|5 0.04394 0.49679 0.088 5|6 0.78551 0.49869 1.575 > > ################################# > ## Check that "."-notation works in formula: > ## December 25th 2014, RHBC > data(wine) > wine2 <- wine[c("rating", "contact", "temp")] > str(wine2) 'data.frame': 72 obs. of 3 variables: $ rating : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 2 3 3 4 4 4 5 5 1 2 ... $ contact: Factor w/ 2 levels "no","yes": 1 1 2 2 1 1 2 2 1 1 ... $ temp : Factor w/ 2 levels "cold","warm": 1 1 1 1 2 2 2 2 1 1 ... > fm0 <- clm(rating ~ ., data=wine2) > fm1 <- clm(rating ~ contact + temp, data=wine2) > keep <- c("coefficients", "logLik", "info") > fun <- function(x, y) stopifnot(isTRUE(all.equal(x, y))) > mapply(fun, fm0[keep], fm1[keep]) $coefficients NULL $logLik NULL $info NULL > ################################# > > proc.time() user system elapsed 1.87 0.26 2.10