[R-es] aov
Carlos Ortega
cof en qualityexcellence.es
Vie Mayo 27 15:00:37 CEST 2016
Hola,
Prueba con el código que te adjunto:
- Son tres funciones.
- El cambio está en la última función (*mysjp.aov1*). Tan sólo hay que
cambiar un par de parámetros en la función "annotate()" y ya sale centrado,
un poco más pequeño y en color rojo...
- Al final aparece la ejecución con tus propios datos.
#----------------------------
get_var_name <- function(x) {
# remove "data frame name"
dollar_pos <- regexpr("$", x, fixed = T)[1]
if (dollar_pos != -1)
x <-
substr(x, start = dollar_pos + 1, stop = nchar(x))
return(x)
}
get_p_stars <- function(pval) {
pan <- ""
if (is.na(pval))
pan <- ""
else if (pval < 0.001)
pan <- "***"
else if (pval < 0.01)
pan <- "**"
else if (pval < 0.05)
pan <- "*"
return(invisible(pan))
}
mysjp.aov1 <- function (var.dep, var.grp, meansums = FALSE, title = NULL,
axis.labels = NULL,
rev.order = FALSE, string.interc = "(Intercept)", axis.title =
"",
axis.lim = NULL, geom.colors = c("#3366a0", "#aa3333"),
geom.size = 3,
wrap.title = 50, wrap.labels = 25, grid.breaks = NULL,
show.values = TRUE,
digits = 2, y.offset = 0.1, show.p = TRUE, show.summary =
FALSE,
prnt.plot = TRUE)
{
var.grp.name <- get_var_name(deparse(substitute(var.grp)))
var.dep.name <- get_var_name(deparse(substitute(var.dep)))
if (is.null(axis.labels))
axis.labels <- sjmisc::get_labels(var.grp, attr.only = F,
include.values = NULL,
include.non.labelled = T)
if (is.null(axis.title))
axis.title <- sjmisc::get_label(var.dep, def.value = var.dep.name)
if (is.null(title)) {
t1 <- sjmisc::get_label(var.grp, def.value = var.grp.name)
t2 <- sjmisc::get_label(var.dep, def.value = var.dep.name)
if (!is.null(t1) && !is.null(t2))
title <- paste0(t1, " by ", t2)
}
if (!is.null(axis.labels) && axis.labels == "")
axis.labels <- NULL
if (!is.null(axis.title) && axis.title == "")
axis.title <- NULL
if (!is.null(title) && title == "")
title <- NULL
if (!is.null(axis.labels)) {
axis.labels[1] <- paste(axis.labels[1], string.interc)
}
if (!is.factor(var.grp))
var.grp <- as.factor(var.grp)
if (!is.null(title))
title <- sjmisc::word_wrap(title, wrap.title)
if (!is.null(axis.title))
axis.title <- sjmisc::word_wrap(axis.title, wrap.title)
if (!is.null(axis.labels))
axis.labels <- sjmisc::word_wrap(axis.labels, wrap.labels)
fit <- stats::aov(var.dep ~ var.grp)
means <- stats::summary.lm(fit)$coefficients[, 1]
means.p <- stats::summary.lm(fit)$coefficients[, 4]
means.lci <- stats::confint(fit)[, 1]
means.uci <- stats::confint(fit)[, 2]
if (meansums) {
for (i in 2:length(means)) {
means[i] <- means[i] + means[1]
means.lci[i] <- means.lci[i] + means[1]
means.uci[i] <- means.uci[i] + means[1]
}
}
if (show.summary) {
ss <- summary(fit)[[1]]["Sum Sq"]
r2 <- stats::summary.lm(fit)$r.squared
r2.adj <- stats::summary.lm(fit)$adj.r.squared
fstat <- stats::summary.lm(fit)$fstatistic[1]
pval <- summary(fit)[[1]]["Pr(>F)"][1, 1]
pan <- get_p_stars(pval)
modsum <- as.character(as.expression(substitute(italic(SS[B]) ==
ssb * "," ~
~italic(SS[W]) == ssw * "," ~ ~R^2 ==
mr2 * "," ~ ~"adj."
* R^2 == ar2 * "," ~ ~"F" ==
f * panval,
list(ssb = sprintf("%.2f", ss[1, ]),
ssw = sprintf("%.2f", ss[2, ]), mr2 = sprintf("%.3f",
r2), ar2 = sprintf("%.3f",
r2.adj), f = sprintf("%.2f",
fstat), panval = pan))))
}
ps <- c(round(means, digits))
if (!show.values)
ps <- rep("", length(ps))
if (show.p) {
for (i in 1:length(means.p)) {
ps[i] <- sjmisc::trim(paste(ps[i], get_p_stars(means.p[i])))
}
}
if (rev.order)
catorder <- c(length(means):1)
else catorder <- c(1:length(means))
df <- data.frame(means, means.lci, means.uci, means.p, ps,
catorder)
if (is.null(axis.labels))
axis.labels <- row.names(df)
axis.labels <- axis.labels[catorder]
names(df) <- c("means", "lower", "upper", "p", "pv", "xv")
df$means <- sjmisc::to_value(df$means, keep.labels = F)
df$lower <- sjmisc::to_value(df$lower, keep.labels = F)
df$upper <- sjmisc::to_value(df$upper, keep.labels = F)
df$p <- sjmisc::to_value(df$p, keep.labels = F)
df$pv <- as.character(df$pv)
df <- cbind(df, geocol = ifelse(df$means >= 0, geom.colors[1],
geom.colors[2]))
if (is.null(axis.lim)) {
maxval <- max(df$upper)
minval <- min(df$lower)
if (maxval > 0)
limfac <- ifelse(abs(maxval) < 5, 5, 10)
else limfac <- ifelse(abs(minval) < 5, 5, 10)
upper_lim <- ifelse(maxval == 0, 0, limfac * ceiling((maxval +
1)/limfac))
lower_lim <- ifelse(minval == 0, 0, limfac * floor(minval/limfac))
}
else {
lower_lim <- axis.lim[1]
upper_lim <- axis.lim[2]
}
if (is.null(grid.breaks))
ticks <- pretty(c(lower_lim, upper_lim))
else ticks <- c(seq(lower_lim, upper_lim, by = grid.breaks))
scaley <- scale_y_continuous(limits = c(lower_lim, upper_lim),
breaks = ticks, labels = ticks)
anovaplot <- ggplot(df, aes(y = means, x = as.factor(xv))) +
geom_point(size = geom.size, colour = df$geocol) +
geom_errorbar(aes(ymin = lower,
ymax = upper), colour = df$geocol, width = 0) + geom_text(aes(label = pv,
y = means),
nudge_x = y.offset, show.legend = FALSE) +
scaley + scale_x_discrete(labels = axis.labels, limits =
1:length(axis.labels)) +
labs(title = title, x = NULL, y = axis.title) + coord_flip()
#---------------Changes: size = 3 / hjust = 1.2 / col = "red"
if (show.summary) {
anovaplot <- anovaplot + annotate("text", label = modsum,
parse = TRUE, x = -Inf, y = Inf,
hjust = 1.2,
vjust = "bottom", size = 3, col =
"red")
}
if (prnt.plot)
graphics::plot(anovaplot)
df <- dplyr::add_rownames(df)
colnames(df) <- c("term", "estimate", "conf.low", "conf.high",
"p.value", "p.string", "xpos", "geom.color")
invisible(structure(class = "sjpaov1", list(plot = anovaplot,
data = df)))
}
mysjp.aov1(
depVar, grpVar, meansums = TRUE, title = NULL,
#axisLabels.y = NULL,
#showAxisLabels.y = TRUE,
#axisTitle.x = "", axisLimits = NULL,
#breakTitleAt = 25, breakLabelsAt = 25,
#expand.grid = TRUE,
rev.order = FALSE,
string.interc = "(Intercept)",
geom.colors = c("#3366a0", "#aa3333"), geom.size = 3,
grid.breaks = NULL,
show.values = TRUE,
digits = 2, y.offset = 0.1, show.p = TRUE,
show.summary = TRUE, prnt.plot = TRUE
)
#----------------------------
Saludos,
Carlos Ortega
www.qualityexcellence.es
El 26 de mayo de 2016, 0:46, josebetancourt.cmw <
josebetancourt.cmw en infomed.sld.cu> escribió:
>
>
> Estimados
>
>
>
> ¿Cómo se podría modificar el script o su source para lograr que en la
> salida se vean de manera clara los números? En este ejemplo quedan
> parcialmente afuera. Adjunto datos y script
>
> Saludos
>
> José
>
>
>
> ##########################################################################
>
>
>
> rm(list = ls())
>
> setwd("D:/Public/Documents/R/r_epidemiología/")
>
>
>
> data<-read.csv("dengueutiuti10.csv",header=TRUE, sep=";", dec=",")
>
> head(data)
>
> depVar=data$plt
>
> grpVar=data$bleed
>
>
>
> library(sjPlot)
>
>
>
> sjp.aov1(depVar, grpVar, meansums = TRUE, title = NULL,
>
> axisLabels.y = NULL, reverseOrder = FALSE,
>
> stringIntercept = "(Intercept)", showAxisLabels.y = TRUE,
>
> axisTitle.x = "", axisLimits = NULL, geom.colors = c("#3366a0",
>
> "#aa3333"), geom.size = 3, breakTitleAt = 25, breakLabelsAt = 25,
>
> gridBreaksAt = NULL, expand.grid = TRUE, showValueLabels = TRUE,
>
> labelDigits = 2, y.offset = 0.1, showPValueLabels = TRUE,
>
> showModelSummary = TRUE, printPlot = TRUE)
>
>
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
--
Saludos,
Carlos Ortega
www.qualityexcellence.es
[[alternative HTML version deleted]]
Más información sobre la lista de distribución R-help-es