[R] sm.options
Mark Difford
mark_difford at yahoo.co.uk
Fri Mar 6 09:22:11 CET 2009
Hi Viviana,
>> I am doing kernel density plots, and am trying to make the lines thicker.
You need to hack the code for sm.density.compare. See the code below. This
uses the same defaults as the original, but you can customize band colour,
line width, and so on using arguments to the function. The easiest way to
use it is to copy it into your environment space. The original code is
untouched and may be accessed as sm:::sm.density.compare().
Regards, Mark.
## Hack on sm.density.compare
## Allows me to change colour of the band, to set ylim, and change lwd (was
only for
## model="equal" option)
sm.density.compare <- function (x, group, h, model = "none", bandcol =
'cyan', lwd = par("lwd"), usePolyg = NULL, asp=NA,
xlab=opt$xlab, ylab=opt$ylab, ...)
{
if (!is.vector(x))
stop("sm.density.compare can handle only 1-d data")
opt <- sm.options(list(...))
sm:::replace.na(opt, ngrid, 50) ## These all changed
from replace.na() --> sm:::
sm:::replace.na(opt, display, "line")
sm:::replace.na(opt, xlab, deparse(substitute(x)))
sm:::replace.na(opt, ylab, "Density")
sm:::replace.na(opt, xlim, c(min(x) - diff(range(x))/4, max(x) +
diff(range(x))/4))
sm:::replace.na(opt, eval.points, seq(opt$xlim[1], opt$xlim[2],
length = opt$ngrid))
if (is.na(opt$band)) {
if (model == "none")
opt$band <- FALSE
else opt$band <- TRUE
}
if ((model == "none") && opt$band)
opt$band <- FALSE
band <- opt$band
ngrid <- opt$ngrid
xlim <- opt$xlim
nboot <- opt$nboot
y <- x
if (is.na(opt$test)) {
if (model == "none")
opt$test <- FALSE
else opt$test <- TRUE
}
if ((model == "none") && opt$test)
opt$test <- FALSE
test <- opt$test
if (opt$display %in% "none")
band <- FALSE
fact <- factor(group)
fact.levels <- levels(fact)
nlev <- length(fact.levels)
ni <- table(fact)
if (band & (nlev > 2)) {
cat("Reference band available to compare two groups only.",
"\n")
band <- FALSE
}
if (length(opt$lty) < nlev)
opt$lty <- 1:nlev
if (length(opt$col) < nlev)
opt$col <- 2:(nlev + 1)
if (missing(h))
h <- h.select(x, y = NA, group = group, ...)
opt$band <- band
opt$test <- test
estimate <- matrix(0, ncol = opt$ngrid, nrow = nlev)
se <- matrix(0, ncol = opt$ngrid, nrow = nlev)
for (i in 1:nlev) {
sm <- sm.density(y[fact == fact.levels[i]], h = h, display = "none",
eval.points = opt$eval.points)
estimate[i, ] <- sm$estimate
se[i, ] <- sm$se
}
eval.points <- sm$eval.points
if (!(opt$display %in% "none" | band)) {
plot(xlim, c(0, 1.1 * max(as.vector(estimate))), xlab = opt$xlab,
ylab = opt$ylab, type = "n")
#for (i in 1:nlev) lines(eval.points, estimate[i, ], lty =
opt$lty[i],
# col = opt$col[i])
for (i in 1:nlev) lines(eval.points, estimate[i, ], lty =
opt$lty[i], ## lwd hacked in
col = opt$col[i], lwd = lwd[i])
}
est <- NULL
p <- NULL
if (model == "equal" & test) {
if (nlev == 2) {
ts <- sum((estimate[1, ] - estimate[2, ])^2)
}
else {
sm.mean <- sm.density(y, h = h, xlim = opt$xlim,
ngrid = opt$ngrid, display = "none")$estimate
ts <- 0
for (i in 1:nlev) ts <- ts + ni[i] * sum((estimate[i,
] - sm.mean)^2)
}
p <- 0
est.star <- matrix(0, ncol = opt$ngrid, nrow = nlev)
for (iboot in 1:nboot) {
ind <- (1:length(y))
for (i in 1:nlev) {
indi <- sample((1:length(ind)), ni[i])
est.star[i, ] <- sm.density(y[ind[indi]], h = h,
ngrid = opt$ngrid, xlim = opt$xlim, display =
"none")$estimate
ind <- ind[-indi]
}
if (nlev == 2) {
ts.star <- sum((est.star[1, ] - est.star[2, ])^2)
}
else {
sm.mean <- sm.density(y, h = h, xlim = opt$xlim,
ngrid = opt$ngrid, display = "none")$estimate
ts.star <- 0
for (i in 1:nlev) {
ts.star <- ts.star + ni[i] * sum((est.star[i,
] - sm.mean)^2)
}
}
if (ts.star > ts)
p <- p + 1
if (opt$verbose > 1) {
cat(iboot)
cat(" ")
}
}
p <- p/nboot
cat("\nTest of equal densities: p-value = ", round(p,
3), "\n")
est <- list(p = p, h = h)
}
if (model == "equal" & band) {
av <- (sqrt(estimate[1, ]) + sqrt(estimate[2, ]))/2
se <- sqrt(se[1, ]^2 + se[2, ]^2)
upper <- (av + se)^2
lower <- pmax(av - se, 0)^2
plot(xlim, c(0, 1.1 * max(as.vector(estimate), upper)),
xlab = xlab, ylab = ylab, type = "n", asp=asp, ...) ## ...
and asp added; was opt$xlab and opt$ylab
polygon(c(eval.points, rev(eval.points)), c(upper, rev(lower)),
col = bandcol, border = 0)
## was col = "cyan"
if (is.null(usePolyg)) {
lines(eval.points, estimate[1, ], lty = opt$lty[1], col =
opt$col[1], lwd = lwd[1])
lines(eval.points, estimate[2, ], lty = opt$lty[2], col =
opt$col[2], lwd = lwd[2])
}
else {
polygon(eval.points, estimate[1, ], lty = opt$lty[1], col =
opt$col[1], lwd = lwd[1])
polygon(eval.points, estimate[2, ], lty = opt$lty[2], col =
opt$col[2], lwd = lwd[2])
}
est <- list(p = p, upper = upper, lower = lower, h = h)
}
invisible(est)
}
Viviana Ruiz wrote:
>
> Hi,
>
> I am doing kernel density plots, and am trying to make the lines thicker.
> I
> comparing three groups, in sm.density.compare. I tried changing lwd to
> make
> the line sthicker right on the density compare call, but was not able to
> do
> it. There is not an option in sm.options to specify line thickness, as
> well
> as cex.ylab or cex.xlab- I tried it and it does not change the thickness
> of
> the lines. Does anyone know how to do this?
>
> Thanks!
>
> VRG
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
--
View this message in context: http://www.nabble.com/sm.options-tp22367057p22367994.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list