support stanmvreg in tidy_stan()

DESCRIPTION

@@ -66,6 +66,7 @@ Suggests:
     pROC,
     psych,
     scales,
+    splines,
     sjPlot,
     survey,
     rstan,

R/mcse.R

@@ -16,6 +16,9 @@ mcse.brmsfit <- function(x, type = c("fixed", "random", "all"), ...) {
 
 #' @export
 mcse.stanmvreg <- function(x, type = c("fixed", "random", "all"), ...) {
+  # check arguments
+  type <- match.arg(type)
+
   s <- summary(x)
   dat <- tibble::tibble(
     term = rownames(s),

R/pred_vars.R

@@ -11,7 +11,7 @@
 #'    model, returns the model frame for fixed effects only.
 #' @param multi.resp Logical, if \code{TRUE} and model is a multivariate response
 #'    model from a \code{brmsfit} object or of class \code{stanmvreg}, then a
-#'    list of values for each regression is returned.
+#'    list of values (one for each regression) is returned.
 #'
 #' @return For \code{pred_vars()} and \code{resp_var()}, the name(s) of the
 #'    response or predictor variables from \code{x} as character vector.

R/tidy_stan.R

@@ -97,27 +97,36 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(
   # compute HDI
   out.hdi <- hdi(x, prob = prob, trans = trans, type = type)
 
-  # we need names of elements, for correct removal
+  # get statistics
   nr <- bayesplot::neff_ratio(x)
 
+  # we need names of elements, for correct removal
+
   if (inherits(x, "brmsfit")) {
     cnames <- make.names(names(nr))
     keep <- cnames %in% out.hdi$term
   } else {
-    keep <- 1:nrow(out.hdi)
+    keep <- names(nr) %in% out.hdi$term
   }
 
-
-  # compute additional statistics, like point estimate, standard errors etc.
-
   nr <- nr[keep]
   ratio <- data.frame(
     term = names(nr),
     ratio = nr,
     stringsAsFactors = FALSE
   )
 
-  rh <- bayesplot::rhat(x)[keep]
+
+  rh <- bayesplot::rhat(x)
+
+  if (inherits(x, "brmsfit")) {
+    cnames <- make.names(names(rh))
+    keep <- cnames %in% out.hdi$term
+  } else {
+    keep <- names(rh) %in% out.hdi$term
+  }
+
+  rh <- rh[keep]
   rhat <- data.frame(
     term = names(rh),
     rhat = rh,
@@ -243,6 +252,9 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(
     }
 
 
+    ## TODO extract Sigma for stanmvreg random effects
+
+
     # find random slopes
 
     rs1 <- grep("b\\[(.*) (.*)\\]", out$term)
@@ -324,7 +336,32 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(
 
   }
 
-  ## TODO add support for multivariate response model for rstanarm
+
+  if (inherits(x, "stanmvreg")) {
+
+    # get response variables
+
+    responses <- resp_var(x)
+    resp.names <- names(responses)
+
+
+    # create "response-level" variable
+
+    out <- tibble::add_column(out, response = "", .before = 1)
+
+
+    # copy name of response into new character variable
+    # and remove response name from term name
+
+    for (i in 1:length(responses)) {
+      pattern <- paste0(resp.names[i], "|")
+      m <- tidyselect::starts_with(pattern, vars = out$term)
+      out$response[intersect(which(out$response == ""), m)] <- responses[i]
+      out$term <- gsub(pattern, "", out$term, fixed = TRUE)
+    }
+
+  }
+
 
   class(out) <- c("tidy_stan", class(out))
 

inst/doc/anova-statistics.html

@@ -12,7 +12,7 @@
 
 <meta name="author" content="Daniel Lüdecke" />
 
-<meta name="date" content="2018-07-08" />
+<meta name="date" content="2018-07-09" />
 
 <title>Statistics for Anova Tables</title>
 
@@ -70,7 +70,7 @@
 
 <h1 class="title toc-ignore">Statistics for Anova Tables</h1>
 <h4 class="author"><em>Daniel Lüdecke</em></h4>
-<h4 class="date"><em>2018-07-08</em></h4>
+<h4 class="date"><em>2018-07-09</em></h4>
 
 
 
@@ -194,9 +194,9 @@ <h1>Confidence Intervals</h1>
 <span class="co">#&gt; # A tibble: 3 x 4</span>
 <span class="co">#&gt;   term     partial.omegasq conf.low conf.high</span>
 <span class="co">#&gt;   &lt;chr&gt;              &lt;dbl&gt;    &lt;dbl&gt;     &lt;dbl&gt;</span>
-<span class="co">#&gt; 1 e42dep           0.278    0.228      0.338 </span>
-<span class="co">#&gt; 2 c172code         0.00547 -0.00610    0.0224</span>
-<span class="co">#&gt; 3 c160age          0.0649   0.0315     0.104</span></code></pre></div>
+<span class="co">#&gt; 1 e42dep           0.278    0.223      0.332 </span>
+<span class="co">#&gt; 2 c172code         0.00547 -0.00671    0.0221</span>
+<span class="co">#&gt; 3 c160age          0.0649   0.0348     0.0997</span></code></pre></div>
 </div>
 <div id="references" class="section level1">
 <h1>References</h1>

inst/doc/mixedmodels-statistics.html

@@ -12,7 +12,7 @@
 
 <meta name="author" content="Daniel Lüdecke" />
 
-<meta name="date" content="2018-07-08" />
+<meta name="date" content="2018-07-09" />
 
 <title>Statistics for Mixed Effects Models</title>
 
@@ -70,7 +70,7 @@
 
 <h1 class="title toc-ignore">Statistics for Mixed Effects Models</h1>
 <h4 class="author"><em>Daniel Lüdecke</em></h4>
-<h4 class="date"><em>2018-07-08</em></h4>
+<h4 class="date"><em>2018-07-09</em></h4>
 
 
 
@@ -180,21 +180,21 @@ <h2>P-Values</h2>
 <div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Using the t-statistics for calculating the p-value</span>
 <span class="kw">p_value</span>(m2)
 <span class="co">#&gt;          term p.value std.error</span>
-<span class="co">#&gt; 1 (Intercept)       0     9.747</span>
-<span class="co">#&gt; 2        Days       0     0.804</span>
+<span class="co">#&gt; 1 (Intercept)       0     9.760</span>
+<span class="co">#&gt; 2        Days       0     0.805</span>
 
 <span class="co"># p-values based on conditional F-tests with </span>
 <span class="co"># Kenward-Roger approximation for the degrees of freedom</span>
 <span class="kw">p_value</span>(m2, <span class="dt">p.kr =</span> <span class="ot">TRUE</span>)
 <span class="co">#&gt;          term p.value std.error</span>
-<span class="co">#&gt; 1 (Intercept)       0     9.769</span>
+<span class="co">#&gt; 1 (Intercept)       0     9.782</span>
 <span class="co">#&gt; 2        Days       0     0.814</span></code></pre></div>
 <p>To see more details on the degrees of freedom when using Kenward-Roger approximation, use the <code>summary()</code>-method:</p>
 <div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">pv &lt;-<span class="st"> </span><span class="kw">p_value</span>(m2, <span class="dt">p.kr =</span> <span class="ot">TRUE</span>)
 <span class="kw">summary</span>(pv)
 <span class="co">#&gt;          term p.value std.error      df statistic</span>
-<span class="co">#&gt; 1 (Intercept)       0     9.769  22.754    25.734</span>
-<span class="co">#&gt; 2        Days       0     0.814 160.934    12.865</span></code></pre></div>
+<span class="co">#&gt; 1 (Intercept)       0     9.782  23.021    25.719</span>
+<span class="co">#&gt; 2        Days       0     0.814 160.966    12.826</span></code></pre></div>
 </div>
 <div id="random-effect-variances" class="section level2">
 <h2>Random Effect Variances</h2>
@@ -246,8 +246,8 @@ <h2>Intraclass-Correlation Coefficient</h2>
 <span class="co">#&gt;  Family: gaussian (identity)</span>
 <span class="co">#&gt; Formula: Reaction ~ Days + (1 | mygrp) + (1 | Subject)</span>
 <span class="co">#&gt; </span>
-<span class="co">#&gt;     ICC (mygrp): 0.000000</span>
-<span class="co">#&gt;   ICC (Subject): 0.589309</span></code></pre></div>
+<span class="co">#&gt;     ICC (mygrp): 0.010332</span>
+<span class="co">#&gt;   ICC (Subject): 0.588298</span></code></pre></div>
 </div>
 </div>
 <div id="references" class="section level1">