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. > # test.clm.flex.link.R > > library(ordinal) > > fm <- clm(rating ~ contact + temp, data=wine, link="log-gamma") Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 9.92952e-05 (tol = 1e-06) > fm formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H log-gamma flexible 72 -85.72 185.44 27(179) 9.93e-05 5.9e+01 Coefficients: contactyes tempwarm 0.8614 1.5072 Link coefficient: lambda 0.1615 Threshold coefficients: 1|2 2|3 3|4 4|5 -0.8798 0.6678 1.9807 2.8593 > summary(fm) formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H log-gamma flexible 72 -85.72 185.44 27(179) 9.93e-05 5.9e+01 Coefficients: Estimate Std. Error z value Pr(>|z|) contactyes 0.8614 0.2675 3.220 0.00128 ** tempwarm 1.5072 0.2947 5.114 3.16e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Link coefficients: Estimate Std. Error z value Pr(>|z|) lambda 0.1615 0.5752 0.281 0.779 Threshold coefficients: Estimate Std. Error z value 1|2 -0.8798 0.5003 -1.759 2|3 0.6678 0.3450 1.936 3|4 1.9807 0.3826 5.177 4|5 2.8593 0.4656 6.141 > vcov(fm) 1|2 2|3 3|4 4|5 contactyes 1|2 0.25029666 0.12858124 0.11888119 0.14341729 0.03339111 2|3 0.12858124 0.11899384 0.10688715 0.12363922 0.04117757 3|4 0.11888119 0.10688715 0.14635447 0.15493865 0.05320442 4|5 0.14341729 0.12363922 0.15493865 0.21676904 0.06051103 contactyes 0.03339111 0.04117757 0.05320442 0.06051103 0.07156735 tempwarm 0.00707258 0.03006561 0.05390802 0.06549176 0.01159715 lambda -0.23173314 -0.13979170 -0.12509774 -0.15581154 -0.01157950 tempwarm lambda 1|2 0.00707258 -0.23173314 2|3 0.03006561 -0.13979170 3|4 0.05390802 -0.12509774 4|5 0.06549176 -0.15581154 contactyes 0.01159715 -0.01157950 tempwarm 0.08687099 0.01970691 lambda 0.01970691 0.33082448 > logLik(fm) 'log Lik.' -85.72141 (df=7) > extractAIC(fm) [1] 7.0000 185.4428 > fm2 <- update(fm, link="probit") > anova(fm, fm2) Likelihood ratio tests of cumulative link models: formula: link: threshold: fm2 rating ~ contact + temp probit flexible fm rating ~ contact + temp log-gamma flexible no.par AIC logLik LR.stat df Pr(>Chisq) fm2 6 183.52 -85.761 fm 7 185.44 -85.721 0.0795 1 0.778 > head(model.matrix(fm)$X) (Intercept) contactyes tempwarm 1 1 0 0 2 1 0 0 3 1 1 0 4 1 1 0 5 1 0 1 6 1 0 1 > head(model.frame(fm)) rating contact temp 1 2 no cold 2 3 no cold 3 3 yes cold 4 4 yes cold 5 4 no warm 6 4 no warm > coef(fm) 1|2 2|3 3|4 4|5 contactyes tempwarm lambda -0.8797774 0.6678307 1.9806826 2.8593125 0.8614400 1.5071960 0.1614741 > coef(summary(fm)) Estimate Std. Error z value Pr(>|z|) 1|2 -0.8797774 0.5002966 -1.7585117 7.866049e-02 2|3 0.6678307 0.3449548 1.9359945 5.286837e-02 3|4 1.9806826 0.3825630 5.1774020 2.249970e-07 4|5 2.8593125 0.4655846 6.1413379 8.182929e-10 contactyes 0.8614400 0.2675208 3.2200868 1.281518e-03 tempwarm 1.5071960 0.2947388 5.1136658 3.159661e-07 lambda 0.1614741 0.5751734 0.2807398 7.789100e-01 > nobs(fm) [1] 72 > terms(fm) rating ~ contact + temp attr(,"variables") list(rating, contact, temp) attr(,"factors") contact temp rating 0 0 contact 1 0 temp 0 1 attr(,"term.labels") [1] "contact" "temp" attr(,"order") [1] 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") attr(,"predvars") list(rating, contact, temp) attr(,"dataClasses") rating contact temp "ordered" "factor" "factor" > # profile(fm) # not implemented > confint(fm) Profile intervals not available for models with flexible link function: reporting Wald intervals instead 2.5 % 97.5 % 1|2 -1.860340614 0.1007859 2|3 -0.008268387 1.3439298 3|4 1.230872849 2.7304924 4|5 1.946783410 3.7718416 contactyes 0.337108999 1.3857711 tempwarm 0.929518434 2.0848735 lambda -0.965845153 1.2887933 > > predict(fm, se=TRUE, interval = TRUE) $fit [1] 0.55792544 0.21565445 0.44140779 0.09846706 0.22428321 0.22428321 [7] 0.29090567 0.29090567 0.21042289 0.55792544 0.05388914 0.44140779 [13] 0.20737471 0.48129483 0.29090567 0.33832254 0.55792544 0.21565445 [19] 0.44140779 0.39099052 0.07269587 0.07269587 0.33832254 0.33832254 [25] 0.21565445 0.55792544 0.44140779 0.39099052 0.48129483 0.20737471 [31] 0.29090567 0.31272339 0.55792544 0.21565445 0.09846706 0.44140779 [37] 0.48129483 0.48129483 0.31272339 0.31272339 0.21565445 0.55792544 [43] 0.44140779 0.39099052 0.20737471 0.22428321 0.29090567 0.33832254 [49] 0.21042289 0.21042289 0.39099052 0.39099052 0.20737471 0.48129483 [55] 0.05636253 0.31272339 0.55792544 0.55792544 0.39099052 0.44140779 [61] 0.48129483 0.48129483 0.31272339 0.33832254 0.21042289 0.55792544 [67] 0.44140779 0.39099052 0.48129483 0.20737471 0.33832254 0.33832254 $se.fit [1] 0.09164058 0.09494135 0.08777656 0.05464347 0.08112461 0.08112461 [7] 0.13494788 0.13494788 0.13693811 0.09164058 0.04809909 0.08777656 [13] 0.09512593 0.07373113 0.13494788 0.08236628 0.09164058 0.09494135 [19] 0.08777656 0.09756681 0.06362571 0.06362571 0.08236628 0.08236628 [25] 0.09494135 0.09164058 0.08777656 0.09756681 0.07373113 0.09512593 [31] 0.13494788 0.09451059 0.09164058 0.09494135 0.05464347 0.08777656 [37] 0.07373113 0.07373113 0.09451059 0.09451059 0.09494135 0.09164058 [43] 0.08777656 0.09756681 0.09512593 0.08112461 0.13494788 0.08236628 [49] 0.13693811 0.13693811 0.09756681 0.09756681 0.09512593 0.07373113 [55] 0.03866211 0.09451059 0.09164058 0.09164058 0.09756681 0.08777656 [61] 0.07373113 0.07373113 0.09451059 0.08236628 0.13693811 0.09164058 [67] 0.08777656 0.09756681 0.07373113 0.09512593 0.08236628 0.08236628 $lwr [1] 0.378604418 0.083839590 0.282278456 0.031650566 0.103885180 0.103885180 [7] 0.102183450 0.102183450 0.050316478 0.378604418 0.008885132 0.282278456 [13] 0.077621509 0.342154617 0.102183450 0.199101469 0.378604418 0.083839590 [19] 0.282278456 0.223352982 0.012177667 0.012177667 0.199101469 0.199101469 [25] 0.083839590 0.378604418 0.282278456 0.223352982 0.342154617 0.077621509 [31] 0.102183450 0.161206607 0.378604418 0.083839590 0.031650566 0.282278456 [37] 0.342154617 0.342154617 0.161206607 0.161206607 0.083839590 0.378604418 [43] 0.282278456 0.223352982 0.077621509 0.103885180 0.102183450 0.199101469 [49] 0.050316478 0.050316478 0.223352982 0.223352982 0.077621509 0.342154617 [55] 0.014165414 0.161206607 0.378604418 0.378604418 0.223352982 0.282278456 [61] 0.342154617 0.342154617 0.161206607 0.199101469 0.050316478 0.378604418 [67] 0.282278456 0.223352982 0.342154617 0.077621509 0.199101469 0.199101469 $upr [1] 0.7233159 0.4523801 0.6135565 0.2673887 0.4189772 0.4189772 0.5965756 [8] 0.5965756 0.5727401 0.7233159 0.2657272 0.6135565 0.4485489 0.6233979 [15] 0.5965756 0.5125880 0.7233159 0.4523801 0.6135565 0.5890239 0.3326786 [22] 0.3326786 0.5125880 0.5125880 0.4523801 0.7233159 0.6135565 0.5890239 [29] 0.6233979 0.4485489 0.5965756 0.5186019 0.7233159 0.4523801 0.2673887 [36] 0.6135565 0.6233979 0.6233979 0.5186019 0.5186019 0.4523801 0.7233159 [43] 0.6135565 0.5890239 0.4485489 0.4189772 0.5965756 0.5125880 0.5727401 [50] 0.5727401 0.5890239 0.5890239 0.4485489 0.6233979 0.1988990 0.5186019 [57] 0.7233159 0.7233159 0.5890239 0.6135565 0.6233979 0.6233979 0.5186019 [64] 0.5125880 0.5727401 0.7233159 0.6135565 0.5890239 0.6233979 0.4485489 [71] 0.5125880 0.5125880 > predict(fm, type="class") $fit [1] 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 [39] 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 Levels: 1 2 3 4 5 > newData <- expand.grid(temp = c("cold", "warm"), + contact = c("no", "yes")) > > ## Predicted probabilities in all five response categories for each of > ## the four cases in newData: > predict(fm, newdata=newData, type="prob") $fit 1 2 3 4 5 1 0.210422885 0.55792544 0.2156545 0.01517801 0.0008192112 2 0.014351371 0.20737471 0.4812948 0.22428321 0.0726958750 3 0.053889141 0.39099052 0.4414078 0.09846706 0.0152454937 4 0.001685878 0.05636253 0.3127234 0.33832254 0.2909056666 > predict(fm, newdata=newData, type="class") $fit [1] 2 3 3 4 Levels: 1 2 3 4 5 > > predict(fm, newdata=newData, type="prob", se.fit = TRUE, interval = TRUE) $fit 1 2 3 4 5 1 0.210422885 0.55792544 0.2156545 0.01517801 0.0008192112 2 0.014351371 0.20737471 0.4812948 0.22428321 0.0726958750 3 0.053889141 0.39099052 0.4414078 0.09846706 0.0152454937 4 0.001685878 0.05636253 0.3127234 0.33832254 0.2909056666 $se.fit 1 2 3 4 5 1 0.136938109 0.09164058 0.09494135 0.01696529 0.002047363 2 0.016247015 0.09512593 0.07373113 0.08112461 0.063625711 3 0.048099092 0.09756681 0.08777656 0.05464347 0.019268560 4 0.002046575 0.03866211 0.09451059 0.08236628 0.134947885 $lwr 1 2 3 4 5 1 0.0503164782 0.37860442 0.08383959 0.001663553 6.091054e-06 2 0.0015305241 0.07762151 0.34215462 0.103885180 1.217767e-02 3 0.0088851324 0.22335298 0.28227846 0.031650566 1.249665e-03 4 0.0001557499 0.01416541 0.16120661 0.199101469 1.021835e-01 $upr 1 2 3 4 5 1 0.57274013 0.7233159 0.4523801 0.1247616 0.09939053 2 0.12150066 0.4485489 0.6233979 0.4189772 0.33267865 3 0.26572719 0.5890239 0.6135565 0.2673887 0.16075951 4 0.01797808 0.1988990 0.5186019 0.5125880 0.59657563 > > > ## Aranda-Ordaz link: > fm <- clm(rating ~ contact + temp, data=wine, link="Aranda-Ordaz") Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 1.21375e-05 (tol = 1e-06) > fm formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H Aranda-Ordaz flexible 72 -86.34 186.68 32(214) 1.21e-05 3.9e+02 Coefficients: contactyes tempwarm 1.193 2.076 Link coefficient: lambda 0.4994 Threshold coefficients: 1|2 2|3 3|4 4|5 -1.5284 0.7669 2.6021 3.7962 > summary(fm) formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H Aranda-Ordaz flexible 72 -86.34 186.68 32(214) 1.21e-05 3.9e+02 Coefficients: Estimate Std. Error z value Pr(>|z|) contactyes 1.1927 0.6526 1.828 0.0676 . tempwarm 2.0758 0.8196 2.533 0.0113 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Link coefficients: Estimate Std. Error z value Pr(>|z|) lambda 0.4994 0.7840 0.637 0.524 Threshold coefficients: Estimate Std. Error z value 1|2 -1.5284 0.5775 -2.647 2|3 0.7669 0.8234 0.931 3|4 2.6021 1.4335 1.815 4|5 3.7962 1.9674 1.930 > vcov(fm) 1|2 2|3 3|4 4|5 contactyes tempwarm 1|2 0.3335060 0.2983831 0.4851478 0.6557867 0.2132245 0.2586680 2|3 0.2983831 0.6780281 1.1261332 1.5330566 0.4699073 0.5953101 3|4 0.4851478 1.1261332 2.0550131 2.7755711 0.8214407 1.0831426 4|5 0.6557867 1.5330566 2.7755711 3.8705149 1.1194947 1.4831338 contactyes 0.2132245 0.4699073 0.8214407 1.1194947 0.4259337 0.3965114 tempwarm 0.2586680 0.5953101 1.0831426 1.4831338 0.3965114 0.6717718 lambda 0.2424390 0.5895695 1.0675177 1.4810300 0.4153800 0.5471764 lambda 1|2 0.2424390 2|3 0.5895695 3|4 1.0675177 4|5 1.4810300 contactyes 0.4153800 tempwarm 0.5471764 lambda 0.6146953 > logLik(fm) 'log Lik.' -86.34171 (df=7) > extractAIC(fm) [1] 7.0000 186.6834 > fm2 <- update(fm, link="logit") > anova(fm, fm2) Likelihood ratio tests of cumulative link models: formula: link: threshold: fm2 rating ~ contact + temp logit flexible fm rating ~ contact + temp Aranda-Ordaz flexible no.par AIC logLik LR.stat df Pr(>Chisq) fm2 6 184.98 -86.492 fm 7 186.68 -86.342 0.3004 1 0.5836 > head(model.matrix(fm)$X) (Intercept) contactyes tempwarm 1 1 0 0 2 1 0 0 3 1 1 0 4 1 1 0 5 1 0 1 6 1 0 1 > head(model.frame(fm)) rating contact temp 1 2 no cold 2 3 no cold 3 3 yes cold 4 4 yes cold 5 4 no warm 6 4 no warm > coef(fm) 1|2 2|3 3|4 4|5 contactyes tempwarm lambda -1.5283880 0.7668896 2.6021283 3.7962201 1.1927250 2.0757935 0.4994416 > coef(summary(fm)) Estimate Std. Error z value Pr(>|z|) 1|2 -1.5283880 0.5774998 -2.6465602 0.008131504 2|3 0.7668896 0.8234246 0.9313416 0.351676870 3|4 2.6021283 1.4335317 1.8151871 0.069495154 4|5 3.7962201 1.9673624 1.9295988 0.053656573 contactyes 1.1927250 0.6526360 1.8275502 0.067617093 tempwarm 2.0757935 0.8196168 2.5326389 0.011320750 lambda 0.4994416 0.7840251 0.6370225 0.524110192 > nobs(fm) [1] 72 > terms(fm) rating ~ contact + temp attr(,"variables") list(rating, contact, temp) attr(,"factors") contact temp rating 0 0 contact 1 0 temp 0 1 attr(,"term.labels") [1] "contact" "temp" attr(,"order") [1] 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") attr(,"predvars") list(rating, contact, temp) attr(,"dataClasses") rating contact temp "ordered" "factor" "factor" > # profile(fm) # not implemented > confint(fm) Profile intervals not available for models with flexible link function: reporting Wald intervals instead 2.5 % 97.5 % 1|2 -2.66026679 -0.3965092 2|3 -0.84699296 2.3807722 3|4 -0.20754218 5.4117988 4|5 -0.05973941 7.6521796 contactyes -0.08641803 2.4718680 tempwarm 0.46937403 3.6822130 lambda -1.03721931 2.0361025 > > predict(fm, se=TRUE, interval = TRUE) $fit [1] 0.58209325 0.21518048 0.46054752 0.09103830 0.22383773 0.22383773 [7] 0.29286879 0.29286879 0.18610442 0.58209325 0.06269201 0.46054752 [13] 0.19716621 0.48293293 0.29286879 0.34021076 0.58209325 0.21518048 [19] 0.46054752 0.36913989 0.06939859 0.06939859 0.34021076 0.34021076 [25] 0.21518048 0.58209325 0.46054752 0.36913989 0.48293293 0.19716621 [31] 0.29286879 0.28974268 0.58209325 0.21518048 0.09103830 0.46054752 [37] 0.48293293 0.48293293 0.28974268 0.28974268 0.21518048 0.58209325 [43] 0.46054752 0.36913989 0.19716621 0.22383773 0.29286879 0.34021076 [49] 0.18610442 0.18610442 0.36913989 0.36913989 0.19716621 0.48293293 [55] 0.06897334 0.28974268 0.58209325 0.58209325 0.36913989 0.46054752 [61] 0.48293293 0.48293293 0.28974268 0.34021076 0.18610442 0.58209325 [67] 0.46054752 0.36913989 0.48293293 0.19716621 0.34021076 0.34021076 $se.fit [1] 0.14151187 0.09458707 0.10284385 0.06173008 0.06089158 0.06089158 [7] 0.17401315 0.17401315 0.08458285 0.14151187 0.03342713 0.10284385 [13] 0.06151483 0.16022476 0.17401315 0.14317218 0.14151187 0.09458707 [19] 0.10284385 0.09671233 0.16165491 0.16165491 0.14317218 0.14317218 [25] 0.09458707 0.14151187 0.10284385 0.09671233 0.16022476 0.06151483 [31] 0.17401315 0.08680258 0.14151187 0.09458707 0.06173008 0.10284385 [37] 0.16022476 0.16022476 0.08680258 0.08680258 0.09458707 0.14151187 [43] 0.10284385 0.09671233 0.06151483 0.06089158 0.17401315 0.14317218 [49] 0.08458285 0.08458285 0.09671233 0.09671233 0.06151483 0.16022476 [55] 0.03861752 0.08680258 0.14151187 0.14151187 0.09671233 0.10284385 [61] 0.16022476 0.16022476 0.08680258 0.14317218 0.08458285 0.14151187 [67] 0.10284385 0.09671233 0.16022476 0.06151483 0.14317218 0.14317218 $lwr [1] 0.308147105 0.083805034 0.274985967 0.022684752 0.126708383 0.126708383 [7] 0.073894277 0.073894277 0.071094501 0.308147105 0.021463145 0.274985967 [13] 0.102866209 0.209838512 0.073894277 0.128702471 0.308147105 0.083805034 [19] 0.274985967 0.205891005 0.000551725 0.000551725 0.128702471 0.128702471 [25] 0.083805034 0.308147105 0.274985967 0.205891005 0.209838512 0.102866209 [31] 0.073894277 0.151441083 0.308147105 0.083805034 0.022684752 0.274985967 [37] 0.209838512 0.209838512 0.151441083 0.151441083 0.083805034 0.308147105 [43] 0.274985967 0.205891005 0.102866209 0.126708383 0.073894277 0.128702471 [49] 0.071094501 0.071094501 0.205891005 0.205891005 0.102866209 0.209838512 [55] 0.022286640 0.151441083 0.308147105 0.308147105 0.205891005 0.274985967 [61] 0.209838512 0.209838512 0.151441083 0.128702471 0.071094501 0.308147105 [67] 0.274985967 0.205891005 0.209838512 0.102866209 0.128702471 0.128702471 $upr [1] 0.8132915 0.4511025 0.6577298 0.3017604 0.3643579 0.3643579 0.6825185 [8] 0.6825185 0.4058734 0.8132915 0.1694071 0.6577298 0.3446985 0.7666175 [15] 0.6825185 0.6428533 0.8132915 0.4511025 0.6577298 0.5690698 0.9097002 [22] 0.9097002 0.6428533 0.6428533 0.4511025 0.8132915 0.6577298 0.5690698 [29] 0.7666175 0.3446985 0.6825185 0.4825258 0.8132915 0.4511025 0.3017604 [36] 0.6577298 0.7666175 0.7666175 0.4825258 0.4825258 0.4511025 0.8132915 [43] 0.6577298 0.5690698 0.3446985 0.3643579 0.6825185 0.6428533 0.4058734 [50] 0.4058734 0.5690698 0.5690698 0.3446985 0.7666175 0.1940499 0.4825258 [57] 0.8132915 0.8132915 0.5690698 0.6577298 0.7666175 0.7666175 0.4825258 [64] 0.6428533 0.4058734 0.8132915 0.6577298 0.5690698 0.7666175 0.3446985 [71] 0.6428533 0.6428533 > predict(fm, type="class") $fit [1] 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 [39] 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 Levels: 1 2 3 4 5 > newData <- expand.grid(temp = c("cold", "warm"), + contact = c("no", "yes")) > > ## Predicted probabilities in all five response categories for each of > ## the four cases in newData: > predict(fm, newdata=newData, type="prob") $fit 1 2 3 4 5 1 0.186104423 0.58209325 0.2151805 0.01478353 0.001838317 2 0.026664532 0.19716621 0.4829329 0.22383773 0.069398592 3 0.062692007 0.36913989 0.4605475 0.09103830 0.016582284 4 0.008204432 0.06897334 0.2897427 0.34021076 0.292868788 > predict(fm, newdata=newData, type="class") $fit [1] 2 3 3 4 Levels: 1 2 3 4 5 > > predict(fm, newdata=newData, type="prob", se.fit = TRUE, interval = TRUE) $fit 1 2 3 4 5 1 0.186104423 0.58209325 0.2151805 0.01478353 0.001838317 2 0.026664532 0.19716621 0.4829329 0.22383773 0.069398592 3 0.062692007 0.36913989 0.4605475 0.09103830 0.016582284 4 0.008204432 0.06897334 0.2897427 0.34021076 0.292868788 $se.fit 1 2 3 4 5 1 0.084582846 0.14151187 0.09458707 0.05105461 0.04225889 2 0.018013151 0.06151483 0.16022476 0.06089158 0.16165491 3 0.033427135 0.09671233 0.10284385 0.06173008 0.09203687 4 0.009189042 0.03861752 0.08680258 0.14317218 0.17401315 $lwr 1 2 3 4 5 1 0.071094501 0.30814710 0.08380503 1.557844e-05 4.591076e-23 2 0.006979909 0.10286621 0.20983851 1.267084e-01 5.517250e-04 3 0.021463145 0.20589100 0.27498597 2.268475e-02 2.647331e-07 4 0.000903636 0.02228664 0.15144108 1.287025e-01 7.389428e-02 $upr 1 2 3 4 5 1 0.40587337 0.8132915 0.4511025 0.9352883 1.0000000 2 0.09647029 0.3446985 0.7666175 0.3643579 0.9097002 3 0.16940715 0.5690698 0.6577298 0.3017604 0.9990698 4 0.07033825 0.1940499 0.4825258 0.6428533 0.6825185 > > ######################################################################## > ### Models with scale + flex link (or cauchit link) > ######################################################################## > > fm <- clm(SURENESS ~ PRODID, scale=~PROD, data = soup, link="Aranda-Ordaz") Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 0.000374031 (tol = 1e-06) > summary(fm) formula: SURENESS ~ PRODID scale: ~PROD data: soup link threshold nobs logLik AIC niter max.grad cond.H Aranda-Ordaz flexible 1847 -2674.65 5373.30 57(719) 3.74e-04 9.0e+03 Coefficients: Estimate Std. Error z value Pr(>|z|) PRODID2 1.3604 0.4137 3.289 0.001007 ** PRODID3 2.2066 0.6659 3.314 0.000921 *** PRODID4 1.2553 0.4086 3.072 0.002124 ** PRODID5 1.9820 0.5886 3.367 0.000759 *** PRODID6 2.4715 0.7713 3.204 0.001355 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 log-scale coefficients: Estimate Std. Error z value Pr(>|z|) PRODTest 0.2807 0.1553 1.808 0.0706 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Link coefficients: Estimate Std. Error z value Pr(>|z|) lambda 2.063 1.136 1.816 0.0693 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Threshold coefficients: Estimate Std. Error z value 1|2 -1.3927 0.1491 -9.344 2|3 -0.1449 0.3533 -0.410 3|4 0.3051 0.4699 0.649 4|5 0.6733 0.5802 1.160 5|6 1.7210 0.9594 1.794 > > fm <- clm(SURENESS ~ PRODID, scale=~PROD, data = soup, link="log-gamma") Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 0.00132216 (tol = 1e-06) > summary(fm) formula: SURENESS ~ PRODID scale: ~PROD data: soup link threshold nobs logLik AIC niter max.grad cond.H log-gamma flexible 1847 -2676.84 5377.68 88(1058) 1.32e-03 8.9e+03 Coefficients: Estimate Std. Error z value Pr(>|z|) PRODID2 0.62876 0.09223 6.817 9.27e-12 *** PRODID3 1.02757 0.13449 7.640 2.16e-14 *** PRODID4 0.59157 0.11435 5.173 2.30e-07 *** PRODID5 0.91935 0.12977 7.085 1.39e-12 *** PRODID6 1.13201 0.13927 8.128 4.35e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 log-scale coefficients: Estimate Std. Error z value Pr(>|z|) PRODTest 0.09066 0.18407 0.493 0.622 Link coefficients: Estimate Std. Error z value Pr(>|z|) lambda 0.4749 0.7032 0.675 0.499 Threshold coefficients: Estimate Std. Error z value 1|2 -1.15648 0.45776 -2.526 2|3 -0.47328 0.29577 -1.600 3|4 -0.24744 0.26546 -0.932 4|5 -0.07223 0.24989 -0.289 5|6 0.37914 0.24138 1.571 > > fm <- clm(SURENESS ~ PRODID, scale=~PROD, data = soup, link="cauchit") > summary(fm) formula: SURENESS ~ PRODID scale: ~PROD data: soup link threshold nobs logLik AIC niter max.grad cond.H cauchit flexible 1847 -2679.34 5380.67 11(1) 3.10e-07 3.0e+02 Coefficients: Estimate Std. Error z value Pr(>|z|) PRODID2 0.70550 0.09611 7.341 2.12e-13 *** PRODID3 1.06186 0.13882 7.649 2.02e-14 *** PRODID4 0.66876 0.11911 5.614 1.97e-08 *** PRODID5 1.01135 0.13277 7.617 2.59e-14 *** PRODID6 1.12784 0.14159 7.966 1.64e-15 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 log-scale coefficients: Estimate Std. Error z value Pr(>|z|) PRODTest -0.20801 0.09524 -2.184 0.029 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Threshold coefficients: Estimate Std. Error z value 1|2 -1.70920 0.14977 -11.412 2|3 -0.36737 0.06302 -5.830 3|4 -0.06736 0.05734 -1.175 4|5 0.14657 0.05708 2.568 5|6 0.67063 0.07072 9.483 > > ######################################################################## > ### clm.fit > ######################################################################## > > ## Example with log-gamma: > fm1 <- clm(rating ~ contact + temp, data=wine, link="log-gamma") Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 9.92952e-05 (tol = 1e-06) > summary(fm1) formula: rating ~ contact + temp data: wine link threshold nobs logLik AIC niter max.grad cond.H log-gamma flexible 72 -85.72 185.44 27(179) 9.93e-05 5.9e+01 Coefficients: Estimate Std. Error z value Pr(>|z|) contactyes 0.8614 0.2675 3.220 0.00128 ** tempwarm 1.5072 0.2947 5.114 3.16e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Link coefficients: Estimate Std. Error z value Pr(>|z|) lambda 0.1615 0.5752 0.281 0.779 Threshold coefficients: Estimate Std. Error z value 1|2 -0.8798 0.5003 -1.759 2|3 0.6678 0.3450 1.936 3|4 1.9807 0.3826 5.177 4|5 2.8593 0.4656 6.141 > ## get the model frame containing y and X: > mf1 <- update(fm1, method="design") > names(mf1) [1] "y" "y.levels" "X" "offset" "terms" "contrasts" [7] "xlevels" "weights" "doFit" "control" "link" "threshold" [13] "start" "formulas" > res <- clm.fit(mf1$y, mf1$X, link="log-gamma") ## invoking the factor method Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 9.92952e-05 (tol = 1e-06) > coef(res) 1|2 2|3 3|4 4|5 contactyes tempwarm lambda -0.8797774 0.6678307 1.9806826 2.8593125 0.8614400 1.5071960 0.1614741 > stopifnot(all.equal(coef(res), coef(fm1))) > > ## Example with Aranda-Ordaz: > fm1 <- clm(rating ~ contact + temp, data=wine, link="Aranda-Ordaz") Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 1.21375e-05 (tol = 1e-06) > mf1 <- update(fm1, method="design") > res <- clm.fit(mf1$y, mf1$X, link="Aranda") ## invoking the factor method Changing to 'nlminb' optimizer for flexible link function Warning message: (-1) Model failed to converge with max|grad| = 1.21375e-05 (tol = 1e-06) > stopifnot(all.equal(coef(res), coef(fm1))) > > > proc.time() user system elapsed 8.15 0.37 8.53