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)
> options(contrasts = c("contr.treatment", "contr.poly"))
> ## library(devtools)
> ## r2path <- "/Users/rhbc/Documents/Rpackages/ordinal/pkg/ordinal"
> ## clean_dll(pkg = r2path)
> ## load_all(r2path)
> 
> ## one zero weight:
> data(wine, package="ordinal")
> wts <- rep(1, nrow(wine))
> wine$rating
 [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
> wts[1] <- 0
> fm1 <- clm(rating ~ contact + temp, data=wine, weights=wts)
> fm1
formula: rating ~ contact + temp
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H 
 logit flexible  71   -85.92 183.84 6(0)  3.43e-12 2.7e+01

Coefficients:
contactyes   tempwarm 
     1.511      2.482 

Threshold coefficients:
   1|2    2|3    3|4    4|5 
-1.295  1.220  3.437  4.975 
> fm1$n ## 72
[1] 72
> fm1$nobs ## 71
[1] 71
> confint(fm1)
               2.5 %   97.5 %
contactyes 0.5951066 2.478924
tempwarm   1.4856829 3.577621
> plot(profile(fm1))
> plot(slice(fm1), 5)
> convergence(fm1)
 nobs logLik niter max.grad cond.H  logLik.Error
 71   -85.92 6(0)  3.43e-12 2.7e+01 <1e-10      

           Estimate Std.Err  Gradient     Error Cor.Dec Sig.Dig
1|2          -1.295  0.5220  7.57e-13  7.84e-14      12      13
2|3           1.220  0.4436  2.65e-12 -1.27e-13      12      13
3|4           3.437  0.6014 -3.43e-12 -7.17e-13      11      12
4|5           4.975  0.7341 -6.06e-14 -7.01e-13      11      12
contactyes    1.511  0.4785  1.53e-13 -1.95e-13      12      13
tempwarm      2.482  0.5303 -4.07e-13 -4.95e-13      12      13

Eigen values of Hessian:
21.7419 18.5476 10.4115  5.2426  4.0755  0.8013 

Convergence message from clm:
(0) successful convergence 
In addition: Absolute and relative convergence criteria were met 

> drop1(fm1, test="Chi")
Single term deletions

Model:
rating ~ contact + temp
        Df    AIC    LRT  Pr(>Chi)    
<none>     183.84                     
contact  1 192.54 10.698  0.001072 ** 
temp     1 208.07 26.234 3.024e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> add1(fm1, scope=~.^2, test="Chi")
Single term additions

Model:
rating ~ contact + temp
             Df    AIC     LRT Pr(>Chi)
<none>          183.84                 
contact:temp  1 185.66 0.18256   0.6692
> ## clm_anova(fm1)
> pred <- predict(fm1, newdata=wine) ## OK
> step.fm1 <- step(fm1, trace=0)
> fitted(fm1)
 [1] 0.19683278 0.44514183 0.09677403 0.20147980 0.20147980 0.27261558
 [7] 0.27261558 0.21497073 0.55705124 0.05699435 0.44514183 0.19814731
[13] 0.50160837 0.27261558 0.36288756 0.55705124 0.19683278 0.44514183
[19] 0.37073110 0.07640020 0.07640020 0.36288756 0.36288756 0.19683278
[25] 0.55705124 0.44514183 0.37073110 0.50160837 0.19814731 0.27261558
[31] 0.30572865 0.55705124 0.19683278 0.09677403 0.44514183 0.50160837
[37] 0.50160837 0.30572865 0.30572865 0.19683278 0.55705124 0.44514183
[43] 0.37073110 0.19814731 0.20147980 0.27261558 0.36288756 0.21497073
[49] 0.21497073 0.37073110 0.37073110 0.19814731 0.50160837 0.05374460
[55] 0.30572865 0.55705124 0.55705124 0.37073110 0.44514183 0.50160837
[61] 0.50160837 0.30572865 0.36288756 0.21497073 0.55705124 0.44514183
[67] 0.37073110 0.50160837 0.19814731 0.36288756 0.36288756
> dim(model.matrix(fm1)$X)
[1] 72  3
> dim(model.matrix(fm1, "B")$B1)
[1] 71  6
> mf <- update(fm1, method="model.frame")
> str(mf)
'data.frame':	72 obs. of  4 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 ...
 $ (weights): num  0 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "terms")=Classes 'terms', 'formula'  language rating ~ contact + temp
  .. ..- attr(*, "variables")= language list(rating, contact, temp)
  .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:3] "rating" "contact" "temp"
  .. .. .. ..$ : chr [1:2] "contact" "temp"
  .. ..- attr(*, "term.labels")= chr [1:2] "contact" "temp"
  .. ..- attr(*, "order")= int [1:2] 1 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(rating, contact, temp)
  .. ..- attr(*, "dataClasses")= Named chr [1:4] "ordered" "factor" "factor" "numeric"
  .. .. ..- attr(*, "names")= chr [1:4] "rating" "contact" "temp" "(weights)"
> wts <- mf$wts
> dim(model.matrix(fm1)$X[wts > 0, , drop=FALSE])
[1] 0 3
> 
> fm1b <- clm(rating ~ temp, scale=~contact, data=wine, weights=wts)
> summary(fm1b)
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
> pr <- profile(fm1b)
> confint(pr)
                    2.5 %    97.5 %
tempwarm        1.2122036 3.4033767
sca.contactyes -0.6721766 0.4156449
> plot(pr, 1)
> fm1c <- clm(rating ~ temp, nominal=~contact, data=wine, weights=wts)
> summary(fm1c)
formula: rating ~ temp
nominal: ~contact
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H 
 logit flexible  72   -86.21 190.42 6(0)  1.64e-10 4.8e+01

Coefficients:
         Estimate Std. Error z value Pr(>|z|)    
tempwarm    2.519      0.535   4.708  2.5e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                Estimate Std. Error z value
1|2.(Intercept)  -1.3230     0.5623  -2.353
2|3.(Intercept)   1.2464     0.4748   2.625
3|4.(Intercept)   3.5500     0.6560   5.411
4|5.(Intercept)   4.6602     0.8604   5.416
1|2.contactyes   -1.6151     1.1618  -1.390
2|3.contactyes   -1.5116     0.5906  -2.559
3|4.contactyes   -1.6748     0.6488  -2.581
4|5.contactyes   -1.0506     0.8965  -1.172
> pr <- profile(fm1c)
> confint(pr)
            2.5 %   97.5 %
tempwarm 1.516862 3.629098
> plot(pr, 1)
> 
> ## nominal.test(fm1)
> ## scale.test(fm1)
> 
> ## zero out an entire response category:
> wts2 <- 1 * with(wine, rating != "2")
> fm2 <- clm(rating ~ contact + temp, data=wine, weights=wts2)
> fm2
formula: rating ~ contact + temp
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H 
 logit flexible  50   -46.76 103.52 7(0)  3.56e-13 4.2e+01

Coefficients:
contactyes   tempwarm 
     1.593      3.062 

Threshold coefficients:
    1|3     3|4     4|5 
-0.5173  3.6238  5.3657 
> fm2$n ## 72
[1] 72
> fm2$nobs ## 50
[1] 50
> ## Dimension of X and B1, B2 differ:
> dim(model.matrix(fm2)$X)
[1] 72  3
> dim(model.matrix(fm2, "B")$B1)
[1] 50  5
> ## Cannot directly evaluate predictions on the original data:
> try(predict(fm2, newdata=wine), silent=TRUE)[1]
[1] "Error in predict.clm(fm2, newdata = wine) : \n  response factor 'rating' has new levels\n"
> confint(fm2)
               2.5 %   97.5 %
contactyes 0.4216272 2.876550
tempwarm   1.5970977 5.026527
> profile(fm2)
$contactyes
        lroot par.vals.1|3 par.vals.3|4 par.vals.4|5 par.vals.contactyes
1   3.4747852 -1.448538968  2.255543944  3.804439509        -0.447374822
2   3.0276372 -1.308493024  2.407453551  3.967827135        -0.192380808
3   2.5822211 -1.174976592  2.565043880  4.140537415         0.062613207
4   2.1395315 -1.048101582  2.728199207  4.322490592         0.317607222
5   1.7005507 -0.927983229  2.896810974  4.513576014         0.572601236
6   1.2662317 -0.814738724  3.070776301  4.713647671         0.827595251
7   0.8374817 -0.708483160  3.249996208  4.922519451         1.082589265
8   0.4151468 -0.609322477  3.434373593  5.139960291         1.337583280
9   0.0000000 -0.517343412  3.623810880  5.365689290         1.592577294
10 -0.4072702 -0.432600941  3.818207180  5.599370818         1.847571309
11 -0.8060673 -0.355104205  4.017454848  5.840609731         2.102565323
12 -1.1958918 -0.284802472  4.221435388  6.088946911         2.357559338
13 -1.5763436 -0.221573075  4.430015019  6.343855637         2.612553352
14 -1.9471215 -0.165213368  4.643040370  6.604739566         2.867547367
15 -2.3080198 -0.115438462  4.860335061  6.870933284         3.122541381
16 -2.6589228 -0.071885700  5.081697900  7.141706504         3.377535396
17 -2.9997972 -0.034125762  5.306903276  7.416272839         3.632529410
18 -3.3306833 -0.001679093  5.535703885  7.693803770         3.887523425
   par.vals.tempwarm
1        2.840981986
2        2.848907417
3        2.863086150
4        2.883227511
5        2.909025226
6        2.940155181
7        2.976272850
8        3.017010511
9        3.061974244
10       3.110740675
11       3.162853507
12       3.217819929
13       3.275107380
14       3.334141312
15       3.394304942
16       3.454942107
17       3.515364293
18       3.574862708

$tempwarm
        lroot par.vals.1|3 par.vals.3|4 par.vals.4|5 par.vals.contactyes
1   3.4300164   -1.3083735    1.7272573    3.1899360           1.3296606
2   2.8813858   -1.1550431    1.9737448    3.4740817           1.3627908
3   2.3491188   -1.0148529    2.2280703    3.7683888           1.3992293
4   1.8354355   -0.8881411    2.4904988    4.0722937           1.4379639
5   1.3422735   -0.7751174    2.7612331    4.3850175           1.4778669
6   0.8712358   -0.6758159    3.0403898    4.7055883           1.5177348
7   0.4235444   -0.5900448    3.3279674    5.0328829           1.5563515
8   0.0000000   -0.5173434    3.6238109    5.3656893           1.5925773
9  -0.3990485   -0.4569595    3.9275814    5.7027863           1.6254526
10 -0.7737184   -0.4078590    4.2387443    6.0430310           1.6542958
11 -1.1245809   -0.3687765    4.5565843    6.3854348           1.6787668
12 -1.4526173   -0.3383014    4.8802522    6.7292117           1.6988730
13 -1.7591420   -0.3149843    5.2088333    7.0737893           1.7149136
14 -2.0457031   -0.2974396    5.5414220    7.4187857           1.7273842
15 -2.3139760   -0.2844232    5.8771833    7.7639689           1.7368704
16 -2.5656664   -0.2748765    6.2153929    8.1092106           1.7439599
17 -2.8024328   -0.2679376    6.5554528    8.4544482           1.7491850
18 -3.0258332   -0.2629290    6.8968876    8.7996559           1.7529951
19 -3.2372944   -0.2593326    7.2393314    9.1448275           1.7557514
20 -3.4380992   -0.2567600    7.5825083    9.4899654           1.7577335
   par.vals.tempwarm
1          0.6468007
2          0.9918255
3          1.3368503
4          1.6818751
5          2.0268999
6          2.3719247
7          2.7169495
8          3.0619742
9          3.4069990
10         3.7520238
11         4.0970486
12         4.4420734
13         4.7870982
14         5.1321230
15         5.4771478
16         5.8221726
17         6.1671974
18         6.5122222
19         6.8572470
20         7.2022718

attr(,"original.fit")
formula: rating ~ contact + temp
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H 
 logit flexible  50   -46.76 103.52 7(0)  3.56e-13 4.2e+01

Coefficients:
contactyes   tempwarm 
     1.593      3.062 

Threshold coefficients:
    1|3     3|4     4|5 
-0.5173  3.6238  5.3657 
attr(,"class")
[1] "profile.clm"
> plot(slice(fm2), 5)
> step.fm2 <- step(fm2, trace=0)
> fitted(fm2)
 [1] 0.60053892 0.77590129 0.09349820 0.27230685 0.27230685 0.32934749
 [7] 0.32934749 0.37347364 0.10813631 0.77590129 0.60973967 0.32934749
[13] 0.40771198 0.60053892 0.77590129 0.09081575 0.09081575 0.40771198
[19] 0.40771198 0.60053892 0.77590129 0.60973967 0.32934749 0.25729874
[25] 0.60053892 0.09349820 0.77590129 0.60973967 0.60973967 0.25729874
[31] 0.25729874 0.60053892 0.77590129 0.27230685 0.32934749 0.40771198
[37] 0.37347364 0.37347364 0.60973967 0.25729874 0.77590129 0.60973967
[43] 0.60973967 0.25729874 0.40771198 0.37347364 0.77590129 0.60973967
[49] 0.40771198 0.40771198
> ## Scale and nominal effects:
> fm2b <- clm(rating ~ temp, scale=~contact, data=wine, weights=wts2)
> summary(fm2b)
formula: rating ~ temp
scale:   ~contact
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H 
 logit flexible  50   -50.35 110.71 9(0)  2.51e-08 1.2e+02

Coefficients:
         Estimate Std. Error z value Pr(>|z|)  
tempwarm    2.667      1.059   2.519   0.0118 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

log-scale coefficients:
           Estimate Std. Error z value Pr(>|z|)
contactyes -0.09812    0.38836  -0.253    0.801

Threshold coefficients:
    Estimate Std. Error z value
1|3  -1.2009     0.5033  -2.386
3|4   2.3219     1.0427   2.227
4|5   3.7986     1.3876   2.738
> pr <- profile(fm2b)
> confint(pr)
                    2.5 %    97.5 %
tempwarm        1.1319929 6.0851931
sca.contactyes -0.8403892 0.7190847
> plot(pr, 1)
> fm2c <- clm(rating ~ temp, nominal=~contact, data=wine, weights=wts2)
> summary(fm2c)
formula: rating ~ temp
nominal: ~contact
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H 
 logit flexible  50   -46.36 106.72 7(0)  1.22e-12 4.8e+01

Coefficients:
         Estimate Std. Error z value Pr(>|z|)    
tempwarm   3.1100     0.8529   3.646 0.000266 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                Estimate Std. Error z value
1|3.(Intercept)  -0.4250     0.6381  -0.666
3|4.(Intercept)   3.6911     0.9505   3.883
4|5.(Intercept)   4.8968     1.1001   4.451
1|3.contactyes   -1.8993     1.2128  -1.566
3|4.contactyes   -1.7003     0.7432  -2.288
4|5.contactyes   -0.9153     0.9243  -0.990
> pr <- profile(fm2c)
> confint(pr)
            2.5 %   97.5 %
tempwarm 1.623761 5.102185
> plot(pr, 1)
> pred <- predict(fm2c, newdata=wine[!names(wine) %in% "rating"])
> pred <- predict(fm2b, newdata=wine[!names(wine) %in% "rating"])
> 
> ## nominal.test(fm2)
> ## scale.test(fm2)
> 
> ## Different data sets (error):
> try(anova(fm1, fm2), silent=TRUE)[1] ## OK
[1] "Error in anova.clm(fm1, fm2) : \n  models were not all fitted to the same dataset\n"
> 
> ## Test clm.fit:
> wts2 <- 1 * with(wine, rating != "2")
> mf2 <- update(fm2, method="design")
> fm3 <- with(mf2, clm.fit(y, X, weights=wts))
> 
> #################################
> 
> proc.time()
   user  system elapsed 
   2.23    0.23    2.45