R Under development (unstable) (2026-03-06 r89553) -- "Unsuffered Consequences" Copyright (C) 2026 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu 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", unload = TRUE) # survival (conflicts) + } Loading required package: survival 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 Concordance= 0.776 (se = 0.044 ) Likelihood ratio test= 53.15 on 2 df, p=3e-12 Wald test = 31.84 on 2 df, p=1e-07 Score (logrank) test = 48.44 on 2 df, p=3e-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) > 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") + } > > # 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 Levels: 1 2 factor> is.na(x)[2] <- TRUE factor> x # [1] 1 [1] 1 Levels: 1 2 factor> is.na(x) [1] FALSE TRUE FALSE factor> # [1] FALSE TRUE FALSE factor> factor> ## More rational, since R 3.4.0 : factor> factor(c(1:2, NA), exclude = "" ) # keeps , as [1] 1 2 Levels: 1 2 factor> factor(c(1:2, NA), exclude = NULL) # always did [1] 1 2 Levels: 1 2 factor> ## exclude = factor> z # ordered levels 'A < B < C' [1] C B A Levels: A < B < C factor> factor(z, exclude = "C") # does exclude [1] B A Levels: A < B factor> factor(z, exclude = "B") # ditto [1] C A Levels: A < C factor> ## Now, labels maybe duplicated: factor> ## factor() with duplicated labels allowing to "merge levels" factor> x <- c("Man", "Male", "Man", "Lady", "Female") factor> ## Map from 4 different values to only two levels: factor> (xf <- factor(x, levels = c("Male", "Man" , "Lady", "Female"), factor+ labels = c("Male", "Male", "Female", "Female"))) [1] Male Male Male Female Female Levels: Male Female factor> #> [1] Male Male Male Female Female factor> #> Levels: Male Female factor> factor> ## Using addNA() factor> Month <- airquality$Month factor> table(addNA(Month)) 5 6 7 8 9 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 > if(requireNamespace("Matrix", lib.loc=.Library)) withAutoprint({ + 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) + }) > 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 > > # free1way() et al. > example(free1way, run.donttest = TRUE) fre1wy> ## No test: fre1wy> fre1wy> ## Kruskal-Wallis test fre1wy> kruskal.test(Ozone ~ Month, data = airquality) Kruskal-Wallis rank sum test data: Ozone by Month Kruskal-Wallis chi-squared = 29.267, df = 4, p-value = 6.901e-06 fre1wy> kt <- free1way(Ozone ~ Month, data = airquality) fre1wy> print(kt) 5-sample Kruskal-Wallis test against proportional odds alternatives data: Ozone by Month (5, 6, 7, 8, 9) Perm chi-squared = 29.267, df = 4, p-value = 6.901e-06 alternative hypothesis: true log-odds ratio is not equal to 0 fre1wy> # log-odds ratios for comparison with control fre1wy> coef(kt) Month6 Month7 Month8 Month9 0.8123639 2.5281594 2.3825979 0.7513235 fre1wy> # Wald inference fre1wy> summary(kt) Call: free1way.formula(formula = Ozone ~ Month, data = airquality) Coefficients: log-odds ratio Std. Error z value P(>|z|) Month6 0.81236 0.64415 1.26113 0.2073 Month7 2.52816 0.53153 4.75639 0.0000 Month8 2.38260 0.53370 4.46432 0.0000 Month9 0.75132 0.47177 1.59256 0.1113 fre1wy> confint(kt, test = "Wald") 2.5 % 97.5 % Month6 -0.4501567 2.074884 Month7 1.4863809 3.569938 Month8 1.3365705 3.428625 Month9 -0.1733322 1.675979 fre1wy> ## Friedman test fre1wy> example(friedman.test, echo = FALSE) fre1wy> me <- colnames(RoundingTimes) fre1wy> d <- expand.grid(me = factor(me, labels = me, levels = me), fre1wy+ id = factor(seq_len(nrow(RoundingTimes)))) fre1wy> d$time <- c(t(RoundingTimes)) fre1wy> # global p-value identical fre1wy> friedman.test(RoundingTimes) Friedman rank sum test data: RoundingTimes Friedman chi-squared = 11.143, df = 2, p-value = 0.003805 fre1wy> ft <- free1way(time ~ me | id, data = d) fre1wy> print(ft) Stratified 3-sample Kruskal-Wallis test against proportional odds alternatives data: time by me (Round Out, Narrow Angle, Wide Angle) stratified by id Perm chi-squared = 11.143, df = 2, p-value = 0.003805 alternative hypothesis: true log-odds ratio is not equal to 0 fre1wy> coef(ft) meNarrow Angle meWide Angle -1.265280 -3.289198 fre1wy> # Wald inference fre1wy> summary(ft) Call: free1way.formula(formula = time ~ me | id, data = d) Coefficients: log-odds ratio Std. Error z value P(>|z|) meNarrow Angle -1.26528 0.68264 -1.85352 0.0638 meWide Angle -3.28920 0.79367 -4.14428 0.0000 fre1wy> confint(ft, test = "Wald") 2.5 % 97.5 % meNarrow Angle -2.603222 0.07266103 meWide Angle -4.844765 -1.73363038 fre1wy> ## McNemar test fre1wy> ## paired binary observations fre1wy> example(mcnemar.test, echo = FALSE) fre1wy> # set-up data frame with survey outcomes for voters fre1wy> s <- gl(2, 1, labels = dimnames(Performance)[[1L]]) fre1wy> survey <- gl(2, 1, labels = c("1st", "2nd")) fre1wy> nvoters <- c(Performance) fre1wy> x <- expand.grid(survey = survey, voter = factor(seq_len(sum(nvoters)))) fre1wy> x$performance <- c(rep(s[c(1, 1)], nvoters[1]), rep(s[c(2, 1)], nvoters[2]), fre1wy+ rep(s[c(1, 2)], nvoters[3]), rep(s[c(2, 2)], nvoters[4])) fre1wy> # note that only those voters changing their minds are relevant fre1wy> mcn <- free1way(xtabs(~ performance + survey + voter, data = x)) fre1wy> # same result as mcnemar.test w/o continuity correction fre1wy> print(mcn) Stratified 2-sample Wilcoxon test against proportional odds alternatives data: performance by survey (c("1st", "2nd")) stratified by voter Perm Z = 4.166, p-value = 3.099e-05 alternative hypothesis: true log-odds ratio is not equal to 0 fre1wy> # X^2 statistic fre1wy> summary(mcn, test = "Permutation")$statistic^2 Perm Z 17.35593 fre1wy> mcnemar.test(Performance, correct = FALSE) McNemar's Chi-squared test data: Performance McNemar's chi-squared = 17.356, df = 1, p-value = 3.099e-05 fre1wy> # Wald inference fre1wy> summary(mcn) Call: free1way.table(y = xtabs(~performance + survey + voter, data = x)) Coefficients: log-odds ratio Std. Error z value P(>|z|) survey2nd 1.11258 0.19128 5.81639 0 fre1wy> confint(mcn, test = "Wald") 2.5 % 97.5 % survey2nd 0.7376684 1.487484 fre1wy> ## Mantel-Haenszel test w/o continuity correction, fre1wy> ## Departments are blocks fre1wy> mantelhaen.test(UCBAdmissions, correct = FALSE) Mantel-Haenszel chi-squared test without continuity correction data: UCBAdmissions Mantel-Haenszel X-squared = 1.5246, df = 1, p-value = 0.2169 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 0.7719074 1.0603298 sample estimates: common odds ratio 0.9046968 fre1wy> mh <- free1way(UCBAdmissions) fre1wy> print(mh) Stratified 2-sample Wilcoxon test against proportional odds alternatives data: Admit by Gender (c("Male", "Female")) stratified by Dept Perm Z = -1.2347, p-value = 0.2169 alternative hypothesis: true log-odds ratio is not equal to 0 fre1wy> # common odds-ratio, with score interval fre1wy> exp(coef(mh)) GenderFemale 0.904955 fre1wy> exp(confint(mh, test = "Rao")) 2.5 % 97.5 % GenderFemale 0.7723922 1.06027 fre1wy> # looking at department-specific fre1wy> # confidence intervals for log-odds ratios fre1wy> # it seems Dept A is out of line fre1wy> apply(UCBAdmissions, MARGIN = 3, fre1wy+ FUN = function(x) confint(free1way(as.table(x)))) Dept A B C D E F [1,] -1.5622307 -1.0556769 -0.1570020 -0.3762303 -0.1913598 -0.7802366 [2,] -0.5415513 0.6167797 0.4068713 0.2122505 0.5918809 0.4025005 fre1wy> ## Mantel-Haenszel test treats variables as fre1wy> ## unordered, free1way allows ordered responses fre1wy> example(mantelhaen.test, echo = FALSE) fre1wy> # Does distribution of job satisfaction (ordered) depend on income fre1wy> # in a stratified proportional odds model? fre1wy> # Job Satisfaction is second in array but needs to be first fre1wy> # for free1way to treat it as ordered response fre1wy> ft <- free1way(aperm(Satisfaction, perm = c(2, 1, 3))) fre1wy> summary(ft) Call: free1way.table(y = aperm(Satisfaction, perm = c(2, 1, 3))) Coefficients: log-odds ratio Std. Error z value P(>|z|) Income5000-15000 0.10794 0.55953 0.19290 0.8470 Income15000-25000 1.45998 0.62561 2.33369 0.0196 Income>25000 1.23282 0.65693 1.87663 0.0606 fre1wy> ## End(No test) fre1wy> fre1wy> > example(power.free1way.test, run.donttest = TRUE) pwr.1.> ## No test: pwr.1.> pwr.1.> ## make example reproducible pwr.1.> set.seed(29) pwr.1.> ## sample from proportional odds model with 1:2 allocation pwr.1.> ## based on odds ratio of 3, with sample sizes (15, 30) pwr.1.> x <- rfree1way(n = 15, delta = log(3), alloc_ratio = 2) pwr.1.> # Wilcoxon-Mann-Whitney rank sum test via classical stats interface pwr.1.> wilcox.test(y ~ groups, data = x, exact = FALSE, correct = FALSE)$p.value [1] 0.008679882 pwr.1.> # Identical p-value obtained from a proportional-odds model pwr.1.> summary(free1way(y ~ groups, data = x), test = "Permutation")$p.value Perm Z 0.008679882 pwr.1.> # approximate power for this test pwr.1.> power.free1way.test(n = 15, delta = log(3), alloc_ratio = 2) 2-sample Wilcoxon test against proportional odds alternatives n = 15 Total sample size = 15 (Control) + 30 (B) = 45 power = 0.4879209 sig.level = 0.05 Standard error = 0.5693603 log-odds ratio = 1.098612 NOTE: 'n' is sample size in control group pwr.1.> ## End(No test) pwr.1.> pwr.1.> > example(ppplot, run.donttest = TRUE) ppplot> ## No test: ppplot> ppplot> ## make example reproducible ppplot> set.seed(29) ppplot> ## well-fitting logistic model ppplot> nd <- data.frame(groups = gl(2, 50, labels = paste0("G", 1:2))) ppplot> nd$y <- rlogis(nrow(nd), location = c(0, 2)[nd$groups]) ppplot> with(with(nd, split(y, groups)), ppplot+ ppplot(G1, G2, conf.level = .95, ppplot+ conf.args = list(link = "logit", type = "Wald", col = 2))) ppplot> # with appropriate Wilcoxon test and log-odds ratio ppplot> coef(ft <- free1way(y ~ groups, data = nd)) groupsG2 1.668854 ppplot> # the model-based probability-probability curve ppplot> prb <- 1:99 / 100 ppplot> points(prb, plogis(qlogis(prb) - coef(ft)), pch = 3) ppplot> ## the corresponding model-based receiver operating characteristic (ROC) ppplot> ## curve, see Sewak and Hothorn (2023) ppplot> plot(prb, plogis(qlogis(1 - prb) - coef(ft), lower.tail = FALSE), ppplot+ xlab = "1 - Specificity", ylab = "Sensitivity", type = "l", ppplot+ main = "ROC Curve") ppplot> abline(a = 0, b = 1, col = "lightgrey") ppplot> # with confidence band ppplot> lines(prb, plogis(qlogis(1 - prb) - confint(ft, test = "Rao")[1], ppplot+ lower.tail = FALSE), lty = 3) ppplot> lines(prb, plogis(qlogis(1 - prb) - confint(ft, test = "Rao")[2], ppplot+ lower.tail = FALSE), lty = 3) ppplot> # and corresponding area under the ROC curve (AUC) ppplot> # with score confidence interval ppplot> coef(ft, what = "AUC") groupsG2 0.75467 ppplot> confint(ft, test = "Rao", what = "AUC") 2.5 % 97.5 % groupsG2 0.6490714 0.8381977 ppplot> ## ill-fitting normal model ppplot> nd$y <- rnorm(nrow(nd), mean = c(0, .5)[nd$groups], sd = c(1, 1.5)[nd$groups]) ppplot> with(with(nd, split(y, groups)), ppplot+ ppplot(G1, G2, conf.level = .95, ppplot+ conf.args = list(link = "probit", type = "Wald", col = 2))) ppplot> # inappropriate probit model ppplot> coef(free1way(y ~ groups, data = nd, link = "probit")) groupsG2 0.6659502 ppplot> ## End(No test) ppplot> ppplot> > > # 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) 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) 'log Lik.' -239.1209 (df=5) > logLik(fm1, REML = TRUE) 'log Lik.' -241.7796 (df=5) > > # 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 Strip S01 S02 S03 S04 S05 S06 S07 S08 S09 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 4 4 4 3 3 3 2 2 2 2 3 2 2 2 2 4 4 3 3 3 S21 3 > > ## 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) Formula: Length ~ cbind(1, exp(-Conc/th)) Parameters: Estimate Std. Error t value Pr(>|t|) th 0.608 0.115 5.31 1.9e-06 *** .lin1 28.963 1.230 23.55 < 2e-16 *** .lin2 -34.227 3.793 -9.02 1.4e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.67 on 57 degrees of freedom Number of iterations to convergence: 5 Achieved convergence tolerance: 9.34e-07 > ## 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) Formula: Length ~ a[Strip] + b[Strip] * exp(-Conc/th) Parameters: Estimate Std. Error t value Pr(>|t|) a1 23.454 0.796 29.46 5.0e-16 *** a2 28.302 0.793 35.70 < 2e-16 *** a3 30.801 1.716 17.95 1.7e-12 *** a4 25.921 3.016 8.60 1.4e-07 *** a5 23.201 2.891 8.02 3.5e-07 *** a6 20.120 2.435 8.26 2.3e-07 *** a7 33.595 1.682 19.98 3.0e-13 *** a8 39.053 3.753 10.41 8.6e-09 *** a9 32.137 3.318 9.69 2.5e-08 *** a10 40.005 3.336 11.99 1.0e-09 *** a11 36.190 3.109 11.64 1.6e-09 *** a12 36.911 1.839 20.07 2.8e-13 *** a13 30.635 1.700 18.02 1.6e-12 *** a14 34.312 3.495 9.82 2.0e-08 *** a15 38.395 3.375 11.38 2.3e-09 *** a16 31.226 0.886 35.26 < 2e-16 *** a17 31.230 0.821 38.02 < 2e-16 *** a18 19.998 1.011 19.78 3.6e-13 *** a19 37.095 1.071 34.65 < 2e-16 *** a20 32.594 1.121 29.07 6.2e-16 *** a21 30.376 1.057 28.74 7.5e-16 *** b1 -27.300 6.873 -3.97 0.00099 *** b2 -26.270 6.754 -3.89 0.00118 ** b3 -30.901 2.270 -13.61 1.4e-10 *** b4 -32.238 3.810 -8.46 1.7e-07 *** b5 -29.941 3.773 -7.94 4.1e-07 *** b6 -20.622 3.647 -5.65 2.9e-05 *** b7 -19.625 8.085 -2.43 0.02661 * b8 -45.780 4.113 -11.13 3.2e-09 *** b9 -31.345 6.352 -4.93 0.00013 *** b10 -38.599 3.955 -9.76 2.2e-08 *** b11 -33.921 3.839 -8.84 9.2e-08 *** b12 -38.268 8.992 -4.26 0.00053 *** b13 -22.568 8.194 -2.75 0.01355 * b14 -36.167 6.358 -5.69 2.7e-05 *** b15 -32.952 6.354 -5.19 7.4e-05 *** b16 -47.207 9.540 -4.95 0.00012 *** b17 -33.875 7.688 -4.41 0.00039 *** b18 -15.896 6.222 -2.55 0.02051 * b19 -28.969 7.235 -4.00 0.00092 *** b20 -36.917 8.033 -4.60 0.00026 *** b21 -26.508 7.012 -3.78 0.00149 ** th 0.797 0.127 6.30 8.0e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.11 on 17 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 2.13e-06 > ## IGNORE_RDIFF_END > options(od) > > # princomp.Rd > ## Robust: > set.seed(1) > (pc.rob <- princomp(stackloss, covmat = MASS::cov.rob(stackloss))) Call: princomp(x = stackloss, covmat = MASS::cov.rob(stackloss)) Standard deviations: Comp.1 Comp.2 Comp.3 Comp.4 7.8322873 4.0077676 1.9114016 0.7624211 4 variables and 21 observations. > > # xtabs.Rd > if(require("Matrix", .Library)) { + ## 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)) + } Loading required package: Matrix Subj Type 1 2 3 4 5 6 7 8 9 T1 4 4 4 4 4 4 4 4 4 T2 4 4 4 4 4 4 4 4 4 T3 4 4 4 4 4 4 4 4 4 T4 4 4 4 4 4 4 4 4 4 4 x 9 sparse Matrix of class "dgCMatrix" Subj Type 1 2 3 4 5 6 7 8 9 T1 . 1 . . . 1 . . . T2 . . . . . . . 1 1 T3 . . 1 1 1 . 1 1 . T4 . 1 . . . . . . . 25 x 4 sparse Matrix of class "dgCMatrix" outer inner A B C E a . . . 3 b . . . 6 c . . 5 . d . 1 . . e . . 3 . f . . 5 . g . . 3 . h 6 . . . i 8 . . . j . 1 . . k . . . 6 l . 2 . . m . . . 7 n . . 3 . o . 3 . . p . 3 . . q . . . 7 r . . . 2 s . 6 . . t . . . 5 u 4 . . . v . . 3 . w . 3 . . x . . 3 . y 2 . . . > > > ## From splines > library(splines) > Matrix::drop0(zapsmall(6*splineDesign(knots = 1:40, x = 4:37, sparse = TRUE))) 34 x 36 sparse Matrix of class "dgCMatrix" [1,] 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [2,] . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [3,] . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [4,] . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [5,] . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [6,] . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . [7,] . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . [8,] . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . . [9,] . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . . [10,] . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . . [11,] . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . . [12,] . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . . [13,] . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . . [14,] . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . . [15,] . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . . [16,] . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . . [17,] . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . . [18,] . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . . [19,] . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . . [20,] . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . . [21,] . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . . [22,] . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . . [23,] . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . . [24,] . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . . [25,] . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . . [26,] . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . . [27,] . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . . [28,] . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . . [29,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . . [30,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . . [31,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . . [32,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . . [33,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 . [34,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 1 > > > ## From tools > library(tools) > example(standard_package_names, run.donttest = TRUE) stnd__> str(stPkgs <- standard_package_names()) List of 2 $ base : chr [1:14] "base" "tools" "utils" "grDevices" ... $ recommended: chr [1:15] "MASS" "lattice" "Matrix" "nlme" ... stnd__> ## No test: stnd__> ## consistency of packageDescription and standard_package_names : stnd__> (pNms <- unlist(stPkgs, FALSE)) base1 base2 base3 base4 base5 "base" "tools" "utils" "grDevices" "graphics" base6 base7 base8 base9 base10 "stats" "datasets" "methods" "grid" "splines" base11 base12 base13 base14 recommended1 "stats4" "tcltk" "compiler" "parallel" "MASS" recommended2 recommended3 recommended4 recommended5 recommended6 "lattice" "Matrix" "nlme" "survival" "boot" recommended7 recommended8 recommended9 recommended10 recommended11 "cluster" "codetools" "foreign" "KernSmooth" "rpart" recommended12 recommended13 recommended14 recommended15 "class" "nnet" "spatial" "mgcv" stnd__> (prio <- sapply(as.vector(pNms), packageDescription, fields = "Priority")) base tools utils grDevices graphics "base" "base" "base" "base" "base" stats datasets methods grid splines "base" "base" "base" "base" "base" stats4 tcltk compiler parallel MASS "base" "base" "base" "base" "recommended" lattice Matrix nlme survival boot "recommended" "recommended" "recommended" "recommended" "recommended" cluster codetools foreign KernSmooth rpart "recommended" "recommended" "recommended" "recommended" "recommended" class nnet spatial mgcv "recommended" "recommended" "recommended" "recommended" stnd__> stopifnot(identical(unname(prio), stnd__+ sub("[0-9]+$", '', names(pNms)))) stnd__> ## End(No test) stnd__> stnd__> stnd__> >