added badges

docs/articles/application.html

@@ -265,6 +265,105 @@
 <a class="sourceLine" id="cb5-31" title="31">legend &lt;-<span class="st"> </span><span class="kw"><a href="https://rdrr.io/r/base/lapply.html">sapply</a></span>(sigma,<span class="cf">function</span>(x) <span class="kw"><a href="https://rdrr.io/r/base/expression.html">as.expression</a></span>(<span class="kw"><a href="https://rdrr.io/r/base/bquote.html">bquote</a></span>(sigma <span class="op">==</span><span class="st"> </span>.(x))))</a>
 <a class="sourceLine" id="cb5-32" title="32">graphics<span class="op">::</span><span class="kw"><a href="https://rdrr.io/r/graphics/legend.html">legend</a></span>(<span class="dt">x=</span><span class="st">"topleft"</span>,<span class="dt">legend=</span>legend,<span class="dt">lty=</span>lty,<span class="dt">bty=</span><span class="st">"n"</span>,<span class="dt">lwd=</span>lwd)</a>
 <a class="sourceLine" id="cb5-33" title="33">grDevices<span class="op">::</span><span class="kw"><a href="https://rdrr.io/r/grDevices/dev.html">dev.off</a></span>()</a></code></pre></div>
+<!--
+# ordinal regression
+
+It is possible to compare combined regression with ordinal regression by slightly modifying the function cv.cornet. Add the argument "ordinal = TRUE/FALSE" and two code chunks (see below). Such a comparison does not make much sense for ordinal response variables with many categories.
+
+
+```r
+cv.cornet <- function (y, cutoff, X, alpha = 1, nfolds.ext = 5, nfolds.int = 10, foldid.ext = NULL, foldid.int = NULL, type.measure = "deviance", ordinal = FALSE,...) {
+    z <- 1 * (y > cutoff)
+    if (is.null(foldid.ext)) {
+        foldid.ext <- palasso:::.folds(y = z, nfolds = nfolds.ext)
+    } else {
+        nfolds.ext <- length(unique(foldid.ext))
+    }
+    cols <- c("intercept", "binomial", "combined")
+    ### trial start ###
+    if(ordinal){cols <- c(cols,"ordinal")}
+    ### trial end ###
+    pred <- matrix(data = NA, nrow = length(y), ncol = length(cols), 
+        dimnames = list(NULL, cols))
+    for (i in seq_len(nfolds.ext)) {
+        y0 <- y[foldid.ext != i]
+        z0 <- z[foldid.ext != i]
+        X0 <- X[foldid.ext != i, ]
+        X1 <- X[foldid.ext == i, ]
+        if (is.null(foldid.int)) {
+            foldid <- palasso:::.folds(y = z0, nfolds = nfolds.int)
+        } else {
+            foldid <- foldid.int[foldid.ext != i]
+        }
+        fit <- cornet::cornet(y = y0, cutoff = cutoff, X = X0, 
+            alpha = alpha, type.measure = type.measure, foldid = foldid, 
+            ...)
+        tryCatch(expr = plot.cornet(fit), error = function(x) NULL)
+        temp <- predict.cornet(fit, newx = X1)
+        if (any(temp < 0 | temp > 1)) {
+            stop("Outside unit interval.", call. = FALSE)
+        }
+        model <- colnames(temp)
+        for (j in seq_along(model)) {
+            pred[foldid.ext == i, model[j]] <- temp[[model[j]]]
+        }
+        if(ordinal){
+          ### trial start ###
+          #browser()
+          y0_ord <- as.factor(y0)
+          fit <- ordinalNet::ordinalNet(x=X0,y=y0_ord,alpha=alpha)
+          pred_ord <- predict(fit,newx=X1)
+          above <- as.numeric(levels(y0_ord))>cutoff
+          pred[foldid.ext ==i,"ordinal"] <- rowSums(pred_ord[,above])
+        } ### trial end ###
+    }
+    type <- c("deviance", "class", "mse", "mae", "auc")
+    loss <- lapply(X = type, FUN = function(x) palasso:::.loss(y = z, 
+        fit = pred, family = "binomial", type.measure = x, 
+        foldid = foldid.ext)[[1]])
+    names(loss) <- type
+    loss <- lapply(loss, function(x) signif(x, digits = 6))
+    return(loss)
+}
+predict.cornet <- cornet:::predict.cornet
+
+loss <- list()
+for(i in seq_along(Y)){
+  load("data/processed_data.RData",verbose=TRUE)
+  cond <- !is.na(Y[[i]])
+  loss[[i]] <- cv.cornet(y=Y[[i]][cond],X=x[[1]][cond,],cutoff=25.5,ordinal=TRUE)
+}
+sapply(loss,function(x) x$deviance)
+
+
+### Ordinal regression with two different packages.
+#
+# load("data/processed_data.RData",verbose=TRUE)
+# cond <- !is.na(Y$moca.V04)
+# y <- as.factor(x=Y$moca.V04[cond])
+# x <- x[[1]][cond,]
+#
+# # ordinalNet
+# model <- ordinalNet::ordinalNet(x=x,y=y)
+# y_hat <- predict(model)
+# below <- as.numeric(levels(y))<=25.5
+# above <- as.numeric(levels(y))>25.5
+# b1 <- rowSums(y_hat[,below])
+# a1 <- rowSums(y_hat[,above])
+# 
+# # glmnetcr
+# model <- glmnetcr::glmnetcr(x=x,y=y)
+# s <- glmnetcr::select.glmnetcr(model)
+# y_hat <- predict(model)$probs[,,s]
+# below <- as.numeric(colnames(y_hat))<=25.5
+# above <- as.numeric(colnames(y_hat))>25.5
+# b2 <- rowSums(y_hat[,below])
+# a2 <- rowSums(y_hat[,above])
+# 
+# cor(a1,a2)
+# cor(b1,b2)
+```
+-->
 </div>
   </div>
 

docs/articles/simulation.html

@@ -132,10 +132,10 @@
 <a class="sourceLine" id="cb2-5" title="5"><span class="cf">for</span>(i <span class="cf">in</span> <span class="kw"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"deviance"</span>,<span class="st">"class"</span>,<span class="st">"mse"</span>,<span class="st">"mae"</span>,<span class="st">"auc"</span>)){</a>
 <a class="sourceLine" id="cb2-6" title="6">  fuse0[[i]] &lt;-<span class="st"> </span><span class="kw"><a href="https://rdrr.io/r/base/lapply.html">sapply</a></span>(ridge,<span class="cf">function</span>(x) (x[[i]][<span class="st">"combined"</span>]<span class="op">-</span>x[[i]][<span class="st">"binomial"</span>]))</a>
 <a class="sourceLine" id="cb2-7" title="7">  fuse1[[i]] &lt;-<span class="st"> </span><span class="kw"><a href="https://rdrr.io/r/base/lapply.html">sapply</a></span>(lasso,<span class="cf">function</span>(x) (x[[i]][<span class="st">"combined"</span>]<span class="op">-</span>x[[i]][<span class="st">"binomial"</span>]))</a>
-<a class="sourceLine" id="cb2-8" title="8">  <span class="cf">if</span>(i<span class="op">==</span><span class="st">"auc"</span>){fuse0[[i]] &lt;-<span class="st"> </span><span class="op">-</span>fuse0[[i]]; fuse1[[i]] &lt;-<span class="st"> </span><span class="op">-</span>fuse1[[i]]}</a>
+<a class="sourceLine" id="cb2-8" title="8">  <span class="co">#if(i=="auc"){fuse0[[i]] &lt;- -fuse0[[i]]; fuse1[[i]] &lt;- -fuse1[[i]]}</span></a>
 <a class="sourceLine" id="cb2-9" title="9">}</a>
 <a class="sourceLine" id="cb2-10" title="10"></a>
-<a class="sourceLine" id="cb2-11" title="11">grDevices<span class="op">::</span><span class="kw"><a href="https://rdrr.io/r/grDevices/pdf.html">pdf</a></span>(<span class="st">"manuscript/figure_BOX.pdf"</span>,<span class="dt">width=</span><span class="dv">6</span>,<span class="dt">height=</span><span class="dv">3</span>)</a>
+<a class="sourceLine" id="cb2-11" title="11">grDevices<span class="op">::</span><span class="kw"><a href="https://rdrr.io/r/grDevices/pdf.html">pdf</a></span>(<span class="st">"manuscript/figure_BOX.pdf"</span>,<span class="dt">width=</span><span class="dv">6</span>,<span class="dt">height=</span><span class="dv">4</span>)</a>
 <a class="sourceLine" id="cb2-12" title="12">graphics<span class="op">::</span><span class="kw"><a href="https://rdrr.io/r/graphics/par.html">par</a></span>(<span class="dt">mar=</span><span class="kw"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fl">1.9</span>,<span class="fl">1.9</span>,<span class="fl">0.1</span>,<span class="fl">0.1</span>))</a>
 <a class="sourceLine" id="cb2-13" title="13">graphics<span class="op">::</span><span class="kw"><a href="https://rdrr.io/r/graphics/frame.html">plot.new</a></span>()</a>
 <a class="sourceLine" id="cb2-14" title="14">ylim &lt;-<span class="st"> </span><span class="kw"><a href="https://rdrr.io/r/base/range.html">range</a></span>(<span class="kw"><a href="https://rdrr.io/r/base/unlist.html">unlist</a></span>(fuse0),<span class="kw"><a href="https://rdrr.io/r/base/unlist.html">unlist</a></span>(fuse1))</a>

vignettes/application.Rmd

@@ -198,3 +198,103 @@ legend <- sapply(sigma,function(x) as.expression(bquote(sigma == .(x))))
 graphics::legend(x="topleft",legend=legend,lty=lty,bty="n",lwd=lwd)
 grDevices::dev.off()
 ```
+
+
+<!--
+# ordinal regression
+
+It is possible to compare combined regression with ordinal regression by slightly modifying the function cv.cornet. Add the argument "ordinal = TRUE/FALSE" and two code chunks (see below). Such a comparison does not make much sense for ordinal response variables with many categories.
+
+```{r ordinal,eval=FALSE}
+cv.cornet <- function (y, cutoff, X, alpha = 1, nfolds.ext = 5, nfolds.int = 10, foldid.ext = NULL, foldid.int = NULL, type.measure = "deviance", ordinal = FALSE,...) {
+    z <- 1 * (y > cutoff)
+    if (is.null(foldid.ext)) {
+        foldid.ext <- palasso:::.folds(y = z, nfolds = nfolds.ext)
+    } else {
+        nfolds.ext <- length(unique(foldid.ext))
+    }
+    cols <- c("intercept", "binomial", "combined")
+    ### trial start ###
+    if(ordinal){cols <- c(cols,"ordinal")}
+    ### trial end ###
+    pred <- matrix(data = NA, nrow = length(y), ncol = length(cols), 
+        dimnames = list(NULL, cols))
+    for (i in seq_len(nfolds.ext)) {
+        y0 <- y[foldid.ext != i]
+        z0 <- z[foldid.ext != i]
+        X0 <- X[foldid.ext != i, ]
+        X1 <- X[foldid.ext == i, ]
+        if (is.null(foldid.int)) {
+            foldid <- palasso:::.folds(y = z0, nfolds = nfolds.int)
+        } else {
+            foldid <- foldid.int[foldid.ext != i]
+        }
+        fit <- cornet::cornet(y = y0, cutoff = cutoff, X = X0, 
+            alpha = alpha, type.measure = type.measure, foldid = foldid, 
+            ...)
+        tryCatch(expr = plot.cornet(fit), error = function(x) NULL)
+        temp <- predict.cornet(fit, newx = X1)
+        if (any(temp < 0 | temp > 1)) {
+            stop("Outside unit interval.", call. = FALSE)
+        }
+        model <- colnames(temp)
+        for (j in seq_along(model)) {
+            pred[foldid.ext == i, model[j]] <- temp[[model[j]]]
+        }
+        if(ordinal){
+          ### trial start ###
+          #browser()
+          y0_ord <- as.factor(y0)
+          fit <- ordinalNet::ordinalNet(x=X0,y=y0_ord,alpha=alpha)
+          pred_ord <- predict(fit,newx=X1)
+          above <- as.numeric(levels(y0_ord))>cutoff
+          pred[foldid.ext ==i,"ordinal"] <- rowSums(pred_ord[,above])
+        } ### trial end ###
+    }
+    type <- c("deviance", "class", "mse", "mae", "auc")
+    loss <- lapply(X = type, FUN = function(x) palasso:::.loss(y = z, 
+        fit = pred, family = "binomial", type.measure = x, 
+        foldid = foldid.ext)[[1]])
+    names(loss) <- type
+    loss <- lapply(loss, function(x) signif(x, digits = 6))
+    return(loss)
+}
+predict.cornet <- cornet:::predict.cornet
+
+loss <- list()
+for(i in seq_along(Y)){
+  load("data/processed_data.RData",verbose=TRUE)
+  cond <- !is.na(Y[[i]])
+  loss[[i]] <- cv.cornet(y=Y[[i]][cond],X=x[[1]][cond,],cutoff=25.5,ordinal=TRUE)
+}
+sapply(loss,function(x) x$deviance)
+
+
+### Ordinal regression with two different packages.
+#
+# load("data/processed_data.RData",verbose=TRUE)
+# cond <- !is.na(Y$moca.V04)
+# y <- as.factor(x=Y$moca.V04[cond])
+# x <- x[[1]][cond,]
+#
+# # ordinalNet
+# model <- ordinalNet::ordinalNet(x=x,y=y)
+# y_hat <- predict(model)
+# below <- as.numeric(levels(y))<=25.5
+# above <- as.numeric(levels(y))>25.5
+# b1 <- rowSums(y_hat[,below])
+# a1 <- rowSums(y_hat[,above])
+# 
+# # glmnetcr
+# model <- glmnetcr::glmnetcr(x=x,y=y)
+# s <- glmnetcr::select.glmnetcr(model)
+# y_hat <- predict(model)$probs[,,s]
+# below <- as.numeric(colnames(y_hat))<=25.5
+# above <- as.numeric(colnames(y_hat))>25.5
+# b2 <- rowSums(y_hat[,below])
+# a2 <- rowSums(y_hat[,above])
+# 
+# cor(a1,a2)
+# cor(b1,b2)
+```
+-->

vignettes/simulation.Rmd

@@ -58,10 +58,10 @@ fuse0 <- fuse1 <- list()
 for(i in c("deviance","class","mse","mae","auc")){
   fuse0[[i]] <- sapply(ridge,function(x) (x[[i]]["combined"]-x[[i]]["binomial"]))
   fuse1[[i]] <- sapply(lasso,function(x) (x[[i]]["combined"]-x[[i]]["binomial"]))
-  if(i=="auc"){fuse0[[i]] <- -fuse0[[i]]; fuse1[[i]] <- -fuse1[[i]]}
+  #if(i=="auc"){fuse0[[i]] <- -fuse0[[i]]; fuse1[[i]] <- -fuse1[[i]]}
 }
 
-grDevices::pdf("manuscript/figure_BOX.pdf",width=6,height=3)
+grDevices::pdf("manuscript/figure_BOX.pdf",width=6,height=4)
 graphics::par(mar=c(1.9,1.9,0.1,0.1))
 graphics::plot.new()
 ylim <- range(unlist(fuse0),unlist(fuse1))