The R Project SVN
R
(root)/branches/R-4-6-branch/tests/reg-examples3.Rout.save – Rev 89719
Rev 89556 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2026-03-06 r89553) -- "Unsuffered Consequences"Copyright (C) 2026 The R Foundation for Statistical ComputingPlatform: x86_64-pc-linux-gnuR 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: survivalCall:coxph(formula = Surv(rep(1, 248L), case) ~ spontaneous + induced +strata(stratum), data = infert, method = "exact")n= 248, number of events= 83coef 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 ' ' 1exp(coef) exp(-coef) lower .95 upper .95spontaneous 7.285 0.1373 3.651 14.536induced 4.092 0.2444 2.018 8.298Concordance= 0.776 (se = 0.044 )Likelihood ratio test= 53.15 on 2 df, p=3e-12Wald test = 31.84 on 2 df, p=1e-07Score (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 * KComplete :(Intercept) block1 block2 block3 block4 block5 N1 P1 K1 N1:P1N1:P1:K1 0 1 1/3 1/6 -3/10 -1/5 0 0 0 0N1:K1 P1:K1N1: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 sLevels: 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 zfactor> as.integer(ff) # the internal codes[1] 19 20 1 20 9 19 20 9 3 19factor> (f. <- factor(ff)) # drops the levels that do not occur[1] s t a t i s t i c sLevels: a c i s tfactor> ff[, drop = TRUE] # the same, more transparently[1] s t a t i s t i c sLevels: a c i s tfactor> 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 letter2020 Levels: letter1 letter2 letter3 letter4 letter5 letter6 letter7 ... letter20factor> 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] <- TRUEfactor> x # [1] 1 <NA> <NA>[1] 1 <NA> <NA>Levels: 1 2 <NA>factor> is.na(x)[1] FALSE TRUE FALSEfactor> # [1] FALSE TRUE FALSEfactor>factor> ## More rational, since R 3.4.0 :factor> factor(c(1:2, NA), exclude = "" ) # keeps <NA> , as[1] 1 2 <NA>Levels: 1 2 <NA>factor> factor(c(1:2, NA), exclude = NULL) # always did[1] 1 2 <NA>Levels: 1 2 <NA>factor> ## exclude = <character>factor> z # ordered levels 'A < B < C'[1] C B ALevels: A < B < Cfactor> factor(z, exclude = "C") # does exclude[1] <NA> B ALevels: A < Bfactor> factor(z, exclude = "B") # ditto[1] C <NA> ALevels: A < Cfactor> ## 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 FemaleLevels: Male Femalefactor> #> [1] Male Male Male Female Femalefactor> #> Levels: Male Femalefactor>factor> ## Using addNA()factor> Month <- airquality$Monthfactor> table(addNA(Month))5 6 7 8 9 <NA>31 30 31 31 30 0factor> table(addNA(Month, ifany = TRUE))5 6 7 8 931 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.2688164c 0 1 -0.2471257 0.2688164i 0 0 -0.1498721 -0.8817814s 0 0 0.8912491 0.0753323t -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.2688164c 0 1 -0.2471257 0.2688164i 0 0 -0.1498721 -0.8817814s 0 0 0.8912491 0.0753323t -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 . . . 1t -1 -1 -1 -1>> # free1way() et al.> example(free1way, run.donttest = TRUE)fre1wy> ## No test:fre1wy>fre1wy> ## Kruskal-Wallis testfre1wy> kruskal.test(Ozone ~ Month, data = airquality)Kruskal-Wallis rank sum testdata: Ozone by MonthKruskal-Wallis chi-squared = 29.267, df = 4, p-value = 6.901e-06fre1wy> kt <- free1way(Ozone ~ Month, data = airquality)fre1wy> print(kt)5-sample Kruskal-Wallis test against proportional odds alternativesdata: Ozone by Month (5, 6, 7, 8, 9)Perm chi-squared = 29.267, df = 4, p-value = 6.901e-06alternative hypothesis: true log-odds ratio is not equal to 0fre1wy> # log-odds ratios for comparison with controlfre1wy> coef(kt)Month6 Month7 Month8 Month90.8123639 2.5281594 2.3825979 0.7513235fre1wy> # Wald inferencefre1wy> 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.2073Month7 2.52816 0.53153 4.75639 0.0000Month8 2.38260 0.53370 4.46432 0.0000Month9 0.75132 0.47177 1.59256 0.1113fre1wy> confint(kt, test = "Wald")2.5 % 97.5 %Month6 -0.4501567 2.074884Month7 1.4863809 3.569938Month8 1.3365705 3.428625Month9 -0.1733322 1.675979fre1wy> ## Friedman testfre1wy> 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 identicalfre1wy> friedman.test(RoundingTimes)Friedman rank sum testdata: RoundingTimesFriedman chi-squared = 11.143, df = 2, p-value = 0.003805fre1wy> ft <- free1way(time ~ me | id, data = d)fre1wy> print(ft)Stratified 3-sample Kruskal-Wallis test against proportional oddsalternativesdata: time by me (Round Out, Narrow Angle, Wide Angle)stratified by idPerm chi-squared = 11.143, df = 2, p-value = 0.003805alternative hypothesis: true log-odds ratio is not equal to 0fre1wy> coef(ft)meNarrow Angle meWide Angle-1.265280 -3.289198fre1wy> # Wald inferencefre1wy> 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.0638meWide Angle -3.28920 0.79367 -4.14428 0.0000fre1wy> confint(ft, test = "Wald")2.5 % 97.5 %meNarrow Angle -2.603222 0.07266103meWide Angle -4.844765 -1.73363038fre1wy> ## McNemar testfre1wy> ## paired binary observationsfre1wy> example(mcnemar.test, echo = FALSE)fre1wy> # set-up data frame with survey outcomes for votersfre1wy> 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 relevantfre1wy> mcn <- free1way(xtabs(~ performance + survey + voter, data = x))fre1wy> # same result as mcnemar.test w/o continuity correctionfre1wy> print(mcn)Stratified 2-sample Wilcoxon test against proportional oddsalternativesdata: performance by survey (c("1st", "2nd"))stratified by voterPerm Z = 4.166, p-value = 3.099e-05alternative hypothesis: true log-odds ratio is not equal to 0fre1wy> # X^2 statisticfre1wy> summary(mcn, test = "Permutation")$statistic^2Perm Z17.35593fre1wy> mcnemar.test(Performance, correct = FALSE)McNemar's Chi-squared testdata: PerformanceMcNemar's chi-squared = 17.356, df = 1, p-value = 3.099e-05fre1wy> # Wald inferencefre1wy> 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 0fre1wy> confint(mcn, test = "Wald")2.5 % 97.5 %survey2nd 0.7376684 1.487484fre1wy> ## Mantel-Haenszel test w/o continuity correction,fre1wy> ## Departments are blocksfre1wy> mantelhaen.test(UCBAdmissions, correct = FALSE)Mantel-Haenszel chi-squared test without continuity correctiondata: UCBAdmissionsMantel-Haenszel X-squared = 1.5246, df = 1, p-value = 0.2169alternative hypothesis: true common odds ratio is not equal to 195 percent confidence interval:0.7719074 1.0603298sample estimates:common odds ratio0.9046968fre1wy> mh <- free1way(UCBAdmissions)fre1wy> print(mh)Stratified 2-sample Wilcoxon test against proportional oddsalternativesdata: Admit by Gender (c("Male", "Female"))stratified by DeptPerm Z = -1.2347, p-value = 0.2169alternative hypothesis: true log-odds ratio is not equal to 0fre1wy> # common odds-ratio, with score intervalfre1wy> exp(coef(mh))GenderFemale0.904955fre1wy> exp(confint(mh, test = "Rao"))2.5 % 97.5 %GenderFemale 0.7723922 1.06027fre1wy> # looking at department-specificfre1wy> # confidence intervals for log-odds ratiosfre1wy> # it seems Dept A is out of linefre1wy> apply(UCBAdmissions, MARGIN = 3,fre1wy+ FUN = function(x) confint(free1way(as.table(x))))DeptA 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.4025005fre1wy> ## Mantel-Haenszel test treats variables asfre1wy> ## unordered, free1way allows ordered responsesfre1wy> example(mantelhaen.test, echo = FALSE)fre1wy> # Does distribution of job satisfaction (ordered) depend on incomefre1wy> # in a stratified proportional odds model?fre1wy> # Job Satisfaction is second in array but needs to be firstfre1wy> # for free1way to treat it as ordered responsefre1wy> 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.8470Income15000-25000 1.45998 0.62561 2.33369 0.0196Income>25000 1.23282 0.65693 1.87663 0.0606fre1wy> ## End(No test)fre1wy>fre1wy>> example(power.free1way.test, run.donttest = TRUE)pwr.1.> ## No test:pwr.1.>pwr.1.> ## make example reproduciblepwr.1.> set.seed(29)pwr.1.> ## sample from proportional odds model with 1:2 allocationpwr.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 interfacepwr.1.> wilcox.test(y ~ groups, data = x, exact = FALSE, correct = FALSE)$p.value[1] 0.008679882pwr.1.> # Identical p-value obtained from a proportional-odds modelpwr.1.> summary(free1way(y ~ groups, data = x), test = "Permutation")$p.valuePerm Z0.008679882pwr.1.> # approximate power for this testpwr.1.> power.free1way.test(n = 15, delta = log(3), alloc_ratio = 2)2-sample Wilcoxon test against proportional odds alternativesn = 15Total sample size = 15 (Control) + 30 (B) = 45power = 0.4879209sig.level = 0.05Standard error = 0.5693603log-odds ratio = 1.098612NOTE: 'n' is sample size in control grouppwr.1.> ## End(No test)pwr.1.>pwr.1.>> example(ppplot, run.donttest = TRUE)ppplot> ## No test:ppplot>ppplot> ## make example reproducibleppplot> set.seed(29)ppplot> ## well-fitting logistic modelppplot> 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 ratioppplot> coef(ft <- free1way(y ~ groups, data = nd))groupsG21.668854ppplot> # the model-based probability-probability curveppplot> prb <- 1:99 / 100ppplot> 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 bandppplot> 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 intervalppplot> coef(ft, what = "AUC")groupsG20.75467ppplot> confint(ft, test = "Rao", what = "AUC")2.5 % 97.5 %groupsG2 0.6490714 0.8381977ppplot> ## ill-fitting normal modelppplot> 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 modelppplot> coef(free1way(y ~ groups, data = nd, link = "probit"))groupsG20.6659502ppplot> ## 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 freedomResidual deviance: 3311.3 on 68 degrees of freedomAIC: 489.97Number 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 stripStripS01 S02 S03 S04 S05 S06 S07 S08 S09 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S204 4 4 3 3 3 2 2 2 2 3 2 2 2 2 4 4 3 3 3S213>> ## 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 ' ' 1Residual standard error: 4.67 on 57 degrees of freedomNumber of iterations to convergence: 5Achieved 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 ' ' 1Residual standard error: 1.11 on 17 degrees of freedomNumber of iterations to convergence: 8Achieved 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.47.8322873 4.0077676 1.9114016 0.76242114 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: MatrixSubjType 1 2 3 4 5 6 7 8 9T1 4 4 4 4 4 4 4 4 4T2 4 4 4 4 4 4 4 4 4T3 4 4 4 4 4 4 4 4 4T4 4 4 4 4 4 4 4 4 44 x 9 sparse Matrix of class "dgCMatrix"SubjType 1 2 3 4 5 6 7 8 9T1 . 1 . . . 1 . . .T2 . . . . . . . 1 1T3 . . 1 1 1 . 1 1 .T4 . 1 . . . . . . .25 x 4 sparse Matrix of class "dgCMatrix"outerinner A B C Ea . . . 3b . . . 6c . . 5 .d . 1 . .e . . 3 .f . . 5 .g . . 3 .h 6 . . .i 8 . . .j . 1 . .k . . . 6l . 2 . .m . . . 7n . . 3 .o . 3 . .p . 3 . .q . . . 7r . . . 2s . 6 . .t . . . 5u 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__>>