R version 3.1.1 RC (2014-07-04 r66081) -- "Sock it to Me"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> ## 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()
Loading required package: MASS
>
> 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") # survival (conflicts)
+ }
Loading required package: survival
Loading required package: splines
Call:
coxph(formula = Surv(rep(1, 248L), case) ~ spontaneous + induced +
strata(stratum), data = infert, method = "exact")
n= 248, number of events= 83
coef exp(coef) se(coef) z Pr(>|z|)
spontaneous 1.9859 7.2854 0.3524 5.635 1.75e-08 ***
induced 1.4090 4.0919 0.3607 3.906 9.38e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
spontaneous 7.285 0.1373 3.651 14.536
induced 4.092 0.2444 2.018 8.298
Rsquare= 0.193 (max possible= 0.519 )
Likelihood ratio test= 53.15 on 2 df, p=2.869e-12
Wald test = 31.84 on 2 df, p=1.221e-07
Score (logrank) test = 48.44 on 2 df, p=3.032e-11
>
>
> ## 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)
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
> 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)
Model :
yield ~ block + N * P * K
Complete :
(Intercept) block1 block2 block3 block4 block5 N1 P1 K1 N1:P1
N1:P1:K1 0 1 1/3 1/6 -3/10 -1/5 0 0 0 0
N1:K1 P1:K1
N1:P1:K1 0 0
> 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)
Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 2.6958215 3.36655581
outcome2 -0.8577018 -0.06255840
outcome3 -0.6753696 0.08244089
treatment2 -0.3932548 0.39325483
treatment3 -0.3932548 0.39325483
> confint.default(glm.D93) # based on asymptotic normality}
2.5 % 97.5 %
(Intercept) 2.7095672 3.37947764
outcome2 -0.8505027 -0.05800787
outcome3 -0.6707552 0.08478093
treatment2 -0.3919928 0.39199279
treatment3 -0.3919928 0.39199279
>
> # contrasts.Rd
> utils::example(factor)
factor> (ff <- factor(substring("statistics", 1:10, 1:10), levels = letters))
[1] s t a t i s t i c s
Levels: a b c d e f g h i j k l m n o p q r s t u v w x y z
factor> as.integer(ff) # the internal codes
[1] 19 20 1 20 9 19 20 9 3 19
factor> (f. <- factor(ff)) # drops the levels that do not occur
[1] s t a t i s t i c s
Levels: a c i s t
factor> ff[, drop = TRUE] # the same, more transparently
[1] s t a t i s t i c s
Levels: a c i s t
factor> factor(letters[1:20], labels = "letter")
[1] letter1 letter2 letter3 letter4 letter5 letter6 letter7 letter8
[9] letter9 letter10 letter11 letter12 letter13 letter14 letter15 letter16
[17] letter17 letter18 letter19 letter20
20 Levels: letter1 letter2 letter3 letter4 letter5 letter6 letter7 ... letter20
factor> class(ordered(4:1)) # "ordered", inheriting from "factor"
[1] "ordered" "factor"
factor> z <- factor(LETTERS[3:1], ordered = TRUE)
factor> ## and "relational" methods work:
factor> stopifnot(sort(z)[c(1,3)] == range(z), min(z) < max(z))
factor> ## Don't show:
factor> of <- ordered(ff)
factor> stopifnot(identical(range(of, rev(of)), of[3:2]),
factor+ identical(max(of), of[2]))
factor> ## End Don't show
factor>
factor> ## suppose you want "NA" as a level, and to allow missing values.
factor> (x <- factor(c(1, 2, NA), exclude = NULL))
[1] 1 2 <NA>
Levels: 1 2 <NA>
factor> is.na(x)[2] <- TRUE
factor> x # [1] 1 <NA> <NA>
[1] 1 <NA> <NA>
Levels: 1 2 <NA>
factor> is.na(x)
[1] FALSE TRUE FALSE
factor> # [1] FALSE TRUE FALSE
factor>
factor> ## Using addNA()
factor> Month <- airquality$Month
factor> table(addNA(Month))
5 6 7 8 9 <NA>
31 30 31 31 30 0
factor> table(addNA(Month, ifany = TRUE))
5 6 7 8 9
31 30 31 31 30
> fff <- ff[, drop = TRUE] # reduce to 5 levels.
> contrasts(fff) <- contr.sum(5)[, 1:2]; contrasts(fff)
[,1] [,2] [,3] [,4]
a 1 0 -0.2471257 0.2688164
c 0 1 -0.2471257 0.2688164
i 0 0 -0.1498721 -0.8817814
s 0 0 0.8912491 0.0753323
t -1 -1 -0.2471257 0.2688164
>
> ## using sparse contrasts: % useful, once model.matrix() works with these :
> ffs <- fff
> contrasts(ffs) <- contr.sum(5, sparse = TRUE)[, 1:2]; contrasts(ffs)
[,1] [,2] [,3] [,4]
a 1 0 -0.2471257 0.2688164
c 0 1 -0.2471257 0.2688164
i 0 0 -0.1498721 -0.8817814
s 0 0 0.8912491 0.0753323
t -1 -1 -0.2471257 0.2688164
> stopifnot(all.equal(ffs, fff))
> contrasts(ffs) <- contr.sum(5, sparse = TRUE); contrasts(ffs)
5 x 4 sparse Matrix of class "dgCMatrix"
a 1 . . .
c . 1 . .
i . . 1 .
s . . . 1
t -1 -1 -1 -1
>
> # glm.Rd
> utils::data(anorexia, package = "MASS")
>
> anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
+ family = gaussian, data = anorexia)
> summary(anorex.1)
Call:
glm(formula = Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian,
data = anorexia)
Deviance Residuals:
Min 1Q Median 3Q Max
-14.1083 -4.2773 -0.5484 5.4838 15.2922
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.7711 13.3910 3.717 0.000410 ***
Prewt -0.5655 0.1612 -3.509 0.000803 ***
TreatCont -4.0971 1.8935 -2.164 0.033999 *
TreatFT 4.5631 2.1333 2.139 0.036035 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 48.69504)
Null deviance: 4525.4 on 71 degrees of freedom
Residual deviance: 3311.3 on 68 degrees of freedom
AIC: 489.97
Number of Fisher Scoring iterations: 2
>
> # 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")
+ summary(musc.1)
+
+ ## 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]))
+ summary(musc.2)
+ options(od)
+
+ # princomp.Rd
+ ## Robust:
+ (pc.rob <- princomp(stackloss, covmat = MASS::cov.rob(stackloss)))
+
+ # 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 tools
+
+ ## This may not be installed
+ gridEx <- system.file("doc", "grid.Rnw", package = "grid")
+ vignetteDepends(gridEx)
+
+
+ ## From utils
+ example(packageDescription)
+
+
+ ## From splines
+ library(splines)
+ Matrix::drop0(zapsmall(6*splineDesign(knots = 1:40, x = 4:37, sparse = TRUE)))
+
+