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 > library(MASS) > > > #data(data_toy_example) > tolerance <- 1e-6 > > mult.obs <- 2 > sigma <- matrix(c(1,0.8,0.8,1), ncol = 2) > betas <- list(c(0.8,-0.5), + c(0.8,-0.5)) > thresholds <- list(c(-1,1),c(-1,1)) > nobs <- 100 > suppressWarnings(RNGversion("3.5.0")) > set.seed(2017) > errors <- mvrnorm(n = nobs, mu = rep(0, mult.obs), Sigma = sigma) > > X1 <- rnorm(nobs, 0, 1) > X2 <- rnorm(nobs, 0, 1) > > pred <- cbind(X1, X2) > > y <- sapply(1:mult.obs, function(j) pred %*% betas[[j]] + errors[, j], simplify = "array") > y.ord <- sapply(1:mult.obs, function(j) cut(y[, , j], c(min(y[, , j]) - 1, + c(thresholds[[j]]), max(y[, , j]) + 1), + labels = FALSE), simplify = "array") > > predictors.fixed <- lapply(1:mult.obs, function(j) pred) > y <- as.data.frame(y.ord) > > for(i in 1:mult.obs){ + y[, i] <- factor(y[, i], levels = sort(unique(y[, i])), + ordered = TRUE) + } > > > > > data_toy_example <- cbind.data.frame(y, predictors.fixed[[1]]) > colnames(data_toy_example) <- c("Y1","Y2", "X1", "X2") > w <- c(rep(1/20, 20), rep(1/30, 30), rep(1/20, 20), rep(1/30, 30)) > > > > # convert data_toy_example into long format > df <- cbind.data.frame("i" = rep(1:100,2), + "j" = rep(1:2,each = 100), + "Y" = c(data_toy_example$Y1, data_toy_example$Y2), + "X1" = rep(data_toy_example$X1,2), + "X2" = rep(data_toy_example$X2,2), + "f1" = factor(sample(rep(data_toy_example$Y2,2)), + ordered = FALSE), + "f2" = factor(rep(data_toy_example$Y1,2), ordered = FALSE), + w = rep(w,2)) > df$X3 <- cut(df$X2, c(-Inf, -0.2, 0.2, Inf)) > > > > # library(ROI) > # ROI_solver <- function(starting.values, objFun, control){ > # n <- length(starting.values) > # op <- OP(objective = F_objective(objFun, n = n), > # bounds = V_bound(li = seq_len(n), lb = rep.int(-Inf, n))) > # optRes <- ROI_solve(op, solver = "nlminb", > # control = c(list(start = starting.values), > # control)) > # list(optpar = optRes$solution, > # objective = optRes$objval) # a vector of length equal to number of parameters to optimize > # } > # > # > # > res <- mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + link = mvprobit(), + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + coef.constraints = c(1,1), + contrasts = list(f1 = function(x) + contr.treatment(nlevels(df$f1), base = 1), + f2 = "contr.sum"), + control= mvord.control(solver="BFGS",se=TRUE)) Warning message: In mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, link = mvprobit(), : Variables f1 and f2 are absent, the contrasts will be ignored. > options(digits = 22) > > res.summary <- summary(res, short = FALSE) Call: mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 1), contrasts = list(f1 = function(x) contr.treatment(nlevels(df$f1), base = 1), f2 = "contr.sum"), control = mvord.control(solver = "BFGS", se = TRUE)) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -134.91 280.34 294.06 32 Thresholds: Estimate Std. Error z value 1 1|2 -0.96257386663499888702 0.16613952738603396386 -5.7937700000000003087 1 2|3 1.03347036873238828925 0.15004482537615168591 6.8877399999999999736 2 1|2 -0.96257386663499888702 0.16613952738603396386 -5.7937700000000003087 2 2|3 1.03347036873238828925 0.15004482537615168591 6.8877399999999999736 Pr(>|z|) 1 1|2 6.8825e-09 *** 1 2|3 5.6684e-12 *** 2 1|2 6.8825e-09 *** 2 2|3 5.6684e-12 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.63801035062442590373 0.13576747301354935571 4.6992900000000004113 X2 1 -0.42672524816265555714 0.13643665802445156809 -3.1276399999999999757 Pr(>|z|) X1 1 2.6107e-06 *** X2 1 0.0017621 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.854266728221043991631 0.062440889990424582046 13.681210000000000093 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.96257386663519672876, 1.03347036873223707687, -0.96257386663519672876, 1.03347036873223707687), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.63801035062404309883, -0.42672524816263474046), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.85426672822122684536), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.16613952738605441972, 0.15004482537617935822, 0.16613952738605441972, 0.15004482537617935822), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.13576747301357894315, 0.13643665802446261481), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.062440889990360744222), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -134.90867383086322207, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 280.3436634512001433, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 294.05508548271637892, tolerance = tolerance)) > > > res2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + coef.constraints = list(matrix(rep(1,4), ncol = 1), matrix(rep(1,4), ncol = 1))) > > mvord:::check(all.equal(res$beta, res2$beta, tolerance = tolerance)) > mvord:::check(all.equal(res$sebeta, res2$sebeta, tolerance = tolerance)) > ######################################################################## > ## No coefficients > res <- mvord::mvord(formula = MMO(Y) ~ -1, + data = df, + link = mvprobit(), + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + control= mvord.control(solver="BFGS",se=TRUE)) > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ -1, data = df, error.structure = cor_general(~1), link = mvprobit(), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS", se = TRUE)) Formula: MMO(Y) ~ -1 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -153.66 313.51 321.57 26 Thresholds: Estimate Std. Error z value 1 1|2 -0.75479706538162727458 0.13105603927552006427 -5.7593500000000004135 1 2|3 0.86086304364912058507 0.13360642334434574829 6.4432799999999996743 2 1|2 -0.75479706538162727458 0.13105603927552006427 -5.7593500000000004135 2 2|3 0.86086304364912058507 0.13360642334434574829 6.4432799999999996743 Pr(>|z|) 1 1|2 8.4440e-09 *** 1 2|3 1.1692e-10 *** 2 1|2 8.4440e-09 *** 2 2|3 1.1692e-10 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value Pr(>|z|) Error Structure: Estimate Std. Error z value corr 1 2 0.905795171446388658332 0.038855435949046335242 23.311930000000000263 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.75479706538110091785, 0.86086304364935783973, -0.75479706538110091785, 0.86086304364935783973), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(numeric(0)), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.90579517144642240911), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.13105603927553965971, 0.13360642334434202905, 0.13105603927553965971, 0.13360642334434202905), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(numeric(0)), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.038855435949034601573), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -153.66397119528727444, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 313.51350940088383368, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 321.57073678022845797, tolerance = tolerance)) > > #polychor > res <- mvord::mvord(formula = MMO(Y) ~ 1, + data = df, + link = mvprobit(), + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + control= mvord.control(solver="BFGS",se=TRUE)) Note: First threshold for each response is fixed to 0 in order to ensure identifiability! > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 1, data = df, error.structure = cor_general(~1), link = mvprobit(), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS", se = TRUE)) Formula: MMO(Y) ~ 1 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit fix1first 100 2 -153.56 315.46 326.32 54 Thresholds: Estimate Std. Error z value 1 1|2 0.00000000000000000000 0.00000000000000000000 NA 1 2|3 1.61576075463128843168 0.15607113685033857653 10.3527199999999997 2 1|2 0.00000000000000000000 0.00000000000000000000 NA 2 2|3 1.61576075463128843168 0.15607113685033857653 10.3527199999999997 Pr(>|z|) 1 1|2 NA 1 2|3 < 2.22e-16 *** 2 1|2 NA 2 2|3 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error (Intercept) 1 0.73786225415091910573 0.13500906210231616855 (Intercept) 2 0.77255634796640459960 0.14063945841907260492 z value Pr(>|z|) (Intercept) 1 5.4652799999999999159 4.6218e-08 *** (Intercept) 2 5.4931700000000001083 3.9478e-08 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.906373932680408223028 0.038916052671567342991 23.290489999999998361 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(0.0000000000000000000, 1.6157607546978449697, 0.0000000000000000000, 1.6157607546978449697), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.73786225358709867095, 0.77255634837212872057), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.90637393266264265623), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.00000000000000000000, 0.15607113684416776267, 0.00000000000000000000, 0.15607113684416776267), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.13500906208906049422, 0.14063945843173727979), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.038916052675185497378), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -153.5640565887688922, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 315.46144651087109878, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 326.31632228582151356, tolerance = tolerance)) > ####################################################################### > ## cor_general(~factor) > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_general(~f2), + threshold.constraints = c(1,1), + coef.constraints = c(1,1), + contrasts = list(f2 = "contr.sum")) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cor_general(~f2), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 1), contrasts = list(f2 = "contr.sum"), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -134.17 283.39 303 38 Thresholds: Estimate Std. Error z value 1 1|2 -0.90552506790523046298 0.18265647646916075009 -4.9575300000000002143 1 2|3 1.00420547521799474744 0.17467164336592064244 5.7491000000000003212 2 1|2 -0.90552506790523046298 0.18265647646916075009 -4.9575300000000002143 2 2|3 1.00420547521799474744 0.17467164336592064244 5.7491000000000003212 Pr(>|z|) 1 1|2 7.1395e-07 *** 1 2|3 8.9718e-09 *** 2 1|2 7.1395e-07 *** 2 2|3 8.9718e-09 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.64793126338544548837 0.14153585035878471499 4.5778600000000002623 X2 1 -0.42893128334062652351 0.13989943737307897198 -3.0659999999999998366 Pr(>|z|) X1 1 4.6976e-06 *** X2 1 0.0021695 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error corr f21 1 2 0.736678597233872745065 0.178238394308662423393 corr f22 1 2 0.918298263052335417278 0.068858656199963313993 corr f23 1 2 0.837106119713062812870 0.131036082998403219380 z value Pr(>|z|) corr f21 1 2 4.1331100000000002836 3.5789e-05 *** corr f22 1 2 13.3359900000000006770 < 2.22e-16 *** corr f23 1 2 6.3883599999999995944 1.6767e-10 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, + c(-0.90552506790519915469, 1.00420547521745429087, -0.90552506790519915469, 1.00420547521745429087), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, + c(0.64793126338493944871, -0.42893128334048874484), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, + c(0.73667859723415496376, 0.91829826305234418804, 0.83710611971260706632), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, + c(0.18265647646915181279, 0.17467164336594759311, 0.18265647646915181279, 0.17467164336594759311), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, + c(0.14153585035881874332, 0.13989943737309898375), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, + c(0.17823839430849069965, 0.06885865619995423792, 0.13103608299872890330), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -134.17046176709035876, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 283.39468697504094052, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 303.00349482656417877, tolerance = tolerance)) > > ################################################################################## > res <- mvord::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, + data = df, + link = mvlogit(df = 8L), + error.structure = cor_general(~1)) > #threshold.constraints = c(1,1), > #coef.constraints = c(1,1)) > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y, i, j) ~ 0 + X1 + X2, data = df, error.structure = cor_general(~1), link = mvlogit(df = 8L)) Formula: MMO(Y, i, j) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvlogit flexible 100 2 -132.87 285.51 311.28 577 Thresholds: Estimate Std. Error z value 1 1|2 -1.61992737593234226168 0.33738755538628556474 -4.8013799999999999812 1 2|3 1.83247486840983819789 0.36742751710892490591 4.9873099999999999099 2 1|2 -1.61590134864211965038 0.33963204146806280637 -4.7577999999999995850 2 2|3 1.74613474535708190771 0.30112277657338160086 5.7987500000000000711 Pr(>|z|) 1 1|2 1.5757e-06 *** 1 2|3 6.1226e-07 *** 2 1|2 1.9571e-06 *** 2 2|3 6.6812e-09 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 1.28924147226571261804 0.30395056000140890928 4.2416200000000001680 X1 2 0.88420331656816741894 0.29740832947245010587 2.9730300000000000615 X2 1 -0.84290937032158830267 0.30080366307813405369 -2.8021899999999999586 X2 2 -0.71568732513894184333 0.25697673972246015683 -2.7850299999999998946 Pr(>|z|) X1 1 2.2192e-05 *** X1 2 0.0029488 ** X2 1 0.0050757 ** X2 2 0.0053523 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.883893587017596016331 0.071174597846466347573 12.418670000000000542 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > # mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-1.6170817306633420, 1.7855897338762188, -1.6170817306633420, 1.7855897338762188), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$coefficients$Estimate, c(1.07242200717115987, 1.07242200717115987, -0.76715925377701732, -0.76715925377701732), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.85317690560688531), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.28997405649500174, 0.27427389826231802, 0.28997405649500174, 0.27427389826231802), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.24111402270822993, 0.24111402270822993, 0.24156664773225886, 0.24156664773225886), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.063316529381183581), tolerance = tolerance)) > # mvord:::check(all.equal(logLik(res), -135.41665313840898, tolerance = tolerance)) > # mvord:::check(all.equal(AIC(res), 281.35962206629165, tolerance = tolerance)) > # mvord:::check(all.equal(BIC(res), 295.07104409780789, tolerance = tolerance)) > > ################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cov_general(~1), + threshold.constraints = c(1,1), + threshold.values = list(c(-1,NA), + c(-1,NA)), + coef.constraints = c(1,1)) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cov_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 1), threshold.values = list(c(-1, NA), c(-1, NA)), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit fix1first 100 2 -134.64 282.04 298.67 37 Thresholds: Estimate Std. Error z value 1 1|2 -1.0000000000000000000 0.0000000000000000000 NA 1 2|3 1.0826679800833745038 0.2492604341815171376 4.3435199999999998255 2 1|2 -1.0000000000000000000 0.0000000000000000000 NA 2 2|3 1.0826679800833745038 0.2492604341815171376 4.3435199999999998255 Pr(>|z|) 1 1|2 NA 1 2|3 1.4022e-05 *** 2 1|2 NA 2 2|3 1.4022e-05 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.68834946142093178523 0.15128036741853495561 4.5501599999999999824 X2 1 -0.45406517655051331639 0.15704134440910930137 -2.8913700000000002177 Pr(>|z|) X1 1 5.3606e-06 *** X2 1 0.0038356 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.856153637803174327736 0.061822221469819631234 13.8486399999999996169 sigma 1 1.001934583868918871019 0.185187960683921543081 5.4103700000000003456 sigma 2 1.089828199315766354971 0.197593289723851150974 5.5155099999999999127 Pr(>|z|) corr 1 2 < 2.22e-16 *** sigma 1 6.2896e-08 *** sigma 2 3.4777e-08 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-1.000000000000000000000, 1.082667980083374503764, -1.000000000000000000000, 1.082667980083374503764), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.6883494614209317852271, -0.4540651765505133163892), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.8561536378031743277361, 1.0019345838689188710191, 1.0898281993157663549709), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.0000000000000000000000, 0.2492604341815170820862, 0.0000000000000000000000, 0.2492604341815170820862), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1512803674185349001036, 0.1570413444091093568833), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.06182222146981963123435, 0.18518796068392143205905, 0.19759328972385112321852), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -134.6391112878561102661, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 282.0441800225207202857, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 298.6729258905298252103, tolerance = tolerance)) > ######################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_equi(~1), + threshold.constraints = c(1,1), + coef.constraints = cbind(c(1,1),c(1,2))) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cor_equi(~1), link = mvprobit(), coef.constraints = cbind(c(1, 1), c(1, 2)), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -134.84 282.45 299.08 35 Thresholds: Estimate Std. Error z value 1 1|2 -0.96275665235849350321 0.16703865743956156087 -5.7636799999999999145 1 2|3 1.03388029135981107665 0.15203091606548349568 6.8004600000000001714 2 1|2 -0.96275665235849350321 0.16703865743956156087 -5.7636799999999999145 2 2|3 1.03388029135981107665 0.15203091606548349568 6.8004600000000001714 Pr(>|z|) 1 1|2 8.2301e-09 *** 1 2|3 1.0429e-11 *** 2 1|2 8.2301e-09 *** 2 2|3 1.0429e-11 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.63820467410984904610 0.13670492809261164391 4.6684799999999997411 X2 1 -0.44676561119570956659 0.15900879662663833236 -2.8096899999999997988 X2 2 -0.40750155519148029137 0.13850979189698683913 -2.9420399999999999885 Pr(>|z|) X1 1 3.0343e-06 *** X2 1 0.0049589 ** X2 2 0.0032606 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value (Intercept) 1.272900571254562907 0.232216118479298822 5.4815300000000002356 Pr(>|z|) (Intercept) 4.2165e-08 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.9627566523584676350112, 1.0338802913596250032668, -0.9627566523584676350112, 1.0338802913596247812222), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.6382046741095841468905, -0.4467656111955374820255, -0.4075015551913261924177), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$error.structure$Estimate, c(1.272900571254629964457), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.1670386574395647250046, 0.1520309160654948199554, 0.1670386574395647250046, 0.1520309160654948199554), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1367049280926301846328, 0.1590087966266460206555, 0.1385097918969951935608), tolerance = tolerance)) > # mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.2322161184793045674013), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -134.8432321319771745038, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 282.4524217107628487611, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 299.0811675787719536856, tolerance = tolerance)) > > res2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_equi(~1), + threshold.constraints = c(1,1), + coef.constraints = list(X2 = cbind(c(1,1,0,0), c(0,0,1,1)), + X1 = matrix(rep(1,4), ncol = 1))) > > mvord:::check(all.equal(res$beta, res2$beta, tolerance = tolerance)) > mvord:::check(all.equal(res$sebeta, res2$sebeta, tolerance = tolerance)) > > ######################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_ar1(~ 1 + X1), + threshold.constraints = c(1,1), + coef.constraints = c(1,1)) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cor_ar1(~1 + X1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -133.94 280.64 297.27 35 Thresholds: Estimate Std. Error z value 1 1|2 -0.95722951152657487217 0.16558288168119045247 -5.7809699999999999420 1 2|3 1.03746793968333328451 0.15306952616698060732 6.7777599999999997848 2 1|2 -0.95722951152657487217 0.16558288168119045247 -5.7809699999999999420 2 2|3 1.03746793968333328451 0.15306952616698060732 6.7777599999999997848 Pr(>|z|) 1 1|2 7.4272e-09 *** 1 2|3 1.2206e-11 *** 2 1|2 7.4272e-09 *** 2 2|3 1.2206e-11 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.65304208674288677372 0.13584016202004495222 4.8074300000000000921 X2 1 -0.42273136793570653680 0.13820217716269447839 -3.0587900000000001199 Pr(>|z|) X1 1 1.5288e-06 *** X2 1 0.0022223 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error (Intercept) 1.29891554197123726411 0.24517466007551522167 X1 0.29350892362879726916 0.29416988646864372647 z value Pr(>|z|) (Intercept) 5.29792000000000040671 1.1713e-07 *** X1 0.99775000000000002576 0.3184 --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.9572295115261755249492, 1.0374679396833161870717, -0.9572295115261755249492, 1.0374679396833161870717), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.6530420867428219366957, -0.4227313679354475217664), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(1.298915541971320308789, 0.293508923629326790028), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.1655828816811988069002, 0.1530695261669812456962, 0.1655828816811988069002, 0.1530695261669812456962), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1358401620200463122412, 0.1382021771627104100855), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.2451746600755327076815, 0.2941698864686516090572), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -133.938827675498004055, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 280.6436127978045078635, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 297.272358665813612788, tolerance = tolerance)) > > res2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = TRUE, + error.structure = cor_ar1(~1 + X1), + threshold.constraints = c(1,1), + coef.constraints = list(matrix(rep(1,4), ncol = 1), matrix(rep(1,4), + ncol = 1))) > > mvord:::check(all.equal(res$beta, res2$beta, tolerance = tolerance)) > mvord:::check(all.equal(res$sebeta, res2$sebeta, tolerance = tolerance)) > > ######################################################################################## > res <- mvord(formula = MMO2(Y1,Y2) ~ 0 + X1 + X2, + data = data_toy_example, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = TRUE, + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + coef.constraints = cbind(c(1,2),c(NA,1))) > > res.summary <- summary(res, short = FALSE) Call: mvord(formula = MMO2(Y1, Y2) ~ 0 + X1 + X2, data = data_toy_example, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = cbind(c(1, 2), c(NA, 1)), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS")) Formula: MMO2(Y1, Y2) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -137.65 288.06 304.69 63 Thresholds: Estimate Std. Error z value Y1 1|2 -0.89892942527381680673 0.16246077762268909539 -5.5332100000000004059 Y1 2|3 0.98413365940576202906 0.15981657213961206199 6.1578900000000000858 Y2 1|2 -0.89892942527381691775 0.16246077762268909539 -5.5332100000000004059 Y2 2|3 0.98413365940576202906 0.15981657213961206199 6.1578900000000000858 Pr(>|z|) Y1 1|2 3.1442e-08 *** Y1 2|3 7.3718e-10 *** Y2 1|2 3.1442e-08 *** Y2 2|3 7.3718e-10 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.681944859288887594850 0.174961892922583350796 3.89767999999999981142 X1 2 0.468376583108680188960 0.157106006693616268022 2.98127999999999993008 X2 1 -0.052179801987867542590 0.089798215416100737873 -0.58108000000000004093 Pr(>|z|) X1 1 9.712e-05 *** X1 2 0.0028705 ** X2 1 0.5611876 --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr Y1 Y2 0.890390657144618291596 0.063175533024045271713 14.09392000000000067 Pr(>|z|) corr Y1 Y2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.8989294252736504953205, 0.9841336594057152886705, -0.8989294252736503842982, 0.9841336594057153996928), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.68194485928859438494953, 0.46837658310856356003171, -0.05217980198794509860694), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.8903906571446290607597), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.1624607776227008359982, 0.1598165721396140603883, 0.1624607776227008359982, 0.1598165721396140603883), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.17496189292260000414103, 0.15710600669362786985239, 0.08979821541610534529898), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.0631755330240382911855), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -137.6494615446075613363, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 288.064880536023622426, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 304.6936264040327273506, tolerance = tolerance)) > > > res2 <- mvord::mvord(formula = MMO2(Y1,Y2) ~ 0 + X1 + X2, + data = data_toy_example, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = TRUE, + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + coef.constraints = list(X1 = cbind(c(1,1,0,0), c(0,0,1,1)), + X2 = matrix(c(rep(0,2),rep(1,2)), ncol = 1))) > > mvord:::check(all.equal(res$beta, res2$beta, tolerance = tolerance)) > mvord:::check(all.equal(res$sebeta, res2$sebeta, tolerance = tolerance)) > > ######################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + offset(X2), + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = TRUE, + error.structure = cor_general(~1), + threshold.constraints = c(1,2), + coef.constraints = list(X1 = cbind(c(1,1,0,0), c(0,0,1,1)))) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + offset(X2), data = df, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = list(X1 = cbind(c(1, 1, 0, 0), c(0, 0, 1, 1))), threshold.constraints = c(1, 2), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + offset(X2) link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -194.33 403.71 423.31 78 Thresholds: Estimate Std. Error z value 1 1|2 -0.76726526230223290881 0.12023341892342831461 -6.3814599999999996882 1 2|3 0.89190483009970544792 0.13058942975202555981 6.8298399999999999110 2 1|2 -0.79675344626816624327 0.11172347359918627374 -7.1314799999999998192 2 2|3 0.87294080476112090317 0.13662238116423594536 6.3894399999999995643 Pr(>|z|) 1 1|2 1.7540e-10 *** 1 2|3 8.5010e-12 *** 2 1|2 9.9297e-13 *** 2 2|3 1.6649e-10 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.53370081959011461681 0.13636430505431412774 3.9137900000000001022 X1 2 0.34385576067169915104 0.10634167794017658493 3.2334999999999998188 Pr(>|z|) X1 1 9.086e-05 *** X1 2 0.0012228 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.940808527820523377017 0.033697206769472437538 27.919480000000000075 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.7672652623027305107684, 0.8919048301000485068357, -0.7967534462687086982413, 0.8729408047614937160574), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.5337008195901541407480, 0.3438557606717420056519), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.9408085278204662005308), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.1202334189234334493879, 0.1305894297520299174309, 0.1117234735991832761393, 0.1366223811642423013879), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1363643050543318635537, 0.1063416779401736012023), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.03369720676949919391241), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -194.3258468555947047207, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 403.705457152049632441, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 423.314265003572927526, tolerance = tolerance)) > > res2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1, + offset = list(df$X2[1:100], df$X2[101:200]), + data = df, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_general(~1), + threshold.constraints = c(1,2), + coef.constraints = list(X1 = cbind(c(1,1,0,0), c(0,0,1,1)))) > > mvord:::check(all.equal(res$beta, res2$beta, tolerance = tolerance)) > mvord:::check(all.equal(res$sebeta, res2$sebeta, tolerance = tolerance)) > > res3 <- mvord::mvord(formula = MMO2(Y1,Y2) ~ 0 + X1 + offset(X2), + data = data_toy_example, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = TRUE, + error.structure = cor_general(~1), + threshold.constraints = c(1,2), + coef.constraints = list(X1 = cbind(c(1,1,0,0), c(0,0,1,1)))) > mvord:::check(all.equal(res$beta, res3$beta, tolerance = tolerance)) > mvord:::check(all.equal(res$sebeta, res3$sebeta, tolerance = tolerance)) > > ######################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + #index = c("i", "j"), + link = mvlogit(), + #solver = "newuoa", + #se = TRUE, + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + coef.constraints = cbind(c(NA, NA), c(1,2)), + coef.values = cbind(c(1, 1), c(NA,NA))) > #rho <- res$rho > > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cor_general(~1), link = mvlogit(), coef.constraints = cbind(c(NA, NA), c(1, 2)), coef.values = cbind(c(1, 1), c(NA, NA)), threshold.constraints = c(1, 1)) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvlogit flexible 100 2 -135.52 281.57 295.28 150 Thresholds: Estimate Std. Error z value 1 1|2 -1.58328892697782519505 0.25851085233511844619 -6.1246499999999999275 1 2|3 1.75864250788365961142 0.24251594148065969136 7.2516600000000002169 2 1|2 -1.58328892697782519505 0.25851085233511844619 -6.1246499999999999275 2 2|3 1.75864250788365961142 0.24251594148065969136 7.2516600000000002169 Pr(>|z|) 1 1|2 9.0882e-10 *** 1 2|3 4.1170e-13 *** 2 1|2 9.0882e-10 *** 2 2|3 4.1170e-13 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X2 1 -0.77306097199511181195 0.26732809082944253287 -2.8918099999999999916 X2 2 -0.72430349472230626251 0.24403984874581674536 -2.9679700000000002191 Pr(>|z|) X2 1 0.0038303 ** X2 2 0.0029977 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.855894467973149053464 0.060599187425615069769 14.123860000000000525 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > tolerance2 <- 1e-4 > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-1.583288915548024533564, 1.758642501421955106622, -1.583288915548024533564, 1.758642501421955106622), tolerance = tolerance2)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(-0.7730608423417825170176, -0.7243032675736262859800), tolerance = tolerance2)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.8558944465934688050623), tolerance = tolerance2)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.2585108482987530376107, 0.2425159370850978601819, 0.2585108482987530376107, 0.2425159370850978601819), tolerance = tolerance2)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.2673280928737744588375, 0.2440398419906193439033), tolerance = tolerance2)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.06059919526643173237623), tolerance = tolerance2)) > mvord:::check(all.equal(logLik(res)[[1]], -135.5210996495467554723, tolerance = tolerance2)) > mvord:::check(all.equal(AIC(res), 281.568515088567210114, tolerance = tolerance2)) > mvord:::check(all.equal(BIC(res), 295.2799371200834457341, tolerance = tolerance2)) > > # res1 <- mvord::mvord(formula = MMO(Y) ~ 1 + X1 + X2, > # data = df, > # index = c("i", "j"), > # link = mvprobit(), > # solver = "BFGS", > # se = TRUE, > # error.structure = cor_equi(~1), > # threshold.constraints = c(1,1), > # coef.constraints = list(Intercept = matrix(rep(1,6), ncol = 1), > # X2 = cbind(c(1,1,1,0,0,0), c(0,0,0,1,1,1)), X1 = matrix(rep(1,6), ncol = 1))) > > > ######################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X3, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = T, + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + coef.constraints = c(1,1)) > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X3, data = df, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X3 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -136.05 284.87 301.5 34 Thresholds: Estimate Std. Error z value 1 1|2 -1.21270698725037129329 0.21131407873814866960 -5.7388799999999999812 1 2|3 0.75845153901176631805 0.18860790313554207098 4.0213099999999997181 2 1|2 -1.21270698725037129329 0.21131407873814866960 -5.7388799999999999812 2 2|3 0.75845153901176631805 0.18860790313554207098 4.0213099999999997181 Pr(>|z|) 1 1|2 9.5302e-09 *** 1 2|3 5.7874e-05 *** 2 1|2 9.5302e-09 *** 2 2|3 5.7874e-05 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error X1 1 0.582346304092164257504 0.133874274950955496344 X3(-0.2,0.2] 1 0.065949564660296139396 0.394250784463485104414 X3(0.2, Inf] 1 -0.639867558058156649103 0.239596299509497845381 z value Pr(>|z|) X1 1 4.34994999999999976126 1.3617e-05 *** X3(-0.2,0.2] 1 0.16728000000000001202 0.8671512 X3(0.2, Inf] 1 -2.67060999999999992838 0.0075714 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.857975765661059597278 0.062391829886098891522 13.751409999999999911 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-1.212706987250618206886, 0.758451539011639419563, -1.212706987250618206886, 0.758451539011639419563), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.58234630409261545214150, 0.06594956465951176682871, -0.63986755805831896370961), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.8579757656611409766256), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.2113140787381464491546, 0.1886079031355630819533, 0.2113140787381464491546, 0.1886079031355630819533), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1338742749509412854891, 0.3942507844634823288565, 0.2395962995094961522913), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.06239182988606906121731), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -136.0526219854373835005, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 284.8712014176832667545, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 301.499947285692371679, tolerance = tolerance)) > ######################################################################################## > res <- mvord::mvord(formula = MMO(Y) ~ 1 + X1 * X2, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + #se = T, + error.structure = cor_general(~1), + threshold.constraints = c(1,1), + threshold.values = list(c(-1,NA), + c(-1,NA)), + coef.constraints = c(1,1)) > res.summary <- summary(res, short = TRUE) Call: mvord::mvord(formula = MMO(Y) ~ 1 + X1 * X2, data = df, error.structure = cor_general(~1), link = mvprobit(), coef.constraints = c(1, 1), threshold.constraints = c(1, 1), threshold.values = list(c(-1, NA), c(-1, NA)), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 1 + X1 * X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit fix1first 100 2 -134.63 282.02 298.65 34 Thresholds: Estimate Std. Error z value 1 1|2 -1.00000000000000000000 0.00000000000000000000 NA 1 2|3 1.00269513578732416548 0.22781659168150072969 4.4013299999999997425 Pr(>|z|) 1 1|2 NA 1 2|3 1.0759e-05 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error (Intercept) 1 -0.021913582738025384755 0.167100624889530507966 X1 1 0.621621046245967145971 0.138540865251019801319 X2 1 -0.427624354284268537452 0.140667605271754703189 X1:X2 1 -0.102567255999721287929 0.150381573286269054623 z value Pr(>|z|) (Intercept) 1 -0.13114000000000000656 0.8956645 X1 1 4.48690999999999995396 7.2262e-06 *** X2 1 -3.03996000000000021757 0.0023661 ** X1:X2 1 -0.68205000000000004512 0.4952094 --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.853574805224538435411 0.063090989221389948138 13.529270000000000351 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-1.000000000000000000000, 1.002695135787324165477), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(-0.02191358273802538475516, 0.62162104624596714597118, -0.42762435428426853745165, -0.10256725599972128792903), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.8535748052245384354109), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.0000000000000000000000, 0.2278165916815008962271), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1671006248895306467439, 0.1385408652510198845853, 0.1406676052717547309445, 0.1503815732862690546234), tolerance = tolerance)) > #mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.063097735777156036), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.06309098922138993426056), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -134.6255603757719541136, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 282.0170781983524079806, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 298.6458240663615129051, tolerance = tolerance)) > ######################################################################################## > df_NA <- df[-c(1,90:110),] > > > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df_NA, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_general(~1), + threshold.constraints = c(1,2)) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), threshold.constraints = c(1, 2), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 99 2 -119.45 258.71 284.4 51 Thresholds: Estimate Std. Error z value 1 1|2 -1.02012637737351052891 0.19996641229615236268 -5.1014900000000000801 1 2|3 1.14180364925250188435 0.22379090693932096601 5.1021000000000000796 2 1|2 -0.90662822300298373523 0.19059651595704599214 -4.7567899999999996297 2 2|3 0.99945117269137251359 0.17403069740199664239 5.7429600000000000648 Pr(>|z|) 1 1|2 3.3699e-07 *** 1 2|3 3.3590e-07 *** 2 1|2 1.9669e-06 *** 2 2|3 9.3036e-09 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.83733616646579622778 0.18898470050446236623 4.4307100000000003703 X1 2 0.49821075302565920317 0.16630384260856054257 2.9957899999999999530 X2 1 -0.44740279449415937263 0.18763575165140966572 -2.3844199999999999839 X2 2 -0.35390389963325241673 0.14489624593607869363 -2.4424600000000000755 Pr(>|z|) X1 1 9.3924e-06 *** X1 2 0.0027374 ** X2 1 0.0171060 * X2 2 0.0145874 * --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.910650307821200821756 0.070553566078427734642 12.907220000000000582 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, + c(-1.0201263773742383911269, 1.1418036492524186176212, -0.9066282230038859024646, 0.9994511726914282467860), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.8373361664658567349306, 0.4982107530258053085248, -0.4474027944949834356692, -0.3539038996339812781500), + tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.9106503078209118307029), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.1999664122961394840949, 0.2237909069393826111405, 0.1905965159570501277209, 0.1740306974019930341679), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.1889847005044749950198, 0.1663038426085674814647, 0.1876357516513920686840, 0.1448962459360401411335), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.07055356607861483497768), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -119.4526280916837492896, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 258.7052561833675099479, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 284.3969426996999345647, tolerance = tolerance)) > > #weights > df_NA$weights <- 0.5 > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df_NA, + link = mvprobit(), + weights.name = "weights", + error.structure = cor_general(~1), + threshold.constraints = c(1,2)) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), threshold.constraints = c(1, 2), weights.name = "weights") Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 99 2 -59.73 139.25 164.94 472 Thresholds: Estimate Std. Error z value 1 1|2 -1.02012825054024824922 0.39993315330717621459 -2.5507499999999998508 1 2|3 1.14180477624979115348 0.44758333888939638712 2.5510399999999999743 2 1|2 -0.90661616506435449558 0.38119351435790660432 -2.3783599999999998076 2 2|3 0.99945807115822449251 0.34806068836082060258 2.8715099999999997848 Pr(>|z|) 1 1|2 0.0107492 * 1 2|3 0.0107401 * 2 1|2 0.0173897 * 2 2|3 0.0040852 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.83733313484785520942 0.37797048927211568392 2.2153399999999998649 X1 2 0.49821352561752080268 0.33260835741097843909 1.4979000000000000092 X2 1 -0.44740651027549038776 0.37527221427704604562 -1.1922200000000000575 X2 2 -0.35390613973431295225 0.28979246443025280522 -1.2212400000000001032 Pr(>|z|) X1 1 0.026737 * X1 2 0.134160 X2 1 0.233176 X2 2 0.221995 --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.91064688734819310145 0.14111078549223804446 6.4534200000000003783 Pr(>|z|) corr 1 2 1.0936e-10 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, + c(-1.0201283793713198377873, 1.1418048707357395521456, -0.9066161533747665313143, 0.9994581570190742558779), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.8373332969208180376341, 0.4982136050008418859392, -0.4474066175154422508875, -0.3539061566757252808024), + tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.9106469196938988819312), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.3999331442401471981007, 0.4475833313547766811880, 0.3811935020996289336104, 0.3480606771313752845209), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.3779704499331774103510, 0.3326083385001772918521, 0.3752721890043577146479, 0.2897924558296670061175), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.1411107414052604758226), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -59.72631403934757088336, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 139.2526280786951531354, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 164.944314595027606174, tolerance = tolerance)) > ########################################################################################## > df_23 <- df > > ##For the first response (first 100 obs), replace the level 3 by 2. > ## Now the first response has 2 levels and the second one 3 levels > df_23$Y[which(df_23$Y[1:100] == 3)] <- 2 > > ## Must be converted to integer. Ordered factor is not OK, as we have diff number of responses. > df_23$Y <- as.integer(df_23$Y) > > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df_23, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_general(~1), + threshold.constraints = c(1,2)) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df_23, error.structure = cor_general(~1), link = mvprobit(), threshold.constraints = c(1, 2), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -109.15 235.69 258.34 53 Thresholds: Estimate Std. Error z value 1 1|2 -0.97243369041066884950 0.20101574846708089583 -4.8376000000000001222 2 1|2 -0.97734151029870541816 0.18730956260704784144 -5.2177899999999999281 2 2|3 1.01224642882145121625 0.15713428754520597508 6.4419199999999996464 Pr(>|z|) 1 1|2 1.3142e-06 *** 2 1|2 1.8107e-07 *** 2 2|3 1.1797e-10 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.82969552084045039564 0.22409294296782589218 3.7024599999999998623 X1 2 0.53332636576263858785 0.14638351173008076755 3.6433499999999998664 X2 1 -0.39479319852743410824 0.18931396166562203254 -2.0853899999999998549 X2 2 -0.38625976505626685720 0.13976399874813491553 -2.7636599999999997834 Pr(>|z|) X1 1 0.00021352 *** X1 2 0.00026911 *** X2 1 0.03703403 * X2 2 0.00571576 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value corr 1 2 0.85295804978378753081 0.10226978918641670135 8.3402700000000002944 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > > mvord:::check(all.equal(res.summary$thresholds$Estimate, + c(-0.9724336904108278334391, -0.9773415102989740921302, 1.0122464288213952610107), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.8296955208405573101160, 0.5333263657626993170524, -0.3947931985275920929723, -0.3862597650564197349077), + tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.8529580497839671648919), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.2010157484670624383760, 0.1873095626070425123721, 0.1571342875452053644558), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.2240929429677874784588, 0.1463835117300714416810, 0.1893139616655838686210, 0.1397639987481304746364), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.1022697891863128538681), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -109.147957798495554016, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 235.6872199448172011671, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 258.3408737360180111864, tolerance = tolerance)) > > ########################################################################################## > ## make three responses > suppressWarnings(RNGversion("3.5.0")) > set.seed(100) > Y3 <- sample(1:3, nobs, replace = TRUE) > dat_3 <- cbind(Y3, data_toy_example) > res <- mvord::mvord(formula = MMO2(Y1, Y2, Y3) ~ 0, + data = dat_3, + link = mvprobit(), + error.structure = cor_general(~1), + control= mvord.control(solver="BFGS",se=TRUE)) Warning message: Responses are unordered. Natural ordering is used. > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO2(Y1, Y2, Y3) ~ 0, data = dat_3, error.structure = cor_general(~1), link = mvprobit(), control = mvord.control(solver = "BFGS", se = TRUE)) Formula: MMO2(Y1, Y2, Y3) ~ 0 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 3 -565.19 1163.2 1205.96 52 Thresholds: Estimate Std. Error z value Y1 1|2 -0.70654808839442972968 0.14413615638980573075 -4.9019500000000002515 Y1 2|3 0.84166389562761123599 0.15094934798831763367 5.5758000000000000895 Y2 1|2 -0.80564017944257926285 0.14777311388303038253 -5.4518700000000004380 Y2 2|3 0.87834342403448395498 0.15259796098861780345 5.7559300000000002129 Y3 1|2 -0.52550677773958742733 0.13805473124275705055 -3.8065099999999998381 Y3 2|3 0.43872841951737284738 0.13625729487946994234 3.2198500000000001009 Pr(>|z|) Y1 1|2 9.4890e-07 *** Y1 2|3 2.4639e-08 *** Y2 1|2 4.9842e-08 *** Y2 2|3 8.6165e-09 *** Y3 1|2 0.00014094 *** Y3 2|3 0.00128256 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value Pr(>|z|) Error Structure: Estimate Std. Error corr Y1 Y2 0.908534877536116503016 0.039421737990092485526 corr Y1 Y3 -0.068886066523438385656 0.133978933777760100821 corr Y2 Y3 -0.150586382344864100347 0.129319546675897428800 z value Pr(>|z|) corr Y1 Y2 23.0465499999999998693 < 2e-16 *** corr Y1 Y3 -0.5141599999999999504 0.60714 corr Y2 Y3 -1.1644499999999999851 0.24424 --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > > mvord:::check(all.equal(res.summary$thresholds$Estimate, + c(-0.70654808839454852354, 0.84166389562778798350 ,-0.80564017944211274713, + 0.87834342403465937021, -0.52550677773899212575, 0.43872841951666374793), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$Estimate, + c(0.908534877536098850470 , -0.068886066520331828977, -0.150586382347683067628), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, + c(0.14413615638981053246, 0.15094934798830694778, 0.14777311388303066009, + 0.15259796098859162994, 0.13805473124275577379, 0.13625729487947765839), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c(0.039421737990101658744,0.133978933777893494117,0.129319546675785962409), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -565.1857975044365502981, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 1163.199625517649110407, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 1205.960928690734590418, tolerance = tolerance)) > ########################################## > ### fixed coefs fixed thresholds > ########################################### > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df_NA, + link = mvprobit(), + error.structure = cor_general(~1), + threshold.values = list(c(-1,1), c(-1,1)), + coef.values = matrix(c(0.8,-0.5, 0.8,-0.5),ncol=2, byrow=TRUE), + control= mvord.control(solver="BFGS",se=TRUE)) We suggest to include an intercept in the model (formula = y ~ 1 + ...)> > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df_NA, error.structure = cor_general(~1), link = mvprobit(), coef.values = matrix(c(0.8, -0.5, 0.8, -0.5), ncol = 2, byrow = TRUE), threshold.values = list(c(-1, 1), c(-1, 1)), control = mvord.control(solver = "BFGS", se = TRUE)) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit fixall 99 2 -124.63 251.28 253.9 20 Thresholds: Estimate Std. Error z value Pr(>|z|) 1 1|2 -1 0 NA NA 1 2|3 1 0 NA NA 2 1|2 -1 0 NA NA 2 2|3 1 0 NA NA Coefficients: Estimate Std. Error z value Pr(>|z|) Error Structure: Estimate Std. Error z value corr 1 2 0.854364908267473466275 0.064212481459284348473 13.305279999999999774 Pr(>|z|) corr 1 2 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 > > > mvord:::check(all.equal(res.summary$error.structure$Estimate, c(0.854364908267473466275), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$error.structure$`Std. Error`, c( 0.064212481459284348473), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -124.6298356200302066554, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 251.280079403325714793, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 253.9016800682576047166, tolerance = tolerance)) > > > ########################################## > ### fixed cor_general > ########################################### > #polychor > res <- mvord::mvord(formula = MMO(Y) ~ 1, + data = df, + link = mvprobit(), + error.structure = cor_general(~1, value = rep(0.5, 100), fixed = TRUE), + threshold.constraints = c(1,1), + control= mvord.control(solver="BFGS",se=TRUE)) Note: First threshold for each response is fixed to 0 in order to ensure identifiability! > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 1, data = df, error.structure = cor_general(~1, value = rep(0.5, 100), fixed = TRUE), link = mvprobit(), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS", se = TRUE)) Formula: MMO(Y) ~ 1 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit fix1first 100 2 -167.98 342.15 350.2 38 Thresholds: Estimate Std. Error z value 1 1|2 0.00000000000000000000 0.00000000000000000000 NA 1 2|3 1.75287953185495593011 0.12271683314006323617 14.283939999999999415 2 1|2 0.00000000000000000000 0.00000000000000000000 NA 2 2|3 1.75287953185495593011 0.12271683314006323617 14.283939999999999415 Pr(>|z|) 1 1|2 NA 1 2|3 < 2.22e-16 *** 2 1|2 NA 2 2|3 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error (Intercept) 1 0.80138487301048844103 0.13627526573459616821 (Intercept) 2 0.83916182795533666994 0.14543826409239238306 z value Pr(>|z|) (Intercept) 1 5.8806300000000000239 4.0870e-09 *** (Intercept) 2 5.7698799999999996757 7.9326e-09 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value Pr(>|z|) > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(0.0000000000000000000, 1.75287953185495593011, 0.0000000000000000000, 1.75287953185495593011), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.80138487301048844103, 0.83916182795533666994), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.00000000000000000000, 0.12271683314006323617, 0.00000000000000000000, 0.12271683314006323617), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.13627526573459616821, 0.14543826409239238306), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -167.9806162341531887705 , tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 342.1467994786156623377, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 350.2040268579602866339, tolerance = tolerance)) > > ########################################## > ### fixed cor_equi > ########################################### > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + #index = c("i", "j"), + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_equi(~0+X1, fixed = TRUE), + threshold.constraints = c(1,1), + coef.constraints = cbind(c(1,1),c(1,2))) > > res.summary <- summary(res, short = FALSE) Call: mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, data = df, error.structure = cor_equi(~0 + X1, fixed = TRUE), link = mvprobit(), coef.constraints = cbind(c(1, 1), c(1, 2)), threshold.constraints = c(1, 1), control = mvord.control(solver = "BFGS")) Formula: MMO(Y) ~ 0 + X1 + X2 link threshold nsubjects ndim logPL CLAIC CLBIC fevals mvprobit flexible 100 2 -161.63 333.79 347.5 32 Thresholds: Estimate Std. Error z value 1 1|2 -0.971882533677594362409 0.106998427476495336363 -9.0831499999999998352 1 2|3 1.049862299350324157388 0.099201905823190800193 10.5830900000000003303 2 1|2 -0.971882533677594362409 0.106998427476495336363 -9.0831499999999998352 2 2|3 1.049862299350324157388 0.099201905823190800193 10.5830900000000003303 Pr(>|z|) 1 1|2 < 2.22e-16 *** 1 2|3 < 2.22e-16 *** 2 1|2 < 2.22e-16 *** 2 2|3 < 2.22e-16 *** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Coefficients: Estimate Std. Error z value X1 1 0.641833309819806796526 0.085721117636041219146 7.4874599999999995603 X2 1 -0.450733949174911885382 0.188620706054736408941 -2.3896299999999999208 X2 2 -0.406790878681185075205 0.157255593575891583491 -2.5868099999999998317 Pr(>|z|) X1 1 7.022e-14 *** X2 1 0.0168653 * X2 2 0.0096868 ** --- Signif. codes: 0 '***' 0.001000000000000000020817 '**' 0.01000000000000000020817 '*' 0.05000000000000000277556 '.' 0.1000000000000000055511 ' ' 1 Error Structure: Estimate Std. Error z value Pr(>|z|) > > options(digits = 22) > > # paste(format(res.summary$thresholds$Estimate), collapse = ",") > # paste(format(res.summary$coefficients$Estimate), collapse = ",") > # paste(format(res.summary$error.structure$Estimate), collapse = ",") > mvord:::check(all.equal(res.summary$thresholds$Estimate, c(-0.971882533677633775326, 1.049862299350415195676, -0.971882533677633775326, 1.049862299350415195676), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$Estimate, c(0.641833309819990649459, -0.450733949174842496443, -0.406790878681064282940), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$thresholds$`Std. Error`, c(0.106998427476494920030 ,0.099201905823191174894, 0.106998427476494920030, 0.099201905823191174894), tolerance = tolerance)) > mvord:::check(all.equal(res.summary$coefficients$`Std. Error`, c(0.085721117636038332566,0.188620706054737241608, 0.157255593575882451907), tolerance = tolerance)) > mvord:::check(all.equal(logLik(res)[[1]], -161.6316436136202128182, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 333.7896030167141248057, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 347.5010250482303604258, tolerance = tolerance)) > ########################################## > ### try all fixed > ########################################### > res <- mvord::mvord(formula = MMO(Y) ~ 0, + data = df, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_equi(~1, value = 0.3, fixed = TRUE), + threshold.values = list(c(-1,1), c(-1,1))) We suggest to include an intercept in the model (formula = y ~ 1 + ...)> > > mvord:::check(all.equal(logLik(res)[[1]], -179.9334335736217553858, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 359.8668671472435107717, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 359.8668671472435107717, tolerance = tolerance)) > > res <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2, + data = df, + link = mvprobit(), + control = mvord.control(solver = "BFGS"), + error.structure = cor_equi(~1, value = 0.8, fixed = TRUE), + threshold.values = list(c(-1,1), c(-1,1)), + coef.values = rbind(c(0.8,-0.5), c(0.8,-0.5))) We suggest to include an intercept in the model (formula = y ~ 1 + ...)> > mvord:::check(all.equal(logLik(res)[[1]], -136.4789503795024359079, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 272.9579007590048718157, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 272.9579007590048718157, tolerance = tolerance)) > > ## make three responses > suppressWarnings(RNGversion("3.5.0")) > set.seed(100) > Y3 <- factor(sample(1:3, nobs, replace = TRUE), ordered = TRUE) > dat_3 <- cbind(Y3, data_toy_example) > Rfixed <- matrix(c(1,0.9,0.2, 0.9, 1, 0.5, 0.2,0.5,1), nrow = 3) > res <- mvord::mvord(formula = MMO2(Y1, Y2, Y3) ~ 0, + data = dat_3, + link = mvprobit(), + error.structure = cor_general(~1, value = Rfixed[lower.tri(Rfixed)], fixed = TRUE), + threshold.values = list(c(-1,1), c(-1,1), c(-1,1)), + control= mvord.control(solver="BFGS",se=TRUE)) We suggest to include an intercept in the model (formula = y ~ 1 + ...)> > mvord:::check(all.equal(logLik(res)[[1]], -647.8423017485074524302, tolerance = tolerance)) > mvord:::check(all.equal(AIC(res), 1295.68460349701490486, tolerance = tolerance)) > mvord:::check(all.equal(BIC(res), 1295.68460349701490486, tolerance = tolerance)) > > ############################### > # tests for predict functions # > ############################### > data_toy_example$X3 <- cut(data_toy_example$X2, c(-Inf, -0.2, 0.2, Inf)) > > dat_train <- data_toy_example[1:90, ] > dat_test <- data_toy_example[-c(1:90), ] > > ## mvord2 > ### example with factor covariates > res_train <- mvord::mvord(formula = MMO2(Y1, Y2) ~ 0 + X1 + X2 + X3, + data = dat_train, + link = mvprobit(), + contrasts = list(X3 = "contr.sum")) > tolerance <- 0.00001 > #### predict in sample > phat_in_mvord2 <- predict(res_train) > mvord:::check(all.equal( + unname(phat_in_mvord2[1:10]), c(0.1464235557313708635530, 0.4321397714843485116099, 0.5968270247741616074677, + 0.0166994091681152423412, 0.5854390261665686212567, 0.1371096523508007203329, + 0.1675153858158427988556, 0.5777761575951960715258, 0.5680673832081621910106, + 0.2403876149104018367098), tolerance = tolerance)) > > pjhat_in_mvord2 <- joint_probabilities(res_train, response.cat = c(1,2)) > pjhat_in_mvord2_1 <- joint_probabilities(res_train, response.cat = c(1,1)) > > p1 <- joint_probabilities(res_train, response.cat = c(1,NA)) > p2 <- joint_probabilities(res_train, response.cat = c(2,NA)) > min(rowSums(cbind(pjhat_in_mvord2/p1, pjhat_in_mvord2_1/p1))) [1] 0.9999166767820104606557 > > mvord:::check(all.equal( + unname(pjhat_in_mvord2[1:10]), c(0.0414494114970503091388565, 0.0007763727128762831775077, 0.0258791356876044154056160, + 0.0166994147132800141442033, 0.0468724891285079159342075, 0.0021647187725574168482012, + 0.0172638911663273442176347, 0.0094837831334130262561644, 0.0190099916431371585012755, + 0.0051418125225030086866695), tolerance = tolerance)) > > pmhat_in_mvord2 <- marginal_predict(res_train) > mvord:::check(all.equal( + unname(c(pmhat_in_mvord2[1:5,1:2])), c(0.16369555938725077748330, 0.51206442964897258551815, 0.70085752018010938346748, + 0.07780827569973108870371, 0.68584056337028265204481, 0.23658538526249328626250, + 0.62104312126876426436439, 0.65391318126703112945108, 0.62425836768740627924501, + 0.64830276740170678095865), tolerance = tolerance)) > > > > #### predict with newdata > phat_new_mvord2 <- predict(res_train, newdata = dat_test) > mvord:::check(all.equal( + unname(phat_new_mvord2), c(0.50788676686680833682885, 0.53701228897073405299523, 0.41931860775657509021741, + 0.56837451159716412263379, 0.46158948077333511461617, 0.02998639798438534898040, + 0.06314648572131303927435, 0.21981855723244725364651, 0.57800150278644024659513, + 0.57873631730860186639376), tolerance = tolerance)) > > pjhat_new_mvord2 <- joint_probabilities(res_train, response.cat = c(1,2), newdata = dat_test) > mvord:::check(all.equal( + unname(pjhat_new_mvord2), c(0.098893552409077684073324, 0.007494474951204033175145, 0.223630828339767884216371, + 0.012572935869544998865877, 0.125624085603978530301106, 0.021774380504568757732642, + 0.063146525118494334360975, 0.000117581041524150720079, 0.017153034070206046868279, + 0.049461049315582095164956), tolerance = tolerance)) > > pmhat_new_mvord2 <- marginal_predict(res_train, newdata = dat_test) > mvord:::check(all.equal( + unname(c(pmhat_new_mvord2[1:5, 1:2])), c(0.5879858728748059704117, 0.6245339414959365509361, 0.6429603932672183219665, + 0.6645168757663386660539, 0.5872139451846044577721, 0.6110788197600012239263, + 0.6185525371757517598681, 0.4254630210115252220149, 0.6355172734548747426331, + 0.4875151835646630571475), tolerance = tolerance)) > > > ### example with offsets > res_train2 <- mvord::mvord(formula = MMO2(Y1, Y2) ~ 0 + X1 + offset(X2) + X3, + data = dat_train, + link = mvprobit(), + contrasts = list(X3 = "contr.sum")) > #### predict in sample > phat_in_offset_mvord2 <- predict(res_train2) > mvord:::check(all.equal( + unname(phat_in_offset_mvord2[1:10]), c(0.04651187641465470007374, 0.37252604618539558734014, 0.52685355820976886853657, + 0.01765062793667723783919, 0.56828750653284465510495, 0.13097924136563171559899, + 0.01252594308629018104995, 0.55895559865178667813268, 0.57808854222406069744977, + 0.47885651300164144839044), tolerance = tolerance)) > > pjhat_in_mvord2 <- joint_probabilities(res_train, response.cat = c(1,2)) > mvord:::check(all.equal( + unname(pjhat_in_mvord2[1:10]), c(0.0414494114970503091388565, 0.0007763727128762831775077, 0.0258791356876044154056160, + 0.0166994147132800141442033, 0.0468724891285079159342075, 0.0021647187725574168482012, + 0.0172638911663273442176347, 0.0094837831334130262561644, 0.0190099916431371585012755, + 0.0051418125225030086866695), tolerance = tolerance)) > > pmhat_in_mvord2 <- marginal_predict(res_train) > mvord:::check(all.equal( + unname(c(pmhat_in_mvord2[1:5,1:2])), c(0.16369555938725077748330, 0.51206442964897258551815, 0.70085752018010938346748, + 0.07780827569973108870371, 0.68584056337028265204481, 0.23658538526249328626250, + 0.62104312126876426436439, 0.65391318126703112945108, 0.62425836768740627924501, + 0.64830276740170678095865), tolerance = tolerance)) > > > > #### predict with newdata > phat_new_offset_mvord2 <- predict(res_train2, newdata = dat_test) > mvord:::check(all.equal( + unname(phat_new_offset_mvord2), c(0.53788568792038271570988, 0.47226651862936941395077, 0.59736906923858912321634, + 0.51730764069679358030385, 0.18750767549430336078586, 0.01993645212535863353587, + 0.05609361865051418205574, 0.15695387935814331115125, 0.57079407463923570453801, + 0.46447453401393973271283), tolerance = tolerance)) > > > pjhat_new_offset_mvord2 <- joint_probabilities(res_train2, response.cat = c(1,2), + newdata = dat_test) > mvord:::check(all.equal( + unname(pjhat_new_offset_mvord2), c(0.0839907843603347747940546, 0.0053737551172480801930931, 0.1657125298469289131908511, + 0.0095039724989203211436006, 0.1110141336202859763115924, 0.0308028085793187883512090, + 0.0560936429791563384572584, 0.0000659201973896017534571, 0.0161492922521473680763648, + 0.0556615607647627450016437), tolerance = tolerance)) > > pmhat_new_offset_mvord2 <- marginal_predict(res_train2, newdata = dat_test) > mvord:::check(all.equal( + unname(c(pmhat_new_offset_mvord2[1:5, 1:2])), c(0.6410297683731078777214, 0.5530240404054447278526, 0.7630836817892157064591, + 0.6103504462453783752096, 0.2985228989322880055468, 0.6269027239883980806567, + 0.5461645309760024824541, 0.6047542519330054711091, 0.5786243958547179211394, + 0.2011002687374978670221), tolerance = tolerance)) > > ## mvord > dat_train <- df[df$i %in% 1:90, ] > dat_test <- df[!(df$i %in% 1:90), ] > > res_train <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + X2 + X3, + data = dat_train, + link = mvprobit(), + contrasts = list(X3 = "contr.sum")) > ### example with factor covariates > #### predictions in sample > phat_in_mvord <- predict(res_train) > mvord:::check(all.equal( + unname(phat_in_mvord[1:10]), c(0.1464235557313708635530, 0.4321397714843485116099, 0.5968270247741616074677, + 0.0166994091681152423412, 0.5854390261665686212567, 0.1371096523508007203329, + 0.1675153858158427988556, 0.5777761575951960715258, 0.5680673832081621910106, + 0.2403876149104018367098), tolerance = tolerance)) > > pjhat_in_mvord <- joint_probabilities(res_train, response.cat = c(1,2)) > mvord:::check(all.equal( + unname(pjhat_in_mvord[1:10]), c(0.0414494114970503091388565, 0.0007763727128762831775077, 0.0258791356876044154056160, + 0.0166994147132800141442033, 0.0468724891285079159342075, 0.0021647187725574168482012, + 0.0172638911663273442176347, 0.0094837831334130262561644, 0.0190099916431371585012755, + 0.0051418125225030086866695), tolerance = tolerance)) > > pmhat_in_mvord <- marginal_predict(res_train) > mvord:::check(all.equal( + unname(c(pmhat_in_mvord[1:5,1:2])), c(0.16369555938725077748330, 0.51206442964897258551815, 0.70085752018010938346748, + 0.07780827569973108870371, 0.68584056337028265204481, 0.23658538526249328626250, + 0.62104312126876426436439, 0.65391318126703112945108, 0.62425836768740627924501, + 0.64830276740170678095865), tolerance = tolerance)) > > > #### predict with newdata > phat_new_mvord <- predict(res_train, newdata = dat_test) > mvord:::check(all.equal( + unname(phat_new_mvord), c(0.50788676686680833682885, 0.53701228897073405299523, 0.41931860775657509021741, + 0.56837451159716412263379, 0.46158948077333511461617, 0.02998639798438534898040, + 0.06314648572131303927435, 0.21981855723244725364651, 0.57800150278644024659513, + 0.57873631730860186639376), tolerance = tolerance)) > > pjhat_new_mvord <- joint_probabilities(res_train, response.cat = c(1,2), newdata = dat_test) > mvord:::check(all.equal( + unname(pjhat_new_mvord), c(0.098893552409077684073324, 0.007494474951204033175145, 0.223630828339767884216371, + 0.012572935869544998865877, 0.125624085603978530301106, 0.021774380504568757732642, + 0.063146525118494334360975, 0.000117581041524150720079, 0.017153034070206046868279, + 0.049461049315582095164956), tolerance = tolerance)) > > pmhat_new_mvord <- marginal_predict(res_train, newdata = dat_test) > mvord:::check(all.equal( + unname(c(pmhat_new_mvord[1:5,1:2])), + c(0.5879858728748059704117, 0.6245339414959365509361, 0.6429603932672183219665, + 0.6645168757663386660539, 0.5872139451846044577721, 0.6110788197600012239263, + 0.6185525371757517598681, 0.4254630210115252220149, 0.6355172734548747426331, + 0.4875151835646630571475), tolerance = tolerance)) > > > > > > > ### example with fixed error structure > res_train2_vec <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cor_ar1(~1, value = rep(0,length(unique(dat_train$i))), fixed = TRUE)) > > mvord:::check(all.equal(logLik(res_train2_vec)[[1]], -148.0213745871528772113, tolerance = tolerance)) > > > res_train2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cor_ar1(~1, value = 0, fixed = TRUE)) > mvord:::check(all.equal(logLik(res_train2_vec)[[1]], -148.0213745871528772113, tolerance = tolerance)) > > phat_mvord <- predict(res_train2) > > phat_mp_mvord <- marginal_predict(res_train2) > > pjhat_mvord <- joint_probabilities(res_train2, response.cat = c(1,2)) > > # res_train2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, > # data = dat_train, > # link = mvprobit(), > # cor_ar1(~1)) > # > # res_train2_free$rho$optpar > # res_train2_free$rho$varGamma > > res_train2_vec <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cor_equi(~1, value = rep(0,length(unique(dat_train$i))), fixed = TRUE)) > mvord:::check(all.equal(logLik(res_train2_vec)[[1]], -148.0213745871528772113, tolerance = tolerance)) > > res_train2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cor_equi(~1, value = 0, fixed = TRUE)) > mvord:::check(all.equal(logLik(res_train2)[[1]], -148.0213745871528772113, tolerance = tolerance)) > > phat_mvord <- predict(res_train2) > > phat_mp_mvord <- marginal_predict(res_train2) > > pjhat_mvord <- joint_probabilities(res_train2, response.cat = c(1,2)) > > > # cor_general > res_train2_vec <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cor_general(~1, value = matrix(rep(0,length(unique(dat_train$i))), ncol = 1), fixed = TRUE)) > mvord:::check(all.equal(logLik(res_train2_vec)[[1]], -148.0213745871528772113, tolerance = tolerance)) > > res_train2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cor_general(~1, value = 0, fixed = TRUE)) > mvord:::check(all.equal(logLik(res_train2)[[1]], -148.0213745871528772113, tolerance = tolerance)) > > phat_mvord <- predict(res_train2) > > phat_mp_mvord <- marginal_predict(res_train2) > > pjhat_mvord <- joint_probabilities(res_train2, response.cat = c(1,2)) > > #cov_general > res_train2_vec <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + cov_general(~1, value = matrix(rep(c(1,0,1),length(unique(dat_train$i))), ncol = 3, byrow= TRUE), fixed = TRUE)) Note: First threshold for each response is fixed to 0 in order to ensure identifiability! > mvord:::check(all.equal(logLik(res_train2_vec)[[1]], -187.3253828363251329847, tolerance = tolerance)) > > res_train2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 +X2, + data = dat_train, + link = mvprobit(), + error.structure = cov_general(~1, value = c(1,0,1), fixed = TRUE)) Note: First threshold for each response is fixed to 0 in order to ensure identifiability! > mvord:::check(all.equal(logLik(res_train2)[[1]], -187.3253828363251329847, tolerance = tolerance)) > > error_structure(res_train2, type ="sigmas") $`1` [,1] [,2] [1,] 1 0 [2,] 0 1 $`2` [,1] [,2] [1,] 1 0 [2,] 0 1 $`3` [,1] [,2] [1,] 1 0 [2,] 0 1 $`4` [,1] [,2] [1,] 1 0 [2,] 0 1 $`5` [,1] [,2] [1,] 1 0 [2,] 0 1 $`6` [,1] [,2] [1,] 1 0 [2,] 0 1 $`7` [,1] [,2] [1,] 1 0 [2,] 0 1 $`8` [,1] [,2] [1,] 1 0 [2,] 0 1 $`9` [,1] [,2] [1,] 1 0 [2,] 0 1 $`10` [,1] [,2] [1,] 1 0 [2,] 0 1 $`11` [,1] [,2] [1,] 1 0 [2,] 0 1 $`12` [,1] [,2] [1,] 1 0 [2,] 0 1 $`13` [,1] [,2] [1,] 1 0 [2,] 0 1 $`14` [,1] [,2] [1,] 1 0 [2,] 0 1 $`15` [,1] [,2] [1,] 1 0 [2,] 0 1 $`16` [,1] [,2] [1,] 1 0 [2,] 0 1 $`17` [,1] [,2] [1,] 1 0 [2,] 0 1 $`18` [,1] [,2] [1,] 1 0 [2,] 0 1 $`19` [,1] [,2] [1,] 1 0 [2,] 0 1 $`20` [,1] [,2] [1,] 1 0 [2,] 0 1 $`21` [,1] [,2] [1,] 1 0 [2,] 0 1 $`22` [,1] [,2] [1,] 1 0 [2,] 0 1 $`23` [,1] [,2] [1,] 1 0 [2,] 0 1 $`24` [,1] [,2] [1,] 1 0 [2,] 0 1 $`25` [,1] [,2] [1,] 1 0 [2,] 0 1 $`26` [,1] [,2] [1,] 1 0 [2,] 0 1 $`27` [,1] [,2] [1,] 1 0 [2,] 0 1 $`28` [,1] [,2] [1,] 1 0 [2,] 0 1 $`29` [,1] [,2] [1,] 1 0 [2,] 0 1 $`30` [,1] [,2] [1,] 1 0 [2,] 0 1 $`31` [,1] [,2] [1,] 1 0 [2,] 0 1 $`32` [,1] [,2] [1,] 1 0 [2,] 0 1 $`33` [,1] [,2] [1,] 1 0 [2,] 0 1 $`34` [,1] [,2] [1,] 1 0 [2,] 0 1 $`35` [,1] [,2] [1,] 1 0 [2,] 0 1 $`36` [,1] [,2] [1,] 1 0 [2,] 0 1 $`37` [,1] [,2] [1,] 1 0 [2,] 0 1 $`38` [,1] [,2] [1,] 1 0 [2,] 0 1 $`39` [,1] [,2] [1,] 1 0 [2,] 0 1 $`40` [,1] [,2] [1,] 1 0 [2,] 0 1 $`41` [,1] [,2] [1,] 1 0 [2,] 0 1 $`42` [,1] [,2] [1,] 1 0 [2,] 0 1 $`43` [,1] [,2] [1,] 1 0 [2,] 0 1 $`44` [,1] [,2] [1,] 1 0 [2,] 0 1 $`45` [,1] [,2] [1,] 1 0 [2,] 0 1 $`46` [,1] [,2] [1,] 1 0 [2,] 0 1 $`47` [,1] [,2] [1,] 1 0 [2,] 0 1 $`48` [,1] [,2] [1,] 1 0 [2,] 0 1 $`49` [,1] [,2] [1,] 1 0 [2,] 0 1 $`50` [,1] [,2] [1,] 1 0 [2,] 0 1 $`51` [,1] [,2] [1,] 1 0 [2,] 0 1 $`52` [,1] [,2] [1,] 1 0 [2,] 0 1 $`53` [,1] [,2] [1,] 1 0 [2,] 0 1 $`54` [,1] [,2] [1,] 1 0 [2,] 0 1 $`55` [,1] [,2] [1,] 1 0 [2,] 0 1 $`56` [,1] [,2] [1,] 1 0 [2,] 0 1 $`57` [,1] [,2] [1,] 1 0 [2,] 0 1 $`58` [,1] [,2] [1,] 1 0 [2,] 0 1 $`59` [,1] [,2] [1,] 1 0 [2,] 0 1 $`60` [,1] [,2] [1,] 1 0 [2,] 0 1 $`61` [,1] [,2] [1,] 1 0 [2,] 0 1 $`62` [,1] [,2] [1,] 1 0 [2,] 0 1 $`63` [,1] [,2] [1,] 1 0 [2,] 0 1 $`64` [,1] [,2] [1,] 1 0 [2,] 0 1 $`65` [,1] [,2] [1,] 1 0 [2,] 0 1 $`66` [,1] [,2] [1,] 1 0 [2,] 0 1 $`67` [,1] [,2] [1,] 1 0 [2,] 0 1 $`68` [,1] [,2] [1,] 1 0 [2,] 0 1 $`69` [,1] [,2] [1,] 1 0 [2,] 0 1 $`70` [,1] [,2] [1,] 1 0 [2,] 0 1 $`71` [,1] [,2] [1,] 1 0 [2,] 0 1 $`72` [,1] [,2] [1,] 1 0 [2,] 0 1 $`73` [,1] [,2] [1,] 1 0 [2,] 0 1 $`74` [,1] [,2] [1,] 1 0 [2,] 0 1 $`75` [,1] [,2] [1,] 1 0 [2,] 0 1 $`76` [,1] [,2] [1,] 1 0 [2,] 0 1 $`77` [,1] [,2] [1,] 1 0 [2,] 0 1 $`78` [,1] [,2] [1,] 1 0 [2,] 0 1 $`79` [,1] [,2] [1,] 1 0 [2,] 0 1 $`80` [,1] [,2] [1,] 1 0 [2,] 0 1 $`81` [,1] [,2] [1,] 1 0 [2,] 0 1 $`82` [,1] [,2] [1,] 1 0 [2,] 0 1 $`83` [,1] [,2] [1,] 1 0 [2,] 0 1 $`84` [,1] [,2] [1,] 1 0 [2,] 0 1 $`85` [,1] [,2] [1,] 1 0 [2,] 0 1 $`86` [,1] [,2] [1,] 1 0 [2,] 0 1 $`87` [,1] [,2] [1,] 1 0 [2,] 0 1 $`88` [,1] [,2] [1,] 1 0 [2,] 0 1 $`89` [,1] [,2] [1,] 1 0 [2,] 0 1 $`90` [,1] [,2] [1,] 1 0 [2,] 0 1 > > > phat_mvord <- predict(res_train2) > #warnings() > > phat_mp_mvord <- marginal_predict(res_train2) > > pjhat_mvord <- joint_probabilities(res_train2, response.cat = c(1,2)) > > ### example with offsets > res_train2 <- mvord::mvord(formula = MMO(Y) ~ 0 + X1 + offset(X2) + X3, + data = dat_train, + link = mvprobit(), + contrasts = list(X3 = "contr.sum")) > #### predict in sample > phat_in_offset_mvord <- predict(res_train2) > mvord:::check(all.equal( + unname(phat_in_offset_mvord[1:10]), c(0.04651187641465470007374, 0.37252604618539558734014, 0.52685355820976886853657, + 0.01765062793667723783919, 0.56828750653284465510495, 0.13097924136563171559899, + 0.01252594308629018104995, 0.55895559865178667813268, 0.57808854222406069744977, + 0.47885651300164144839044), tolerance = tolerance)) > > pjhat_in_offset_mvord <- joint_probabilities(res_train2, response.cat = c(1,2)) > mvord:::check(all.equal( + unname(pjhat_in_offset_mvord[1:10]), c(0.0472419001589241827065990, 0.0005045621079537432329687, 0.0160654878216455188066902, + 0.0176506444778379567583926, 0.0349580018370856671072744, 0.0017616969113226077503498, + 0.0043762331129306719645911, 0.0079872301725372740754949, 0.0211556317124290127473785, + 0.0023057907016680312395351), tolerance = tolerance)) > > pmhat_in_offset_mvord <- marginal_predict(res_train2) > mvord:::check(all.equal( + unname(c(pmhat_in_offset_mvord[1:5, 1:2])), c(0.05687801748367393717132, 0.43314394646591181103901, 0.63410255320208053220199, + 0.11659868247577427624595, 0.68287284549406501721336, 0.07832694375364634975512, + 0.55065457504773962504885, 0.57624707519535223188001, 0.62708167599782682621878, + 0.62089949532970512002805), tolerance = tolerance)) > > > ## predict with newdata > phat_new_offset_mvord <- predict(res_train2, newdata = dat_test) > mvord:::check(all.equal( + unname(phat_new_offset_mvord), c(0.53788568792038271570988, 0.47226651862936941395077, 0.59736906923858912321634, + 0.51730764069679358030385, 0.18750767549430336078586, 0.01993645212535863353587, + 0.05609361865051418205574, 0.15695387935814331115125, 0.57079407463923570453801, + 0.46447453401393973271283), tolerance = tolerance)) > > pjhat_new_offset_mvord <- joint_probabilities(res_train2, response.cat = c(1,2), newdata = dat_test) > mvord:::check(all.equal( + unname(pjhat_new_offset_mvord), c(0.0839907843603347747940546, 0.0053737551172480801930931, 0.1657125298469289131908511, + 0.0095039724989203211436006, 0.1110141336202859763115924, 0.0308028085793187883512090, + 0.0560936429791563384572584, 0.0000659201973896017534571, 0.0161492922521473680763648, + 0.0556615607647627450016437), tolerance = tolerance)) > > pmhat_new_offset_mvord <- marginal_predict(res_train2, newdata = dat_test) > mvord:::check(all.equal( + unname(c(pmhat_new_offset_mvord[1:5,1:2])), c( + 0.6410297683731078777214, 0.5530240404054447278526, 0.7630836817892157064591, + 0.6103504462453783752096, 0.2985228989322880055468, 0.6269027239883980806567, + 0.5461645309760024824541, 0.6047542519330054711091, 0.5786243958547179211394, + 0.2011002687374978670221), tolerance = tolerance)) > > > > proc.time() user system elapsed 35.5700000000000002842171 0.7099999999999999644729 36.2800000000000011368684