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