1## For examples skipped in testing because they need recommended packages. 2 3## This is skipped entirely on a Unix-alike if recommended packages are, 4## so for Windows 5if(!require("MASS")) q() 6 7pdf("reg-examples-3.pdf", encoding = "ISOLatin1.enc") 8 9## From datasets 10if(require("survival")) { 11 model3 <- clogit(case ~ spontaneous+induced+strata(stratum), data = infert) 12 print(summary(model3)) 13 detach("package:survival", unload = TRUE) # survival (conflicts) 14} 15 16 17## From grDevices 18x1 <- matrix(rnorm(1e3), ncol = 2) 19x2 <- matrix(rnorm(1e3, mean = 3, sd = 1.5), ncol = 2) 20x <- rbind(x1, x2) 21 22dcols <- densCols(x) 23graphics::plot(x, col = dcols, pch = 20, main = "n = 1000") 24 25 26## From graphics: 27## A largish data set 28set.seed(123) 29n <- 10000 30x1 <- matrix(rnorm(n), ncol = 2) 31x2 <- matrix(rnorm(n, mean = 3, sd = 1.5), ncol = 2) 32x <- rbind(x1, x2) 33 34oldpar <- par(mfrow = c(2, 2)) 35smoothScatter(x, nrpoints = 0) 36smoothScatter(x) 37 38## a different color scheme: 39Lab.palette <- colorRampPalette(c("blue", "orange", "red"), space = "Lab") 40smoothScatter(x, colramp = Lab.palette) 41 42## somewhat similar, using identical smoothing computations, 43## but considerably *less* efficient for really large data: 44plot(x, col = densCols(x), pch = 20) 45 46## use with pairs: 47par(mfrow = c(1, 1)) 48y <- matrix(rnorm(40000), ncol = 4) + 3*rnorm(10000) 49y[, c(2,4)] <- -y[, c(2,4)] 50pairs(y, panel = function(...) smoothScatter(..., nrpoints = 0, add = TRUE)) 51 52par(oldpar) 53 54 55## From stats 56# alias.Rd 57op <- options(contrasts = c("contr.helmert", "contr.poly")) 58npk.aov <- aov(yield ~ block + N*P*K, npk) 59alias(npk.aov) 60options(op) # reset 61 62# as.hclust.Rd 63if(require("cluster", quietly = TRUE)) {# is a recommended package 64 set.seed(123) 65 x <- matrix(rnorm(30), ncol = 3) 66 hc <- hclust(dist(x), method = "complete") 67 ag <- agnes(x, method = "complete") 68 hcag <- as.hclust(ag) 69 ## The dendrograms order slightly differently: 70 op <- par(mfrow = c(1,2)) 71 plot(hc) ; mtext("hclust", side = 1) 72 plot(hcag); mtext("agnes", side = 1) 73 detach("package:cluster") 74} 75 76# confint.Rd 77counts <- c(18,17,15,20,10,20,25,13,12) 78outcome <- gl(3, 1, 9); treatment <- gl(3, 3) 79glm.D93 <- glm(counts ~ outcome + treatment, family = poisson()) 80confint(glm.D93) 81confint.default(glm.D93) # based on asymptotic normality} 82 83# contrasts.Rd 84utils::example(factor) 85fff <- ff[, drop = TRUE] # reduce to 5 levels. 86contrasts(fff) <- contr.sum(5)[, 1:2]; contrasts(fff) 87 88## using sparse contrasts: % useful, once model.matrix() works with these : 89ffs <- fff 90contrasts(ffs) <- contr.sum(5, sparse = TRUE)[, 1:2]; contrasts(ffs) 91stopifnot(all.equal(ffs, fff)) 92contrasts(ffs) <- contr.sum(5, sparse = TRUE); contrasts(ffs) 93 94# glm.Rd 95utils::data(anorexia, package = "MASS") 96 97anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt), 98 family = gaussian, data = anorexia) 99summary(anorex.1) 100 101# logLik.Rd 102utils::data(Orthodont, package = "nlme") 103fm1 <- lm(distance ~ Sex * age, Orthodont) 104logLik(fm1) 105logLik(fm1, REML = TRUE) 106 107# nls.Rd 108od <- options(digits=5) 109## The muscle dataset in MASS is from an experiment on muscle 110## contraction on 21 animals. The observed variables are Strip 111## (identifier of muscle), Conc (Cacl concentration) and Length 112## (resulting length of muscle section). 113utils::data(muscle, package = "MASS") 114 115## The non linear model considered is 116## Length = alpha + beta*exp(-Conc/theta) + error 117## where theta is constant but alpha and beta may vary with Strip. 118 119with(muscle, table(Strip)) # 2, 3 or 4 obs per strip 120 121## We first use the plinear algorithm to fit an overall model, 122## ignoring that alpha and beta might vary with Strip. 123 124musc.1 <- nls(Length ~ cbind(1, exp(-Conc/th)), muscle, 125 start = list(th = 1), algorithm = "plinear") 126## IGNORE_RDIFF_BEGIN 127summary(musc.1) 128## IGNORE_RDIFF_END 129 130## Then we use nls' indexing feature for parameters in non-linear 131## models to use the conventional algorithm to fit a model in which 132## alpha and beta vary with Strip. The starting values are provided 133## by the previously fitted model. 134## Note that with indexed parameters, the starting values must be 135## given in a list (with names): 136b <- coef(musc.1) 137musc.2 <- nls(Length ~ a[Strip] + b[Strip]*exp(-Conc/th), muscle, 138 start = list(a = rep(b[2], 21), b = rep(b[3], 21), th = b[1])) 139## IGNORE_RDIFF_BEGIN 140summary(musc.2) 141## IGNORE_RDIFF_END 142options(od) 143 144# princomp.Rd 145## Robust: 146(pc.rob <- princomp(stackloss, covmat = MASS::cov.rob(stackloss))) 147 148# termplot.R 149library(MASS) 150hills.lm <- lm(log(time) ~ log(climb)+log(dist), data = hills) 151termplot(hills.lm, partial.resid = TRUE, smooth = panel.smooth, 152 terms = "log(dist)", main = "Original") 153termplot(hills.lm, transform.x = TRUE, 154 partial.resid = TRUE, smooth = panel.smooth, 155 terms = "log(dist)", main = "Transformed") 156 157# xtabs.Rd 158if(require("Matrix")) { 159 ## similar to "nlme"s 'ergoStool' : 160 d.ergo <- data.frame(Type = paste0("T", rep(1:4, 9*4)), 161 Subj = gl(9, 4, 36*4)) 162 print(xtabs(~ Type + Subj, data = d.ergo)) # 4 replicates each 163 set.seed(15) # a subset of cases: 164 print(xtabs(~ Type + Subj, data = d.ergo[sample(36, 10), ], sparse = TRUE)) 165 166 ## Hypothetical two-level setup: 167 inner <- factor(sample(letters[1:25], 100, replace = TRUE)) 168 inout <- factor(sample(LETTERS[1:5], 25, replace = TRUE)) 169 fr <- data.frame(inner = inner, outer = inout[as.integer(inner)]) 170 print(xtabs(~ inner + outer, fr, sparse = TRUE)) 171} 172 173## From utils 174example(packageDescription) 175 176 177## From splines 178library(splines) 179Matrix::drop0(zapsmall(6*splineDesign(knots = 1:40, x = 4:37, sparse = TRUE))) 180