R Under development (unstable) (2024-02-26 r85990 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 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(mvord) Loading required package: minqa Loading required package: BB Loading required package: ucminf Loading required package: dfoptim > data(data_mvord_toy) > tolerance <- 1e-6 > > # convert data_mvord_toy into long format > df <- cbind.data.frame("i" = rep(1:100,2), "j" = rep(1:2,each = 100), + "Y" = c(data_mvord_toy$Y1,data_mvord_toy$Y2), + "X1" = rep(data_mvord_toy$X1,2), + "X2" = rep(data_mvord_toy$X2,2), + "f1" = factor(sample(rep(data_mvord_toy$Y2,2)), ordered =F), + "f2" = factor(rep(data_mvord_toy$Y1,2), ordered=F)) > > df_NA <- df[-c(1,90:110),] > > > res <- mvord:::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, + data = df_NA, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "newuoa"), + #se = T, + error.structure = cor_general(~1), + threshold.constraints = c(1,2), + coef.constraints = c(1,1)) > > > > print(res) Call: mvord:::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 2), control = mvord.control(solver = "newuoa")) Thresholds: $`1` 1|2 2|3 -2.243410 1.664118 $`2` 1|2 2|3 -1.856055 1.655492 Coefficients: X1 1 X2 1 1.398855 -2.970139 Parameters of the error structure: [1] 0.8117169 > summary(res) Call: mvord:::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 2), control = mvord.control(solver = "newuoa")) Formula: MMO(Y, i, j) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 99 2 -59.09 133.24 152.78 578 Thresholds: Estimate Std. Error z value Pr(>|z|) 1 1|2 -2.24341 0.49643 -4.5191 6.212e-06 *** 1 2|3 1.66412 0.35086 4.7429 2.107e-06 *** 2 1|2 -1.85605 0.48449 -3.8309 0.0001277 *** 2 2|3 1.65549 0.35398 4.6768 2.913e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Coefficients: Estimate Std. Error z value Pr(>|z|) X1 1 1.39885 0.27219 5.1392 2.759e-07 *** X2 1 -2.97014 0.54436 -5.4562 4.865e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error Structure: Estimate Std. Error z value Pr(>|z|) corr 1 2 0.81172 0.11925 6.8067 9.988e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > print(summary(res)) Call: mvord:::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 2), control = mvord.control(solver = "newuoa")) Formula: MMO(Y, i, j) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 99 2 -59.09 133.24 152.78 578 Thresholds: Estimate Std. Error z value Pr(>|z|) 1 1|2 -2.24341 0.49643 -4.5191 6.212e-06 *** 1 2|3 1.66412 0.35086 4.7429 2.107e-06 *** 2 1|2 -1.85605 0.48449 -3.8309 0.0001277 *** 2 2|3 1.65549 0.35398 4.6768 2.913e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Coefficients: Estimate Std. Error z value Pr(>|z|) X1 1 1.39885 0.27219 5.1392 2.759e-07 *** X2 1 -2.97014 0.54436 -5.4562 4.865e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error Structure: Estimate Std. Error z value Pr(>|z|) corr 1 2 0.81172 0.11925 6.8067 9.988e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 $call mvord:::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 2), control = mvord.control(solver = "newuoa")) $formula MMO(Y, i, j) ~ 0 + X1 + X2 $info c("link", object$rho$link$name) c("threshold", object$rho$threshold) link threshold value mvprobit flexible c("nsubjects", object$rho$n) c("ndim", object$rho$ndim) nsubjects ndim value 99 2 c("logPL", round(-object$rho$objective, 2)) logPL value -59.09 c("CLAIC", ifelse(object$rho$se, round(object$rho$claic, 2), CLAIC value 133.24 c("CLBIC", ifelse(object$rho$se, round(object$rho$clbic, 2), CLBIC value 152.78 c("fevals", if (is.null(object$rho$optRes$fevals)) NA else object$rho$optRes$fevals) fevals value 578 $thresholds Estimate Std. Error z value Pr(>|z|) 1 1|2 -2.243410 0.4964336 -4.519053 6.211673e-06 1 2|3 1.664118 0.3508646 4.742906 2.106744e-06 2 1|2 -1.856055 0.4844930 -3.830922 1.276640e-04 2 2|3 1.655492 0.3539762 4.676845 2.913218e-06 $coefficients Estimate Std. Error z value Pr(>|z|) X1 1 1.398855 0.2721936 5.139190 2.759246e-07 X2 1 -2.970139 0.5443612 -5.456192 4.864539e-08 $error.structure Estimate Std. Error z value Pr(>|z|) corr 1 2 0.8117169 0.1192531 6.806673 9.98816e-12 attr(,"class") [1] "summary.mvord" > > coef(res) X1 1 X2 1 1.398855 -2.970139 > thresholds(res) $`1` 1|2 2|3 -2.243410 1.664118 $`2` 1|2 2|3 -1.856055 1.655492 > AIC(res) [1] 133.2368 > BIC(res) [1] 152.7849 > > > logLik(res) 'log Lik.' -59.08581 (df=7.532609) > nobs(res) [1] 99 > vcov(res) 1.1|2 1.2|3 2.1|2 2.2|3 X1 1 1.1|2 0.246446297 -0.062302524 0.204551432 -0.047162902 -0.0717943408 1.2|3 -0.062302524 0.123105975 -0.043330129 0.092190967 0.0678538108 2.1|2 0.204551432 -0.043330129 0.234733454 -0.028964504 -0.0494469233 2.2|3 -0.047162902 0.092190967 -0.028964504 0.125299118 0.0620299787 X1 1 -0.071794341 0.067853811 -0.049446923 0.062029979 0.0740893439 X2 1 0.228883912 -0.095211934 0.215063901 -0.076358773 -0.1059414930 corr 1 2 0.002456206 -0.004787881 0.008222221 -0.005408435 0.0008225634 X2 1 corr 1 2 1.1|2 0.228883912 0.0024562058 1.2|3 -0.095211934 -0.0047878805 2.1|2 0.215063901 0.0082222212 2.2|3 -0.076358773 -0.0054084352 X1 1 -0.105941493 0.0008225634 X2 1 0.296329064 0.0064684246 corr 1 2 0.006468425 0.0142213048 > terms(res) MMO(Y, i, j) ~ 0 + X1 + X2 attr(,"variables") list(MMO(Y, i, j), X1, X2) attr(,"factors") X1 X2 MMO(Y, i, j) 0 0 X1 1 0 X2 0 1 attr(,"term.labels") [1] "X1" "X2" attr(,"order") [1] 1 1 attr(,"intercept") [1] 0 attr(,"response") [1] 1 attr(,".Environment") > model.matrix(res) [[1]] X1 X2 2 0.92950908 -0.009815203 3 2.80414176 -0.258776324 4 1.44463712 3.901872032 5 -0.19094010 0.049584010 6 1.24822910 2.229198147 7 0.88102184 -0.218488291 8 1.68698898 1.033619596 9 -0.40777779 0.429813068 10 1.57532204 -1.050563512 11 1.63770153 0.773427523 12 0.05187631 2.385234977 13 -0.51553197 1.343198997 14 1.26653819 0.834533920 15 1.29153226 0.520385466 16 2.70892244 -0.507291538 17 1.47670185 0.912693957 18 2.30755180 -0.267032339 19 2.80601269 1.424991878 20 2.03538599 0.809149448 21 1.51651830 2.466845363 22 0.60589648 -0.437537800 23 -0.51024023 0.920473976 24 1.00598984 -0.039202237 25 1.85500699 1.368593878 26 1.72985145 1.650157519 27 1.34875887 0.517696000 28 2.36136747 -0.427844335 29 0.71133787 -0.005731945 30 2.88840366 1.725733179 31 0.98262423 1.710377729 32 1.54108291 -1.198273701 33 0.45995889 0.594696370 34 0.91915060 1.262139340 35 0.97036604 2.867203589 36 0.85474952 0.456297180 37 0.98913740 0.567098045 38 -0.37871876 0.655296921 39 1.67153568 -0.032872886 40 -0.15420545 -1.328973544 41 1.42259779 1.100554626 42 0.61301027 0.585340083 43 -0.24052598 0.103426439 44 1.45696215 -1.117414314 45 -0.41079335 0.477094746 46 2.93390753 0.906281025 47 1.56956065 1.009045870 48 1.93277975 2.086671077 49 -0.54534734 1.546311776 50 1.60161839 -0.547938308 51 0.04923384 -0.290600527 52 2.36036594 0.967298191 53 0.73347641 1.053826861 54 -2.03494577 -0.180978405 55 0.76118009 -2.294959603 56 0.72781328 -0.584586126 57 2.19039680 -0.157097799 58 1.05119609 -0.469290748 59 1.33104803 0.912763202 60 1.53161910 -0.415410286 61 1.61403147 0.051684001 62 1.07671589 -0.298502012 63 0.73912074 -0.517287856 64 0.07812959 1.514988922 65 1.11760616 -0.151507685 66 -1.13841022 0.538215849 67 0.68058040 0.684078667 68 0.25930559 0.951722586 69 0.07420857 0.807880065 70 1.75793379 0.784858233 71 1.35956026 1.164551039 72 0.17076432 0.169025212 73 0.55424788 1.991598133 74 0.32908795 0.672090210 75 0.19908310 0.691160283 76 2.61094670 0.220384992 77 1.01851833 1.275019562 78 -0.08339934 0.164154709 79 2.28944976 0.489470507 80 0.85608726 0.360864997 81 -1.28800884 0.446534053 82 0.91528318 1.829496165 83 1.19905983 1.489769380 84 2.22982583 2.164565787 85 1.00455253 -1.619883648 86 -0.29460764 0.641480817 87 0.60283862 0.322193797 88 2.73948199 0.797883930 89 2.02302095 1.639433348 90 NA NA 91 NA NA 92 NA NA 93 NA NA 94 NA NA 95 NA NA 96 NA NA 97 NA NA 98 NA NA 99 NA NA 100 NA NA attr(,"assign") [1] 1 2 [[2]] X1 X2 2 NA NA 3 NA NA 4 NA NA 5 NA NA 6 NA NA 7 NA NA 8 NA NA 9 NA NA 10 NA NA 11 1.63770153 0.773427523 12 0.05187631 2.385234977 13 -0.51553197 1.343198997 14 1.26653819 0.834533920 15 1.29153226 0.520385466 16 2.70892244 -0.507291538 17 1.47670185 0.912693957 18 2.30755180 -0.267032339 19 2.80601269 1.424991878 20 2.03538599 0.809149448 21 1.51651830 2.466845363 22 0.60589648 -0.437537800 23 -0.51024023 0.920473976 24 1.00598984 -0.039202237 25 1.85500699 1.368593878 26 1.72985145 1.650157519 27 1.34875887 0.517696000 28 2.36136747 -0.427844335 29 0.71133787 -0.005731945 30 2.88840366 1.725733179 31 0.98262423 1.710377729 32 1.54108291 -1.198273701 33 0.45995889 0.594696370 34 0.91915060 1.262139340 35 0.97036604 2.867203589 36 0.85474952 0.456297180 37 0.98913740 0.567098045 38 -0.37871876 0.655296921 39 1.67153568 -0.032872886 40 -0.15420545 -1.328973544 41 1.42259779 1.100554626 42 0.61301027 0.585340083 43 -0.24052598 0.103426439 44 1.45696215 -1.117414314 45 -0.41079335 0.477094746 46 2.93390753 0.906281025 47 1.56956065 1.009045870 48 1.93277975 2.086671077 49 -0.54534734 1.546311776 50 1.60161839 -0.547938308 51 0.04923384 -0.290600527 52 2.36036594 0.967298191 53 0.73347641 1.053826861 54 -2.03494577 -0.180978405 55 0.76118009 -2.294959603 56 0.72781328 -0.584586126 57 2.19039680 -0.157097799 58 1.05119609 -0.469290748 59 1.33104803 0.912763202 60 1.53161910 -0.415410286 61 1.61403147 0.051684001 62 1.07671589 -0.298502012 63 0.73912074 -0.517287856 64 0.07812959 1.514988922 65 1.11760616 -0.151507685 66 -1.13841022 0.538215849 67 0.68058040 0.684078667 68 0.25930559 0.951722586 69 0.07420857 0.807880065 70 1.75793379 0.784858233 71 1.35956026 1.164551039 72 0.17076432 0.169025212 73 0.55424788 1.991598133 74 0.32908795 0.672090210 75 0.19908310 0.691160283 76 2.61094670 0.220384992 77 1.01851833 1.275019562 78 -0.08339934 0.164154709 79 2.28944976 0.489470507 80 0.85608726 0.360864997 81 -1.28800884 0.446534053 82 0.91528318 1.829496165 83 1.19905983 1.489769380 84 2.22982583 2.164565787 85 1.00455253 -1.619883648 86 -0.29460764 0.641480817 87 0.60283862 0.322193797 88 2.73948199 0.797883930 89 2.02302095 1.639433348 90 0.47243275 0.626737533 91 1.84599618 0.014049743 92 1.82154940 -0.503831691 93 -0.33576939 2.284203198 94 2.04311804 -0.367144559 95 -0.12299092 -0.470574785 96 1.69092329 -0.702874783 97 1.36860557 0.132686504 98 1.64851336 0.513773615 99 0.56016969 1.820088710 100 0.13357915 0.989273674 attr(,"assign") [1] 1 2 > > > fitted(res) 2 3 4 5 6 7 8 0.36891923 0.99876532 1.00000000 0.94747098 0.99574990 0.41399042 0.92860303 9 10 11 12 13 14 15 0.65386581 0.99987382 0.89185313 0.99999898 0.99222658 0.84468948 0.86338611 16 17 18 19 20 21 22 0.99975274 0.85582442 0.98530957 0.89311183 0.82948889 0.99825982 0.59916738 23 24 25 26 27 28 29 0.07288144 0.09662462 0.61500280 0.55752716 0.84849279 0.99697454 0.65913835 30 31 32 33 34 35 36 0.12167508 0.91865259 0.99995409 0.73812888 0.23509453 0.99999951 0.89658754 37 38 39 40 41 42 43 0.89343136 0.55538642 0.14612431 0.97003042 0.68586662 0.80654628 0.85607577 44 45 46 47 48 49 50 0.99980434 0.36039939 0.30608862 0.82521396 0.88105298 0.99894407 0.97837027 51 52 53 54 55 56 57 0.15786571 0.83272182 0.36170948 0.48728127 1.00000000 0.05439768 0.95302972 58 59 60 61 62 63 64 0.83889357 0.81427751 0.93529495 0.08808047 0.69299690 0.75404637 0.98186103 65 66 67 68 69 70 71 0.26870870 0.80807378 0.09775575 0.55049358 0.48182468 0.88074557 0.17167737 72 73 74 75 76 77 78 0.89497340 0.99785060 0.16985295 0.49331756 0.87017143 0.26804280 0.86247234 79 80 81 82 83 84 85 0.10058395 0.88098323 0.78999119 0.96803321 0.15516503 0.83945237 0.99999513 86 87 88 89 90 91 92 0.49068966 0.01851734 0.48047059 0.04029408 0.25609806 0.81193757 0.99155380 93 94 95 96 97 98 99 0.99999997 0.98907611 0.33364670 0.99742501 0.55336992 0.80513961 0.99716488 100 0.81470472 > > constraints(res) $X1 X1 1 1|2 1 2|3 1 1|2 1 2|3 1 $X2 X2 1 1|2 1 2|3 1 1|2 1 2|3 1 > > names_constraints(Y ~ 0 + X1 + X2, df_NA) [1] "X1" "X2" > > #predict functions > marginal_predict(res, type = "prob", subjectID = c(2,5,32,88)) 1 2 2 0.3689192 NA 5 0.9474710 NA 32 0.9999745 0.9999754 88 0.5798607 0.5761386 > marginal_predict(res, type = "class", subjectID = c(2,5,32,88)+2) 1 2 34 1 1 NA NA.1 NA.2 > marginal_predict(res, type = "class") 1 2 2 2 3 3 4 1 5 2 6 1 7 3 8 2 9 2 10 3 11 2 2 12 1 1 13 1 1 14 2 2 15 2 2 16 3 3 17 2 2 18 3 3 19 2 2 20 2 2 21 1 1 22 3 3 23 1 1 24 2 2 25 2 2 26 1 1 27 2 2 28 3 3 29 2 2 30 2 2 31 1 1 32 3 3 33 2 2 34 1 1 35 1 1 36 2 2 37 2 2 38 1 1 39 3 3 40 3 3 41 2 2 42 2 2 43 2 2 44 3 3 45 2 1 46 2 2 47 2 2 48 1 1 49 1 1 50 3 3 51 2 2 52 2 2 53 2 1 54 1 1 55 3 3 56 3 3 57 3 3 58 3 3 59 2 2 60 3 3 61 3 3 62 3 3 63 3 3 64 1 1 65 3 3 66 1 1 67 2 2 68 1 1 69 1 1 70 2 2 71 2 2 72 2 2 73 1 1 74 2 2 75 2 2 76 3 3 77 1 1 78 2 2 79 3 3 80 2 2 81 1 1 82 1 1 83 1 1 84 1 1 85 3 3 86 1 1 87 2 2 88 2 2 89 2 1 90 2 91 3 92 3 93 1 94 3 95 2 96 3 97 2 98 2 99 1 100 1 > marginal_predict(res, type = "cum.prob", subjectID = c(2,5,32,88)+8) $`1` 1 2 3 10 1.904192e-14 0.0001261822 1 13 9.931919e-01 0.9999999999 1 40 1.150910e-09 0.0193478930 1 96 NA NA NA $`2` 1 2 3 10 NA NA NA 13 9.978454e-01 1.000000000 1 40 1.151279e-08 0.018945431 1 96 1.403799e-10 0.002574989 1 > predict(res, type = "prob", subjectID = c(3,6,33,55,90)) 3 6 33 55 90 0.9987653 0.9957499 0.7381289 1.0000000 0.2560981 > predict(res, type = "class", subjectID = c(3,6,33,55,90)+1) 1 2 4 1 1 7 3 1 34 1 1 56 3 3 91 1 3 > predict(res, type = "cum.prob", subjectID = c(3,6,33,55,90)+2) 5 8 35 57 92 0.9811677 0.9912079 0.9999995 1.0000000 1.0000000 > joint_probabilities(res, response.cat = c(1,2)) 2 3 4 5 6 7 1.765852e-04 2.036926e-12 1.000000e+00 3.369673e-02 9.957499e-01 1.855512e-05 8 9 10 11 12 13 6.260487e-02 3.459111e-01 1.909584e-14 4.457066e-03 9.678837e-08 9.653155e-04 14 15 16 17 18 19 1.556766e-02 2.413987e-03 2.209344e-14 1.423357e-02 1.625171e-10 8.124942e-03 20 21 22 23 24 25 1.519807e-03 2.123098e-04 3.810048e-06 1.437393e-02 4.846157e-05 3.374894e-02 26 27 28 29 30 31 3.652322e-02 1.940424e-03 4.197753e-12 2.886535e-04 2.466682e-02 9.490683e-03 32 33 34 35 36 37 8.881784e-16 2.561238e-02 3.682428e-02 4.468147e-08 6.124460e-03 8.028036e-03 38 39 40 41 42 43 3.661081e-02 1.023284e-06 9.735815e-10 2.944963e-02 1.960374e-02 1.420228e-02 44 45 46 47 48 49 1.387779e-14 4.022474e-02 7.314049e-05 1.771196e-02 1.340959e-02 1.274785e-04 50 51 52 53 54 55 4.234856e-10 3.723502e-04 1.581193e-03 4.023380e-02 3.889148e-02 0.000000e+00 56 57 58 59 60 61 2.171408e-07 3.205610e-09 1.241225e-07 1.883676e-02 7.785401e-09 4.596234e-06 62 63 64 65 66 67 1.252033e-06 5.386927e-07 2.243356e-03 6.802065e-06 2.026894e-02 2.453040e-02 68 69 70 71 72 73 3.680746e-02 3.903019e-02 3.306715e-03 3.545221e-02 7.485498e-03 2.632552e-04 74 75 76 77 78 79 3.505538e-02 3.871578e-02 6.190749e-08 3.827261e-02 1.338481e-02 2.026696e-05 80 81 82 83 84 85 3.323694e-03 2.182140e-02 3.906614e-03 3.116727e-02 1.743731e-02 0.000000e+00 86 87 88 89 90 91 3.880141e-02 5.582849e-03 6.089830e-05 4.029408e-02 7.417577e-01 1.880569e-01 92 93 94 95 96 97 8.446201e-03 3.368666e-08 1.092389e-02 6.653241e-01 2.574988e-03 5.533699e-01 98 99 100 8.051396e-01 2.835116e-03 1.852900e-01 > joint_probabilities(res, response.cat = res$rho$y) 2 3 4 5 6 7 8 0.36891923 0.99876532 1.00000000 0.94747098 0.99574990 0.41399042 0.92860303 9 10 11 12 13 14 15 0.65386581 0.99987382 0.89185313 0.99999898 0.99222658 0.84468948 0.86338611 16 17 18 19 20 21 22 0.99975274 0.85582442 0.98530957 0.89311183 0.82948889 0.99825982 0.59916738 23 24 25 26 27 28 29 0.07288144 0.09662462 0.61500280 0.55752716 0.84849279 0.99697454 0.65913835 30 31 32 33 34 35 36 0.12167508 0.91865259 0.99995409 0.73812888 0.23509453 0.99999951 0.89658754 37 38 39 40 41 42 43 0.89343136 0.55538642 0.14612431 0.97003042 0.68586662 0.80654628 0.85607577 44 45 46 47 48 49 50 0.99980434 0.36039939 0.30608862 0.82521396 0.88105298 0.99894407 0.97837027 51 52 53 54 55 56 57 0.15786571 0.83272182 0.36170948 0.48728127 1.00000000 0.05439768 0.95302972 58 59 60 61 62 63 64 0.83889357 0.81427751 0.93529495 0.08808047 0.69299690 0.75404637 0.98186103 65 66 67 68 69 70 71 0.26870870 0.80807378 0.09775575 0.55049358 0.48182468 0.88074557 0.17167737 72 73 74 75 76 77 78 0.89497340 0.99785060 0.16985295 0.49331756 0.87017143 0.26804280 0.86247234 79 80 81 82 83 84 85 0.10058395 0.88098323 0.78999119 0.96803321 0.15516503 0.83945237 0.99999513 86 87 88 89 90 91 92 0.49068966 0.01851734 0.48047059 0.04029408 0.25609806 0.81193757 0.99155380 93 94 95 96 97 98 99 0.99999997 0.98907611 0.33364670 0.99742501 0.55336992 0.80513961 0.99716488 100 0.81470472 > joint_probabilities(res, type = "cum.prob", response.cat = c(1,2)) 2 3 4 5 6 7 1.765852e-04 2.036926e-12 1.000000e+00 3.369673e-02 9.957499e-01 1.855512e-05 8 9 10 11 12 13 6.260487e-02 3.459111e-01 1.909584e-14 1.263896e-02 9.999991e-01 9.931919e-01 14 15 16 17 18 19 6.221636e-02 6.131948e-03 2.353673e-14 5.499110e-02 1.870566e-10 2.642223e-02 20 21 22 23 24 25 3.601259e-03 9.984721e-01 5.654021e-06 8.857589e-01 8.258421e-05 2.196485e-01 26 27 28 29 30 31 5.940504e-01 4.764082e-03 4.635847e-12 5.659788e-04 1.233914e-01 9.281433e-01 32 33 34 35 36 37 8.881784e-16 1.312515e-01 5.868936e-01 9.999996e-01 1.858854e-02 2.602565e-02 38 39 40 41 42 43 5.919972e-01 1.439404e-06 1.150909e-09 1.673683e-01 8.653888e-02 5.482613e-02 44 45 46 47 48 49 1.476597e-14 4.006241e-01 1.282209e-04 7.465252e-02 8.944626e-01 9.990715e-01 50 51 52 53 54 55 4.941321e-10 7.481381e-04 3.767689e-03 4.445574e-01 5.261727e-01 0.000000e+00 56 57 58 59 60 61 2.899150e-07 3.869218e-09 1.630182e-07 8.161065e-02 9.561773e-09 6.878737e-06 62 63 64 65 66 67 1.774609e-06 7.406236e-07 9.841044e-01 1.036974e-05 8.283427e-01 1.222861e-01 68 69 70 71 72 73 5.873010e-01 5.208549e-01 8.861244e-03 2.462435e-01 2.383739e-02 9.981139e-01 74 75 76 77 78 79 2.396109e-01 3.195148e-01 7.978050e-08 5.472896e-01 5.058747e-02 3.272855e-05 80 81 82 83 84 85 8.914937e-03 8.118126e-01 9.719398e-01 6.929051e-01 8.568897e-01 0.000000e+00 86 87 88 89 90 91 5.294911e-01 1.659683e-02 1.053952e-04 4.191857e-01 9.978558e-01 1.880624e-01 92 93 94 95 96 97 8.446203e-03 1.000000e+00 1.092389e-02 6.663533e-01 2.574989e-03 5.537371e-01 98 99 100 8.093328e-01 1.000000e+00 9.999948e-01 > > error_structure(res) $`2` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`3` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`4` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`5` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`6` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`7` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`8` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`9` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`10` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`11` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`12` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`13` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`14` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`15` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`16` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`17` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`18` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`19` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`20` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`21` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`22` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`23` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`24` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`25` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`26` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`27` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`28` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`29` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`30` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`31` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`32` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`33` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`34` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`35` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`36` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`37` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`38` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`39` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`40` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`41` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`42` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`43` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`44` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`45` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`46` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`47` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`48` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`49` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`50` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`51` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`52` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`53` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`54` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`55` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`56` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`57` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`58` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`59` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`60` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`61` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`62` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`63` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`64` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`65` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`66` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`67` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`68` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`69` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`70` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`71` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`72` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`73` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`74` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`75` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`76` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`77` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`78` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`79` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`80` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`81` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`82` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`83` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`84` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`85` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`86` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`87` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`88` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`89` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`90` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`91` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`92` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`93` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`94` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`95` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`96` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`97` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`98` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`99` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 $`100` 1 2 1 1.0000000 0.8117169 2 0.8117169 1.0000000 > > newdata <- df[1:20,] > newdata[,"i"] <- rep(1:10,2) > newdata[,"j"] <- rep(1:2, each=10) > #newdata[,"Y"] <- NA > > marginal_predict(res, type = "prob", newdata = newdata)#subjectID = c(2,5,32,88)) 1 2 1 0.9982266 0.9195477 2 0.3689192 0.9999999 3 0.9987653 0.9978454 4 1.0000000 0.8656609 5 0.9474710 0.9012825 6 0.9957499 0.9998640 7 0.4139904 0.8763303 8 0.9286030 0.9909986 9 0.6538658 0.9144494 10 0.9998738 0.8764365 > marginal_predict(res, type = "class", newdata = newdata)#subjectID = c(2,5,32,88)+2) 1 2 1 1 2 2 2 1 3 3 1 4 1 2 5 2 2 6 1 3 7 3 2 8 2 3 9 2 2 10 3 2 > marginal_predict(res, type = "linpred", newdata = newdata) $U 1 2 1 2.9158820 1.661773 2 10000.0000000 5.155857 3 10000.0000000 2.854587 4 7.3248550 2.362471 5 2.0784866 1.394443 6 2.6315274 10000.000000 7 -0.2172443 2.300628 8 2.3742592 10000.000000 9 3.5111442 1.962712 10 10000.0000000 1.211569 $L 1 2 1 -1.000000e+04 -1.849774 2 3.347171e-01 -10000.000000 3 -3.027071e+00 -10000.000000 4 -1.000000e+04 -1.149076 5 -1.829041e+00 -2.117104 6 -1.000000e+04 -3.640623 7 -4.124772e+00 -1.210918 8 -1.533268e+00 -2.365561 9 -3.963834e-01 -1.548835 10 -3.659848e+00 -2.299978 > marginal_predict(res, type = "cum.prob", newdata = newdata) $`1` 1 2 3 1 9.982266e-01 1.0000000000 1 2 1.765852e-04 0.6310807658 1 3 2.036887e-12 0.0012346814 1 4 1.000000e+00 1.0000000000 1 5 3.369673e-02 0.9811677183 1 6 9.957499e-01 1.0000000000 1 7 1.855512e-05 0.4140089727 1 8 6.260487e-02 0.9912079002 1 9 3.459111e-01 0.9997769088 1 10 1.904192e-14 0.0001261822 1 $`2` 1 2 3 1 3.217307e-02 0.9517208072 1 2 9.999999e-01 1.0000000000 1 3 9.978454e-01 0.9999999999 1 4 1.252623e-01 0.9909232098 1 5 1.712553e-02 0.9184080609 1 6 4.270838e-13 0.0001359894 1 7 1.129634e-01 0.9892936772 1 8 2.087485e-09 0.0090013887 1 9 6.071072e-02 0.9751601624 1 10 1.072474e-02 0.8871612751 1 > marginal_predict(res, type = "all.prob", newdata = newdata) $`1` 1 2 3 1 9.982266e-01 1.773423e-03 4.445222e-12 2 1.765852e-04 6.309042e-01 3.689192e-01 3 2.036887e-12 1.234681e-03 9.987653e-01 4 1.000000e+00 1.195710e-13 0.000000e+00 5 3.369673e-02 9.474710e-01 1.883228e-02 6 9.957499e-01 4.250101e-03 3.095435e-11 7 1.855512e-05 4.139904e-01 5.859910e-01 8 6.260487e-02 9.286030e-01 8.792100e-03 9 3.459111e-01 6.538658e-01 2.230912e-04 10 1.904192e-14 1.261822e-04 9.998738e-01 $`2` 1 2 3 1 3.217307e-02 9.195477e-01 4.827919e-02 2 9.999999e-01 1.262367e-07 0.000000e+00 3 9.978454e-01 2.154640e-03 9.692624e-11 4 1.252623e-01 8.656609e-01 9.076790e-03 5 1.712553e-02 9.012825e-01 8.159194e-02 6 4.270838e-13 1.359894e-04 9.998640e-01 7 1.129634e-01 8.763303e-01 1.070632e-02 8 2.087485e-09 9.001387e-03 9.909986e-01 9 6.071072e-02 9.144494e-01 2.483984e-02 10 1.072474e-02 8.764365e-01 1.128387e-01 > > > predict(res, type = "prob", newdata=newdata) 1 2 3 4 5 6 7 8 0.9194344 0.3689191 0.9966107 0.8656609 0.8745944 0.9956139 0.3053300 0.9272670 9 10 0.6274574 0.8764268 > predict(res, type = "class", newdata=newdata, subjectID = c(1:5)) 1 2 1 1 2 2 2 1 3 3 1 4 1 2 5 2 2 > predict(res, type = "cum.prob", newdata=newdata) 1 2 3 4 5 6 7 8 0.9516075 0.9999999 0.9978454 0.9909232 0.9150325 0.9957499 0.4140083 0.9912079 9 10 0.9751520 0.8871613 > > joint_probabilities(res, response.cat = c(1,2), newdata=newdata) 1 2 3 4 5 9.194344e-01 0.000000e+00 -1.292470e-26 8.656609e-01 2.331255e-02 6 7 8 9 10 1.359894e-04 1.037137e-09 7.665328e-03 2.869839e-01 0.000000e+00 > mat <- matrix(sample(1:3, 20, replace=T),ncol=2) > joint_probabilities(res, response.cat = mat, + newdata=newdata) 1 2 3 4 5 6 3.217307e-02 9.902609e-19 1.234681e-03 1.413747e-29 6.741346e-03 3.567632e-19 7 8 9 10 4.285068e-03 1.804112e-16 2.149242e-04 9.742626e-06 > joint_probabilities(res, response.cat = mat, + type = "cum.prob", newdata=newdata) 1 2 3 4 5 6 0.0321730682 0.6310807658 0.0012346814 1.0000000000 0.0171255306 0.0001359894 7 8 9 10 0.1129633905 0.0090013887 1.0000000000 0.0001261822 > joint_probabilities(res, response.cat = mat, + type = "cum.prob", newdata=newdata, subjectID = c(1,3)) 1 3 0.032173068 0.001234681 > > > proc.time() user system elapsed 2.03 0.29 2.31