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)
> data(wine)
> 
> #################################
> ## Appropriate evaluation of formulas:
> 
> ## These all work as intended with no warnings or errors:
> fm1 <- clmm(rating ~ contact + (1|judge), data=wine)
> fm1
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -98.80 209.59 228(686) 3.67e-06

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 0.4428   0.6654  
Number of groups:  judge 9 

Coefficients:
contactyes 
       1.3 

Thresholds:
     1|2      2|3      3|4      4|5 
-2.28331  0.04325  1.86062  3.20298 
> fm1 <- clmm("rating ~ contact + (1|judge)", data=wine)
> fm1
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -98.80 209.59 228(686) 3.67e-06

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 0.4428   0.6654  
Number of groups:  judge 9 

Coefficients:
contactyes 
       1.3 

Thresholds:
     1|2      2|3      3|4      4|5 
-2.28331  0.04325  1.86062  3.20298 
> fm1 <- clmm(as.formula("rating ~ contact + (1|judge)"), data=wine)
> fm1 
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -98.80 209.59 228(686) 3.67e-06

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 0.4428   0.6654  
Number of groups:  judge 9 

Coefficients:
contactyes 
       1.3 

Thresholds:
     1|2      2|3      3|4      4|5 
-2.28331  0.04325  1.86062  3.20298 
> fm1 <- clmm(as.formula(rating ~ contact + (1|judge)), data=wine) 
> fm1
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -98.80 209.59 228(686) 3.67e-06

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 0.4428   0.6654  
Number of groups:  judge 9 

Coefficients:
contactyes 
       1.3 

Thresholds:
     1|2      2|3      3|4      4|5 
-2.28331  0.04325  1.86062  3.20298 
> 
> #################################
> 
> ### finding variables in the environment of the formula:
> makeform <- function() {
+   f1 <- as.formula(rating ~ temp + contact + (1|judge))
+   rating <- wine$rating
+   temp <- wine$temp
+   contact <- wine$contact
+   judge <- wine$judge
+   f1
+ }
> ## 'makeform' makes are formula object in the environment of the
> ## function makeform:
> f1 <- makeform()
> f1 # print
rating ~ temp + contact + (1 | judge)
<environment: 0x000001e46128a650>
> class(f1)
[1] "formula"
> ## If we give the data, we can evaluate the model:
> fm1 <- clmm(f1, data=wine)
> ## We can also evaluate the model because the data are available in
> ## the environment associated with the formula:
> fm1 <- clmm(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")
<simpleError: object 'rating' not found>
> 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 + (1|judge))
> (try(fm2 <- clmm(f2), silent=TRUE))
[1] "Error in eval(predvars, data, env) : object 'rating' not found\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in eval(predvars, data, env): object 'rating' not found>
> environment(f2) <- environment(f1)
> fm2 <- clmm(f2)
> 
> #################################
> ## Use of formula-objects
> f <- formula(rating ~ temp + contact + (1|judge))
> m2 <- clmm(f, data = wine)
> summary(m2)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad cond.H 
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05 2.8e+01

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
           Estimate Std. Error z value Pr(>|z|)    
tempwarm     3.0630     0.5954   5.145 2.68e-07 ***
contactyes   1.8349     0.5125   3.580 0.000344 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
    Estimate Std. Error z value
1|2  -1.6237     0.6824  -2.379
2|3   1.5134     0.6038   2.507
3|4   4.2285     0.8090   5.227
4|5   6.0888     0.9725   6.261
> 
> #################################
> ## Other ways to construct formulas:
> set.seed(12345)
> y <- factor(sample(1:4,20,replace=TRUE))
> x <- rnorm(20)
> b <- gl(5, 4, labels=letters[1:5])
> data <- data.frame(y=y, x=x, b=b)
> rm(x, y, b)
> clmm(y ~ x + (1|b), data=data)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: y ~ x + (1 | b)
data:    data

 link  threshold nobs logLik AIC   niter    max.grad
 logit flexible  20   -25.37 60.74 338(288) 7.12e-08

Random effects:
 Groups Name        Variance  Std.Dev. 
 b      (Intercept) 9.736e-09 9.867e-05
Number of groups:  b 5 

Coefficients:
     x 
0.2527 

Thresholds:
     1|2      2|3      3|4 
-2.18234  0.06059  0.92746 
> fit <- clmm(data$y ~ data$x + (1|data$b))
> fit
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: data$y ~ data$x + (1 | data$b)

 link  threshold nobs logLik AIC   niter    max.grad
 logit flexible  20   -25.37 60.74 338(288) 7.12e-08

Random effects:
 Groups Name        Variance  Std.Dev. 
 data$b (Intercept) 9.736e-09 9.867e-05
Number of groups:  data$b 5 

Coefficients:
data$x 
0.2527 

Thresholds:
     1|2      2|3      3|4 
-2.18234  0.06059  0.92746 
> fit <- clmm(data[, 1] ~ data[, 2] + (1|data[, 3]))
> fit
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: data[, 1] ~ data[, 2] + (1 | data[, 3])

 link  threshold nobs logLik AIC   niter    max.grad
 logit flexible  20   -25.37 60.74 338(288) 7.12e-08

Random effects:
 Groups    Name        Variance  Std.Dev. 
 data[, 3] (Intercept) 9.736e-09 9.867e-05
Number of groups:  data[, 3] 5 

Coefficients:
data[, 2] 
   0.2527 

Thresholds:
     1|2      2|3      3|4 
-2.18234  0.06059  0.92746 
> 
> #################################
> ## 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 clmm work with glmulti.
> 
> fun.clmm <- function(formula, data)
+ ### This only works because clmm via eclmm.model.frame is careful to
+ ### evaluate the 'formula' in the parent environment such it is not the
+ ### character "formula" that is attempted evaluated.
+   clmm(formula, data = data)
> 
> fun2.clmm <- function(formula, data, weights, subset) {
+ ### This should be the safe way to ensure evaluation of clmm in the
+ ### right environment.
+   mc <- match.call()
+   mc[[1]] <- as.name("clmm")
+   eval.parent(mc)
+ }
> 
> fun.clmm(rating ~ temp + contact + (1|judge), data=wine) ## works
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    data

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
  tempwarm contactyes 
     3.063      1.835 

Thresholds:
   1|2    2|3    3|4    4|5 
-1.624  1.513  4.229  6.089 
> fun2.clmm(rating ~ temp + contact + (1|judge), data=wine) ## works
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
  tempwarm contactyes 
     3.063      1.835 

Thresholds:
   1|2    2|3    3|4    4|5 
-1.624  1.513  4.229  6.089 
> 
> form1 <- "rating ~ temp + contact + (1|judge)"
> fun.clmm(form1, data=wine) ## works
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    data

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
  tempwarm contactyes 
     3.063      1.835 

Thresholds:
   1|2    2|3    3|4    4|5 
-1.624  1.513  4.229  6.089 
> fun2.clmm(form1, data=wine) ## works
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
  tempwarm contactyes 
     3.063      1.835 

Thresholds:
   1|2    2|3    3|4    4|5 
-1.624  1.513  4.229  6.089 
> 
> form2 <- formula(rating ~ temp + contact + (1|judge))
> fun.clmm(form2, data=wine) ## works
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    data

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
  tempwarm contactyes 
     3.063      1.835 

Thresholds:
   1|2    2|3    3|4    4|5 
-1.624  1.513  4.229  6.089 
> fun2.clmm(form2, data=wine) ## works
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ temp + contact + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad
 logit flexible  72   -81.57 177.13 332(999) 1.03e-05

Random effects:
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131   
Number of groups:  judge 9 

Coefficients:
  tempwarm contactyes 
     3.063      1.835 

Thresholds:
   1|2    2|3    3|4    4|5 
-1.624  1.513  4.229  6.089 
> ## Notice that clmm is not able to get the name of the data (wine)
> ## correct when using fun.clmm.
> 
> #################################
> 
> ##    ## Example 2: using clmm function
> ##    #
> ##    ## Now I want to consider judge as a random effect to account for
> ##    ## grouping structure of data 
> ##    mod2 <- clmm(rating ~ temp + contact + (1|judge), data=wine)
> ##    
> ##    ##Again, I started by using my own code to run all potential models:
> ##    ## put names of all your variables in this vector:
> ##    vl2 <- c("temp", "contact")
> ##    ## generate list of possible combinations of variables:
> ##    combos2 <- NULL
> ##    for(i in 1:length(vl2)) {
> ##      combos2 <- c(combos2, combn(vl2, i, simplify = F))
> ##    }
> ##    ## create formulae and run models one by one, saving them as model1,
> ##    ## model2 etc... 
> ##    for (i in 1:length(combos2)) {
> ##      vs2 <- paste(combos2[[i]], collapse=" + ")
> ##      f2 <- formula(paste("rating ~ ", vs2, "+(1|judge)", sep=""))
> ##      print(f2)
> ##      assign(paste("model", i, sep=""), clmm(f2, data=wine))
> ##    }
> ##    summary(model1) # etc
> ##    summary(model2) # etc
> ##    summary(model3) # etc
> ##    
> ##    models <- vector("list", length(combos2))
> ##    for(i in 1:length(combos2)) {
> ##      vs2 <- paste(combos2[[i]], collapse=" + ")
> ##      f2 <- formula(paste("rating ~ ", vs2, "+(1|judge)", sep=""))
> ##      print(f2)
> ##      models[[i]] <- clmm(f2, data=wine)
> ##      ## assign(paste("model", i, sep=""), clmm(f2, data=wine))
> ##    }
> ##    
> ##    ## Coefficients, AIC and BIC:
> ##    lapply(models, function(m) coef(summary(m)))
> ##    lapply(models, AIC)
> ##    lapply(models, BIC)
> ##    
> ##    ## library(MuMIn)
> ##    ## dd2 <- dredge(mod2) ## does not work
> ##    ## ?dredge
> ##    ## traceback()
> ##    ## mod2$formula
> ##    ## terms(as.formula(formula(mod2)))
> ##    ## 
> ##    ## library(lme4)
> ##    ## fmm1 <- lmer(response ~ temp + contact + (1|judge), data=wine)
> ##    ## fmm1
> ##    ## terms(as.formula(lme4:::formula(fmm1)))
> ##    ## terms(as.formula(formula(fmm1)))
> 
> proc.time()
   user  system elapsed 
   3.20    0.35    3.50