## For examples skipped in testing because they need recommended packages. ## This is skipped entirely on a Unix-alike if recommended packages are, ## so for Windows if(!require("MASS")) q() pdf("reg-examples-3.pdf", encoding = "ISOLatin1.enc") ## From datasets if(require("survival")) { model3 <- clogit(case ~ spontaneous+induced+strata(stratum), data = infert) print(summary(model3)) detach("package:survival", unload = TRUE) # survival (conflicts) } ## From grDevices x1 <- matrix(rnorm(1e3), ncol = 2) x2 <- matrix(rnorm(1e3, mean = 3, sd = 1.5), ncol = 2) x <- rbind(x1, x2) dcols <- densCols(x) graphics::plot(x, col = dcols, pch = 20, main = "n = 1000") ## From graphics: ## A largish data set set.seed(123) n <- 10000 x1 <- matrix(rnorm(n), ncol = 2) x2 <- matrix(rnorm(n, mean = 3, sd = 1.5), ncol = 2) x <- rbind(x1, x2) oldpar <- par(mfrow = c(2, 2)) smoothScatter(x, nrpoints = 0) smoothScatter(x) ## a different color scheme: Lab.palette <- colorRampPalette(c("blue", "orange", "red"), space = "Lab") smoothScatter(x, colramp = Lab.palette) ## somewhat similar, using identical smoothing computations, ## but considerably *less* efficient for really large data: plot(x, col = densCols(x), pch = 20) ## use with pairs: par(mfrow = c(1, 1)) y <- matrix(rnorm(40000), ncol = 4) + 3*rnorm(10000) y[, c(2,4)] <- -y[, c(2,4)] pairs(y, panel = function(...) smoothScatter(..., nrpoints = 0, add = TRUE)) par(oldpar) ## From stats # alias.Rd op <- options(contrasts = c("contr.helmert", "contr.poly")) npk.aov <- aov(yield ~ block + N*P*K, npk) alias(npk.aov) options(op) # reset # as.hclust.Rd if(require("cluster", quietly = TRUE)) {# is a recommended package set.seed(123) x <- matrix(rnorm(30), ncol = 3) hc <- hclust(dist(x), method = "complete") ag <- agnes(x, method = "complete") hcag <- as.hclust(ag) ## The dendrograms order slightly differently: op <- par(mfrow = c(1,2)) plot(hc) ; mtext("hclust", side = 1) plot(hcag); mtext("agnes", side = 1) detach("package:cluster") } # confint.Rd counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3, 1, 9); treatment <- gl(3, 3) glm.D93 <- glm(counts ~ outcome + treatment, family = poisson()) confint(glm.D93) confint.default(glm.D93) # based on asymptotic normality} # contrasts.Rd utils::example(factor) fff <- ff[, drop = TRUE] # reduce to 5 levels. contrasts(fff) <- contr.sum(5)[, 1:2]; contrasts(fff) ## using sparse contrasts: % useful, once model.matrix() works with these : ffs <- fff contrasts(ffs) <- contr.sum(5, sparse = TRUE)[, 1:2]; contrasts(ffs) stopifnot(all.equal(ffs, fff)) contrasts(ffs) <- contr.sum(5, sparse = TRUE); contrasts(ffs) # glm.Rd utils::data(anorexia, package = "MASS") anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian, data = anorexia) summary(anorex.1) # logLik.Rd utils::data(Orthodont, package = "nlme") fm1 <- lm(distance ~ Sex * age, Orthodont) logLik(fm1) logLik(fm1, REML = TRUE) # nls.Rd od <- options(digits=5) ## The muscle dataset in MASS is from an experiment on muscle ## contraction on 21 animals. The observed variables are Strip ## (identifier of muscle), Conc (Cacl concentration) and Length ## (resulting length of muscle section). utils::data(muscle, package = "MASS") ## The non linear model considered is ## Length = alpha + beta*exp(-Conc/theta) + error ## where theta is constant but alpha and beta may vary with Strip. with(muscle, table(Strip)) # 2, 3 or 4 obs per strip ## We first use the plinear algorithm to fit an overall model, ## ignoring that alpha and beta might vary with Strip. musc.1 <- nls(Length ~ cbind(1, exp(-Conc/th)), muscle, start = list(th = 1), algorithm = "plinear") ## IGNORE_RDIFF_BEGIN summary(musc.1) ## IGNORE_RDIFF_END ## Then we use nls' indexing feature for parameters in non-linear ## models to use the conventional algorithm to fit a model in which ## alpha and beta vary with Strip. The starting values are provided ## by the previously fitted model. ## Note that with indexed parameters, the starting values must be ## given in a list (with names): b <- coef(musc.1) musc.2 <- nls(Length ~ a[Strip] + b[Strip]*exp(-Conc/th), muscle, start = list(a = rep(b[2], 21), b = rep(b[3], 21), th = b[1])) ## IGNORE_RDIFF_BEGIN summary(musc.2) ## IGNORE_RDIFF_END options(od) # princomp.Rd ## Robust: (pc.rob <- princomp(stackloss, covmat = MASS::cov.rob(stackloss))) # termplot.R library(MASS) hills.lm <- lm(log(time) ~ log(climb)+log(dist), data = hills) termplot(hills.lm, partial.resid = TRUE, smooth = panel.smooth, terms = "log(dist)", main = "Original") termplot(hills.lm, transform.x = TRUE, partial.resid = TRUE, smooth = panel.smooth, terms = "log(dist)", main = "Transformed") # xtabs.Rd if(require("Matrix")) { ## similar to "nlme"s 'ergoStool' : d.ergo <- data.frame(Type = paste0("T", rep(1:4, 9*4)), Subj = gl(9, 4, 36*4)) print(xtabs(~ Type + Subj, data = d.ergo)) # 4 replicates each set.seed(15) # a subset of cases: print(xtabs(~ Type + Subj, data = d.ergo[sample(36, 10), ], sparse = TRUE)) ## Hypothetical two-level setup: inner <- factor(sample(letters[1:25], 100, replace = TRUE)) inout <- factor(sample(LETTERS[1:5], 25, replace = TRUE)) fr <- data.frame(inner = inner, outer = inout[as.integer(inner)]) print(xtabs(~ inner + outer, fr, sparse = TRUE)) } ## From splines library(splines) Matrix::drop0(zapsmall(6*splineDesign(knots = 1:40, x = 4:37, sparse = TRUE))) ## From tools library(tools) example(standard_package_names, run.donttest = TRUE)