Skip to content

Commit

Permalink
revise
Browse files Browse the repository at this point in the history
  • Loading branch information
@strengejacke
strengejacke committed May 17, 2019
1 parent 0b2b8c3 commit fcdd9a9
Show file tree
Hide file tree
Showing 18 changed files with 56 additions and 3,127 deletions.
10 changes: 2 additions & 8 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,13 @@
S3method(AIC,svyglm.nb)
S3method(as.data.frame,sj_resample)
S3method(as.integer,sj_resample)
S3method(cramer,formula)
S3method(cramer,table)
S3method(deviance,svyglm.nb)
S3method(family,svyglm.nb)
S3method(formula,svyglm.nb)
S3method(getME,brmsfit)
S3method(getME,stanreg)
S3method(knit_print,tidy_stan)
S3method(mcse,brmsfit)
S3method(mcse,stanmvreg)
S3method(mcse,stanreg)
Expand Down Expand Up @@ -77,8 +78,6 @@ S3method(se,lm)
S3method(se,lmerMod)
S3method(se,merModLmerTest)
S3method(se,nlmerMod)
S3method(se,sj_icc)
S3method(se,sj_icc_merMod)
S3method(se,stanfit)
S3method(se,stanreg)
S3method(se,svyglm.nb)
Expand Down Expand Up @@ -243,13 +242,9 @@ importFrom(insight,get_response)
importFrom(insight,link_inverse)
importFrom(insight,model_info)
importFrom(insight,print_color)
importFrom(knitr,asis_output)
importFrom(knitr,knit_print)
importFrom(lme4,findbars)
importFrom(lme4,fixef)
importFrom(lme4,getME)
importFrom(lme4,glmer)
importFrom(lme4,lmer)
importFrom(lme4,ngrps)
importFrom(lme4,ranef)
importFrom(magrittr,"%>%")
Expand Down Expand Up @@ -364,4 +359,3 @@ importFrom(stats,xtabs)
importFrom(tidyr,gather)
importFrom(tidyr,nest)
importFrom(tidyr,unnest)
importFrom(utils,txtProgressBar)
21 changes: 7 additions & 14 deletions R/S3-methods.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,13 +188,6 @@ print.tidy_stan <- function(x, ...) {
}


#' @importFrom knitr knit_print asis_output
#' @export
knit_print.tidy_stan <- function(input, ...) {
knitr::asis_output(.print_tidy_stan(input, .color = FALSE, ...))
}


.print_tidy_stan <- function(x, .color, ...) {
if (.color)
insight::print_color("\n# Summary Statistics of Stan-Model\n\n", "blue")
Expand Down Expand Up @@ -333,7 +326,7 @@ knit_print.tidy_stan <- function(input, ...) {
} else {

if (.color)
insight::print_color(sprintf("## Response%s: %s\n\n", resp.string, insight::print_color(resp, "red")), "blue")
insight::print_color(sprintf("## Response%s: %s\n\n", resp.string, resp), "blue")
else
cat(sprintf("## Response%s: %s\n\n", resp.string, resp))

Expand Down Expand Up @@ -488,12 +481,12 @@ print_stan_ranef <- function(x, zeroinf = FALSE, .color) {

if (!zeroinf) {
if (.color)
insight::print_color(sprintf("## Random effect %s\n\n", insight::print_color(r, "red")), "blue")
insight::print_color(sprintf("## Random effect %s\n\n", r), "blue")
else
cat(sprintf("## Random effect %s\n\n", r))
} else {
if (.color)
insight::print_color(sprintf("## Conditional Model: Random effect %s\n\n", insight::print_color(r, "red")), "blue")
insight::print_color(sprintf("## Conditional Model: Random effect %s\n\n", r), "blue")
else
cat(sprintf("## Conditional Model: Random effect %s\n\n", r))
}
Expand Down Expand Up @@ -548,7 +541,7 @@ print_stan_mv_re <- function(x, resp, .color) {
if (!sjmisc::is_empty(fe)) {

if (.color)
insight::print_color(sprintf("## Fixed effects for response: %s\n\n", insight::print_color(resp, "red")), "blue")
insight::print_color(sprintf("## Fixed effects for response: %s\n\n", resp), "blue")
else
cat(sprintf("## Fixed effects for response: %s\n\n", resp))

Expand Down Expand Up @@ -577,8 +570,8 @@ print_stan_mv_re <- function(x, resp, .color) {

if (!sjmisc::is_empty(find.re)) {
if (.color) {
insight::print_color(sprintf("## Random effect %s", insight::print_color(resp, "red")), "blue")
insight::print_color(sprintf(" for response %s\n\n", insight::print_color(resp, "red")), "blue")
insight::print_color(sprintf("## Random effect %s", resp), "blue")
insight::print_color(sprintf(" for response %s\n\n", resp), "blue")
} else {
cat(sprintf("## Random effect %s", resp))
cat(sprintf(" for response %s\n\n", resp))
Expand Down Expand Up @@ -650,7 +643,7 @@ print_stan_zeroinf_ranef <- function(x, .color) {
for (r in re) {

if (.color)
insight::print_color(sprintf("## Zero-Inflated Model: Random effect %s\n\n", insight::print_color(r, "red")), "blue")
insight::print_color(sprintf("## Zero-Inflated Model: Random effect %s\n\n", r), "blue")
else
cat(sprintf("## Zero-Inflated Model: Random effect %s\n\n", r))

Expand Down
20 changes: 19 additions & 1 deletion R/phi.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,25 @@ phi <- function(tab) {

#' @rdname xtab_statistics
#' @export
cramer <- function(tab) {
cramer <- function(tab, ...) {
UseMethod("cramer")
}


#' @export
cramer.table <- function(tab) {
.cramer(tab)
}

#' @rdname xtab_statistics
#' @export
cramer.formula <- function(formula, data, ...) {
terms <- all.vars(formula)
tab <- table(data[[terms[1]]], data[[terms[2]]])
.cramer(tab)
}

.cramer <- function(tab) {
# convert to flat table
if (!inherits(tab, "ftable")) tab <- stats::ftable(tab)
sqrt(phi(tab)^2 / min(dim(tab) - 1))
Expand Down
104 changes: 0 additions & 104 deletions R/se.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,18 +259,6 @@ se.svyglm.nb <- function(x, ...) {
}


#' @rdname se
#' @export
se.sj_icc_merMod <- function(x, nsim = 100, ...) {
std_e_icc(x, nsim, adjusted = FALSE)
}

#' @rdname se
#' @export
se.sj_icc <- function(x, nsim = 100, ...) {
std_e_icc(x, nsim, adjusted = TRUE)
}

#' @export
se.glmerMod <- function(x, ...) {
std_merMod(x)
Expand Down Expand Up @@ -360,98 +348,6 @@ std_merMod <- function(fit) {
se.merMod
}


std_e_icc <- function(x, nsim, adjusted = FALSE) {
# check whether model is still in environment?
obj.name <- attr(x, ".obj.name", exact = T)
if (!exists(obj.name, envir = globalenv()))
stop(sprintf("Can't find merMod-object `%s` (that was used to compute the ICC) in the environment.", obj.name), call. = F)

# get object, see whether formulas match
fitted.model <- globalenv()[[obj.name]]
model.formula <- attr(x, "formula", exact = T)

if (!identical(model.formula, formula(fitted.model)))
stop(sprintf("merMod-object `%s` was fitted with a different formula than ICC-model", obj.name), call. = F)

# get model family, we may have glmer
model.family <- attr(x, "family", exact = T)

# check for all required arguments
if (missing(nsim) || is.null(nsim)) nsim <- 100

# get ICC, and compute bootstrapped SE, than return both
bstr <- bootstr_icc_se(
model_frame(fitted.model, fe.only = FALSE),
nsim,
model.formula,
model.family,
adjusted
)

# get adjusted ICC only
if (adjusted) x <- x[[1]]

# now compute SE and p-values for the bootstrapped ICC
res <- data_frame(
model = obj.name,
icc = as.vector(x),
std.err = boot_se(bstr)[["std.err"]],
p.value = boot_p(bstr)[["p.value"]]
)

structure(
class = "sj_se_icc",
list(result = res, bootstrap_data = bstr, adjusted = adjusted)
)
}



#' @importFrom dplyr mutate
#' @importFrom lme4 lmer glmer
#' @importFrom utils txtProgressBar
#' @importFrom purrr map map_dbl
#' @importFrom rlang .data
bootstr_icc_se <- function(dd, nsim, formula, model.family, adjusted) {
# create progress bar
pb <- utils::txtProgressBar(min = 1, max = nsim, style = 3)

# generate bootstraps
dummy <- dd %>%
bootstrap(nsim) %>%
dplyr::mutate(
models = purrr::map(.data$strap, function(x) {
# update progress bar
utils::setTxtProgressBar(pb, x$resample.id)
# check model family, then compute mixed model
if (model.family == "gaussian")
lme4::lmer(formula, data = x)
else
lme4::glmer(formula, data = x, family = model.family)
}),
# compute ICC(s) for each "bootstrapped" regression
icc = suppressMessages(purrr::map(.data$models, ~ icc(.x, adjusted = adjusted)))
)

# we may have more than one random term in the model, so we might have
# multiple ICC values. In this case, we need to split the multiple icc-values
# into multiple columns, i.e. one column per ICC value

if (adjusted) {
icc_data <-
as.data.frame(matrix(unlist(purrr::map(dummy$icc, ~ .x[[1]])), nrow = nrow(dummy)))
} else {
icc_data <-
as.data.frame(matrix(unlist(purrr::map(dummy$icc, ~ .x)), nrow = nrow(dummy)))
}

# close progresss bar
close(pb)

icc_data
}

## TODO se for poisson etc

# # for poisson family, we need a different delta method approach
Expand Down
30 changes: 16 additions & 14 deletions R/tidy_stan.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,11 +100,11 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(

# compute HDI / ci
if (!is.null(trans)) {
out.hdi <- bayestestR::ci(x, ci = prob, effects = type)
colnames(out.hdi) <- c("term", "ci.lvl", "ci.low", "ci.high")
out.hdi <- bayestestR::ci(x, ci = prob, effects = type, component = "all")
colnames(out.hdi)[1:4] <- c("term", "ci.lvl", "ci.low", "ci.high")
} else {
out.hdi <- bayestestR::hdi(x, ci = prob, effects = type)
colnames(out.hdi) <- c("term", "ci.lvl", "hdi.low", "hdi.high")
out.hdi <- bayestestR::hdi(x, ci = prob, effects = type, component = "all")
colnames(out.hdi)[1:4] <- c("term", "ci.lvl", "hdi.low", "hdi.high")
}

# transform data frame for multiple ci-levels
Expand All @@ -127,9 +127,9 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(

if (inherits(x, "brmsfit")) {
cnames <- make.names(names(nr))
keep <- cnames %in% out.hdi$term
keep <- cnames %in% names(mod.dat)
} else {
keep <- names(nr) %in% out.hdi$term
keep <- names(nr) %in% names(mod.dat)
}

nr <- nr[keep]
Expand All @@ -143,9 +143,9 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(

if (inherits(x, "brmsfit")) {
cnames <- make.names(names(rh))
keep <- cnames %in% out.hdi$term
keep <- cnames %in% names(mod.dat)
} else {
keep <- names(rh) %in% out.hdi$term
keep <- names(rh) %in% names(mod.dat)
}

rh <- rh[keep]
Expand All @@ -160,7 +160,7 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(
}

se <- mcse(x, type = type)
se <- se[se$term %in% out.hdi$term, ]
se <- se[se$term %in% names(mod.dat), ]


# transform estimate, if requested
Expand All @@ -179,12 +179,9 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(
out <- data_frame(
term = names(est),
estimate = est,
std.error = purrr::map_dbl(mod.dat, stats::mad)
std.error = purrr::map_dbl(mod.dat, stats::mad),
out.hdi[, -1]
) %>%
dplyr::inner_join(
out.hdi,
by = "term"
) %>%
dplyr::inner_join(
ratio,
by = "term"
Expand Down Expand Up @@ -404,6 +401,9 @@ tidy_stan <- function(x, prob = .89, typical = "median", trans = NULL, type = c(
}


if (obj_has_name(out, "Component")) out[, "Component"] <- NULL
if (obj_has_name(out, "Group")) out[, "Group"] <- NULL

class(out) <- c("tidy_stan", class(out))

attr(out, "digits") <- digits
Expand Down Expand Up @@ -560,6 +560,8 @@ n_of_chains <- function(x) {
ratio <- ess/tss
ratio <- ratio[!names(ratio) %in% c("mean_PPD", "log-posterior")]
}

ratio
}


Expand Down
42 changes: 0 additions & 42 deletions inst/doc/anova-statistics.R
View file

This file was deleted.

Loading

0 comments on commit fcdd9a9

Please sign in to comment.