Re: [R] sm.options

From: Mark Difford <mark_difford_at_yahoo.co.uk>
Date: Fri, 06 Mar 2009 00:22:11 -0800 (PST)

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.

______________________________________________
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.
Received on Fri 06 Mar 2009 - 07:24:53 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 06 Mar 2009 - 07:30:23 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive