diff -Nru r-cran-mass-7.3-59/ChangeLog r-cran-mass-7.3-60/ChangeLog --- r-cran-mass-7.3-59/ChangeLog 2023-03-07 08:04:00.000000000 +0000 +++ r-cran-mass-7.3-60/ChangeLog 2023-04-07 08:40:16.000000000 +0000 @@ -1,3 +1,12 @@ +7.3-59 + man/{bacteria,epil,gamma.shape,housing,summaey.negbin,theta.md}.Rd: + Use IGNORE_RDIFF markup for change to summary.glm in R 4.3.0. + Back to depending on R (>= 4.0). + +7.3-58.4 (released 2023-04-05) + Emergency fix for summary.glm late changes in 4.3.0. + Depends on R(>= 4.3) and has updated reference output. + 7.3-58.3 (released 2023-03-07) Rename internal function biplot.bdr to biplotBDR Update inst/CITATION. diff -Nru r-cran-mass-7.3-59/debian/changelog r-cran-mass-7.3-60/debian/changelog --- r-cran-mass-7.3-59/debian/changelog 2023-04-29 15:12:42.000000000 +0000 +++ r-cran-mass-7.3-60/debian/changelog 2023-06-29 02:06:23.000000000 +0000 @@ -1,8 +1,22 @@ -r-cran-mass (7.3-59-1.2004.0) focal; urgency=medium +r-cran-mass (7.3-60-2.2004.0) focal; urgency=medium * Compilation for Ubuntu 20.04.6 LTS - -- Michael Rutter Sat, 29 Apr 2023 15:12:42 +0000 + -- Michael Rutter Thu, 29 Jun 2023 02:06:23 +0000 + +r-cran-mass (7.3-60-2) unstable; urgency=medium + + * Rebuilding for unstable following bookworm release + + * debian/control: Set Build-Depends: to current R version + + -- Dirk Eddelbuettel Wed, 21 Jun 2023 06:50:23 -0500 + +r-cran-mass (7.3-60-1) experimental; urgency=medium + + * New upstream release (into 'experimental' while Debian is frozen) + + -- Dirk Eddelbuettel Mon, 08 May 2023 17:35:19 -0500 r-cran-mass (7.3-59-1) experimental; urgency=medium diff -Nru r-cran-mass-7.3-59/debian/control r-cran-mass-7.3-60/debian/control --- r-cran-mass-7.3-59/debian/control 2023-04-29 15:12:42.000000000 +0000 +++ r-cran-mass-7.3-60/debian/control 2023-06-29 02:06:23.000000000 +0000 @@ -2,7 +2,7 @@ Section: gnu-r Priority: optional Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper, r-base-dev (>= 4.3.0), dh-r +Build-Depends: debhelper, r-base-dev (>= 4.3.1), dh-r Standards-Version: 4.6.2 Vcs-Browser: https://salsa.debian.org/edd/r-cran-mass Vcs-Git: https://salsa.debian.org/edd/r-cran-mass.git diff -Nru r-cran-mass-7.3-59/DESCRIPTION r-cran-mass-7.3-60/DESCRIPTION --- r-cran-mass-7.3-59/DESCRIPTION 2023-04-21 10:54:23.000000000 +0000 +++ r-cran-mass-7.3-60/DESCRIPTION 2023-05-04 07:32:21.000000000 +0000 @@ -1,9 +1,9 @@ Package: MASS Priority: recommended -Version: 7.3-59 -Date: 2023-04-03 -Revision: $Rev: 3618 $ -Depends: R (>= 4.2.0), grDevices, graphics, stats, utils +Version: 7.3-60 +Date: 2023-05-02 +Revision: $Rev: 3621 $ +Depends: R (>= 4.0), grDevices, graphics, stats, utils Imports: methods Suggests: lattice, nlme, nnet, survival Authors@R: c(person("Brian", "Ripley", role = c("aut", "cre", "cph"), @@ -24,7 +24,7 @@ URL: http://www.stats.ox.ac.uk/pub/MASS4/ Contact: NeedsCompilation: yes -Packaged: 2023-04-03 15:09:54 UTC; ripley +Packaged: 2023-05-02 16:42:41 UTC; ripley Author: Brian Ripley [aut, cre, cph], Bill Venables [ctb], Douglas M. Bates [ctb], @@ -33,4 +33,4 @@ David Firth [ctb] Maintainer: Brian Ripley Repository: CRAN -Date/Publication: 2023-04-21 10:54:22 UTC +Date/Publication: 2023-05-04 07:32:21 UTC diff -Nru r-cran-mass-7.3-59/man/beav1.Rd r-cran-mass-7.3-60/man/beav1.Rd --- r-cran-mass-7.3-59/man/beav1.Rd 2013-07-18 09:32:11.000000000 +0000 +++ r-cran-mass-7.3-60/man/beav1.Rd 2023-05-02 16:26:42.000000000 +0000 @@ -74,7 +74,7 @@ stemp <- as.vector(temp - alpha*lag(temp, -1)) sX <- X[-1, ] - alpha * X[-115,] beav1.ls <- lm(stemp ~ -1 + sX, na.action = na.omit) -summary(beav1.ls, cor = FALSE) +summary(beav1.ls, correlation = FALSE) rm(temp, activ) } \keyword{datasets} diff -Nru r-cran-mass-7.3-59/man/beav2.Rd r-cran-mass-7.3-60/man/beav2.Rd --- r-cran-mass-7.3-59/man/beav2.Rd 2022-07-27 06:56:41.000000000 +0000 +++ r-cran-mass-7.3-60/man/beav2.Rd 2023-05-02 16:35:53.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/beav2.Rd -% copyright (C) 1994-9 W. N. Venables and B. D. Ripley +% copyright (C) 1994-2023 W. N. Venables and B. D. Ripley % \name{beav2} \alias{beav2} @@ -68,7 +68,7 @@ ## IGNORE_RDIFF_BEGIN library(nlme) # for gls and corAR1 -beav2.gls <- gls(temp ~ activ, data = beav2, corr = corAR1(0.8), +beav2.gls <- gls(temp ~ activ, data = beav2, correlation = corAR1(0.8), method = "ML") summary(beav2.gls) summary(update(beav2.gls, subset = 6:100)) diff -Nru r-cran-mass-7.3-59/man/boxcox.Rd r-cran-mass-7.3-60/man/boxcox.Rd --- r-cran-mass-7.3-59/man/boxcox.Rd 2012-05-15 23:17:48.000000000 +0000 +++ r-cran-mass-7.3-60/man/boxcox.Rd 2023-05-02 16:35:36.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/boxcox.Rd -% copyright (C) 1994-2005 W. N. Venables and B. D. Ripley +% copyright (C) 1994-2023 W. N. Venables and B. D. Ripley % \name{boxcox} \alias{boxcox} @@ -62,10 +62,10 @@ } \examples{ boxcox(Volume ~ log(Height) + log(Girth), data = trees, - lambda = seq(-0.25, 0.25, length = 10)) + lambda = seq(-0.25, 0.25, length.out = 10)) boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, - lambda = seq(-0.05, 0.45, len = 20)) + lambda = seq(-0.05, 0.45, length.out = 20)) } \keyword{regression} \keyword{models} diff -Nru r-cran-mass-7.3-59/man/epil.Rd r-cran-mass-7.3-60/man/epil.Rd --- r-cran-mass-7.3-59/man/epil.Rd 2023-04-03 15:03:10.000000000 +0000 +++ r-cran-mass-7.3-60/man/epil.Rd 2023-05-02 16:29:21.000000000 +0000 @@ -58,7 +58,7 @@ \examples{ ## IGNORE_RDIFF_BEGIN summary(glm(y ~ lbase*trt + lage + V4, family = poisson, - data = epil), cor = FALSE) + data = epil), correlation = FALSE) ## IGNORE_RDIFF_END epil2 <- epil[epil$period == 1, ] epil2["period"] <- rep(0, 59); epil2["y"] <- epil2["base"] @@ -75,7 +75,7 @@ dimnames = list(NULL, c("placebo-base", "drug-placebo"))) ## IGNORE_RDIFF_BEGIN summary(glm(y ~ pred + factor(subject) + offset(log(time)), - family = poisson, data = epil3), cor = FALSE) + family = poisson, data = epil3), correlation = FALSE) ## IGNORE_RDIFF_END summary(glmmPQL(y ~ lbase*trt + lage + V4, diff -Nru r-cran-mass-7.3-59/man/housing.Rd r-cran-mass-7.3-60/man/housing.Rd --- r-cran-mass-7.3-59/man/housing.Rd 2023-04-03 15:04:21.000000000 +0000 +++ r-cran-mass-7.3-60/man/housing.Rd 2023-05-02 16:35:22.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/housing.Rd -% copyright (C) 1999 W. N. Venables and B. D. Ripley +% copyright (C) 1999-2023 W. N. Venables and B. D. Ripley % \name{housing} \alias{housing} @@ -53,14 +53,14 @@ house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family = poisson, data = housing) ## IGNORE_RDIFF_BEGIN -summary(house.glm0, cor = FALSE) +summary(house.glm0, correlation = FALSE) ## IGNORE_RDIFF_END addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test = "Chisq") house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) ## IGNORE_RDIFF_BEGIN -summary(house.glm1, cor = FALSE) +summary(house.glm1, correlation = FALSE) ## IGNORE_RDIFF_END 1 - pchisq(deviance(house.glm1), house.glm1$df.residual) diff -Nru r-cran-mass-7.3-59/man/logtrans.Rd r-cran-mass-7.3-60/man/logtrans.Rd --- r-cran-mass-7.3-59/man/logtrans.Rd 2012-05-15 23:17:48.000000000 +0000 +++ r-cran-mass-7.3-60/man/logtrans.Rd 2023-05-02 16:35:14.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/logtrans.Rd -% copyright (C) 1994-2004 W. N. Venables and B. D. Ripley +% copyright (C) 1994-2023 W. N. Venables and B. D. Ripley % \name{logtrans} \alias{logtrans} @@ -71,7 +71,7 @@ } \examples{ logtrans(Days ~ Age*Sex*Eth*Lrn, data = quine, - alpha = seq(0.75, 6.5, len=20)) + alpha = seq(0.75, 6.5, length.out = 20)) } \keyword{regression} \keyword{models} diff -Nru r-cran-mass-7.3-59/man/OME.Rd r-cran-mass-7.3-60/man/OME.Rd --- r-cran-mass-7.3-59/man/OME.Rd 2022-07-27 06:56:41.000000000 +0000 +++ r-cran-mass-7.3-60/man/OME.Rd 2023-05-02 16:36:01.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/OME.Rd -% copyright (C) 1994-2014 W. N. Venables and B. D. Ripley +% copyright (C) 1994-2023 W. N. Venables and B. D. Ripley % \name{OME} \alias{OME} @@ -121,7 +121,7 @@ summary(lm(L75 ~ Noise/Age, data = OMEif, na.action = na.omit)) summary(lm(L75 ~ Noise/(Age + OME), data = OMEif, subset = (Age >= 30 & Age <= 60), - na.action = na.omit), cor = FALSE) + na.action = na.omit), correlation = FALSE) # Or fit by weighted least squares fpl75 <- deriv(~ sqrt(n)*(r/n - 0.5 - 0.5/(1 + exp(-(x-L75)/scal))), diff -Nru r-cran-mass-7.3-59/man/petrol.Rd r-cran-mass-7.3-60/man/petrol.Rd --- r-cran-mass-7.3-59/man/petrol.Rd 2022-07-27 06:56:41.000000000 +0000 +++ r-cran-mass-7.3-60/man/petrol.Rd 2023-05-02 16:34:49.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/petrol.Rd -% copyright (C) 1994-9 W. N. Venables and B. D. Ripley +% copyright (C) 1994-2023 W. N. Venables and B. D. Ripley % \name{petrol} \alias{petrol} @@ -59,7 +59,7 @@ pet3.lme <- lme(Y ~ SG + VP + V10 + EP, random = ~ 1 | No, data = Petrol) pet3.lme <- update(pet3.lme, method = "ML") -pet4.lme <- update(pet3.lme, fixed = Y ~ V10 + EP) +pet4.lme <- update(pet3.lme, fixed. = Y ~ V10 + EP) anova(pet4.lme, pet3.lme) } \keyword{datasets} diff -Nru r-cran-mass-7.3-59/man/plot.profile.Rd r-cran-mass-7.3-60/man/plot.profile.Rd --- r-cran-mass-7.3-59/man/plot.profile.Rd 2012-05-15 23:17:48.000000000 +0000 +++ r-cran-mass-7.3-60/man/plot.profile.Rd 2023-05-02 16:34:39.000000000 +0000 @@ -1,5 +1,5 @@ % file MASS/man/plot.profile.Rd -% copyright (C) 1999-2008 W. N. Venables and B. D. Ripley +% copyright (C) 1999-2023 W. N. Venables and B. D. Ripley % \name{plot.profile} \alias{plot.profile} @@ -45,7 +45,7 @@ ## a version of example(profile.nls) from R >= 2.8.0 fm1 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD) -pr1 <- profile(fm1, alpha = 0.1) +pr1 <- profile(fm1, alphamax = 0.1) MASS:::plot.profile(pr1) pairs(pr1) # a little odd since the parameters are highly correlated diff -Nru r-cran-mass-7.3-59/MD5 r-cran-mass-7.3-60/MD5 --- r-cran-mass-7.3-59/MD5 2023-04-21 10:54:23.000000000 +0000 +++ r-cran-mass-7.3-60/MD5 2023-05-04 07:32:21.000000000 +0000 @@ -1,9 +1,9 @@ -5a8536600acf8f1f058152f15a0e440c *ChangeLog -ff57f171b0da21d94198ff2501793d8a *DESCRIPTION +c32327b0d8b470b18771ad69bbec57fb *ChangeLog +4a77c1020501bbff680b80e1d4523554 *DESCRIPTION 35aff05a505ecf7e81e0473767794ca9 *INDEX c7acdc0fa828f781a0a5586ab9d4fa1b *LICENCE.note 5ae9f8060ce744604b563ddc755782f5 *NAMESPACE -aed71ed40bb9bfdd7d19edfdbd163087 *R/add.R +f1baa4375c9dab6d182ede47c4ea433b *R/add.R 667178114b2cc60142fed7fc16c00f90 *R/area.R 2b449dd3fde83489799ef5f72bc8ead0 *R/boxcox.R a22a1823498627d344110a4deb10c0ac *R/confint.R @@ -22,22 +22,22 @@ feed05160fcb53f8d07c559e1afa9873 *R/hubers.R b714f739c6c9dcb33729bb2baa532a0c *R/isoMDS.R 0e20d824adf45263385c0245060ecc14 *R/kde2d.R -6ea9e10aca09bbe90c746c613aaa5cd8 *R/lda.R +26e63a4814e5daea7e5140f586b7aea8 *R/lda.R 797d778deac7bbe14bbc77b9ad96379d *R/lm.gls.R 4cb17b20ccacf80d37cc726736244907 *R/lm.ridge.R 411b677b9b0b478d204518b50f363332 *R/loglm.R a03b55f77943de6230a6e7e61b650c56 *R/logtrans.R -a70b5f159d36abb0741a144779232d92 *R/lqs.R +355f8fd440cc127d599efdf1373a9391 *R/lqs.R 865f29d0ba0995ca13543bf44f391e19 *R/mca.R f7d208bdf6bd06e047c5ffff547dac7f *R/misc.R c00a7143684b5b3bca4c0d0288259c21 *R/mvrnorm.R 1de4b445224bc4bc1c977414bf744f7e *R/neg.bin.R -a2717441255135a2687efee2f29e6ea2 *R/negbin.R +2416ff9c2618ebc7758a075fd0856cfc *R/negbin.R f4deec8fd9e6c10ee006abce09a7b246 *R/negexp.R d46385a583ad1e374fe3470d9c9b4a7e *R/parcoord.R -c741f8f80f1d5574a17703fee6465f1c *R/polr.R +3a4c3754f9f46b4bef65b3d83b0beef0 *R/polr.R 54a2b2d504391f36f4d9211fb361a26d *R/profiles.R -a5771ffbdd84fc93ad1993f4ca3af668 *R/qda.R +dc943bfe2e1f589a6cfc1193da709248 *R/qda.R 4d8329c5c96150724dadebf17fb0ff3f *R/rlm.R 7d35a77423de76adcd49a981816c30eb *R/rms.curv.R c301b0a87d4188e4e9dd9f361aefb3f9 *R/sammon.R @@ -170,7 +170,7 @@ a6a2cea128307c1138bde4c046fcfb5c *man/MASS-internal.Rd fa8c2c8e00a52e7731604be22fafcec6 *man/Melanoma.Rd 9535c4de30b076e0ba2ec4d1f868b9ac *man/Null.Rd -9071e3f3ca15b70da7dfed264675d965 *man/OME.Rd +c2730a9f8d03251c4acc3cd7e2ffc3f9 *man/OME.Rd b3ee04fa917bb291c47d52830d77ffdc *man/Pima.tr.Rd a7ba31f8b57db3cbb58b636ec45ad19d *man/Rabbit.Rd c542b9d655977d4279f267e714581841 *man/Rubber.Rd @@ -191,11 +191,11 @@ 776675b37ae67af99bd7aa16baa121d4 *man/bacteria.Rd 5472d70fb96f35a995e891ad9ef9a777 *man/bandwidth.nrd.Rd e79ea25d60ba827b71a57f7458f733d2 *man/bcv.Rd -61fbfe12f436f23b478400e0c085731d *man/beav1.Rd -9c16e1752eba19b699081c3f7f82db01 *man/beav2.Rd +e448793ca52ceee0bce4c823ad2ccfc9 *man/beav1.Rd +086cdd567cff1d6cb92b8606c4c94b4f *man/beav2.Rd 55b0db236ec4393bae165ab49c92ee48 *man/biopsy.Rd 2f6ba17e14d2daa4afc541d3445190ad *man/birthwt.Rd -bebc2ad312668d479d67adbb5fd3777d *man/boxcox.Rd +e61cee7038d3c48c7d530f2e6e5e37bf *man/boxcox.Rd 1a1233168369b0abcf4e2d207d99dc72 *man/cabbages.Rd 9a9a2f5323a13a701a74a514785cd771 *man/caith.Rd 726868057d0e31cad05e01c7b379573c *man/cats.Rd @@ -216,7 +216,7 @@ 32d46414b81cc9c21663b07a5b3e24a1 *man/drivers.Rd 989cb309d961aca2b8d9d54df183223d *man/dropterm.Rd 57e48204af505d9d279440f9a75fe46d *man/eagles.Rd -2546dca2ae4c53fe9e89b67bebfcfcdb *man/epil.Rd +d7271f3dc835de0337cbc886ff587e62 *man/epil.Rd df9ff3b3ed9842a6c2d767dfb14233d8 *man/eqscplot.Rd d6c9433b3ee2a9fe920dcedecbefdf06 *man/farms.Rd 7d46e561298348eafcb2464adcfa66b2 *man/fgl.Rd @@ -236,7 +236,7 @@ b07145bdfcdcf3d621483562aab95335 *man/glmmPQL.Rd e65737579ca1cd2810198d72f0260e21 *man/hills.Rd 5440204f0acdca119adda063b91d9e17 *man/hist.scott.Rd -16e4f5a728d37e3d423566bba03d63f3 *man/housing.Rd +05e046d8196bef63e85ab93406219262 *man/housing.Rd 3ff54e8dfc1c1b949d3a33fcbd916fee *man/huber.Rd 7cb24fac3e55b03ac894c3be943287bf *man/hubers.Rd 606105b6a9c526833a3116619a07d9c9 *man/immer.Rd @@ -249,7 +249,7 @@ e852b476aaca59a57b3a5deb4523977e *man/lm.ridge.Rd 9a82b09d10ce38560bc5abc66b5f5dfe *man/loglm.Rd 82b7f4c063c79b490fe7e02d261ddbb7 *man/loglm1.Rd -4ccf847b74d8eabfa2608f644a1a5f5c *man/logtrans.Rd +caebc0a782ad8bd34fe2038e13ccff77 *man/logtrans.Rd 59c78564e95b8d79447af0f21ed60bc9 *man/lqs.Rd fbec6bc6908ad65e919c06945c4d6b42 *man/mammals.Rd 9fbb7d21860ae43c710d73e78d68741f *man/mca.Rd @@ -269,11 +269,11 @@ 00698a2eb415969b83999b0ff995cced *man/painters.Rd 71a3eccdc2ae99c4b1120d6a6855b14a *man/pairs.lda.Rd 0020b25127cad664be91aef4c72996b4 *man/parcoord.Rd -f65ceb9db76de192862127a7a63daab3 *man/petrol.Rd +d33e861e12b7a626a82be4d99b0313d1 *man/petrol.Rd 61889526fa86b85807a2d64d15c06098 *man/phones.Rd 815f79defe5a3e762f3b45ce1da2b677 *man/plot.lda.Rd 376f478df44a0d407e103f23bff67857 *man/plot.mca.Rd -3acff5e8ecf57dc83a632ad19f645597 *man/plot.profile.Rd +4181ccf11c5a5b1e649183eef39d1c53 *man/plot.profile.Rd 3a80cfc9c47920bd19b77000253ade06 *man/polr.Rd f10d097bc2a4fbb9a4904db9fa2d512e *man/predict.glmmPQL.Rd 6eca3725d70f18fd1a5d257d0d1e50db *man/predict.lda.Rd @@ -324,7 +324,8 @@ f47fd24f72996b91d2a1618cd105c908 *src/MASS.c 2c7b72ca038d349935bcb4f60235bcd7 *src/lqs.c 7b1360150b8598b576e82534b12345a2 *tests/BankWages.rda -26860fec833859bcf689559b80c5deda *tests/Examples/MASS-Ex.Rout.save +42554314ed77c341a224e20a018906ca *tests/Examples/MASS-Ex.Rout.save +15150ef489298cb9d9fdf0c94b6474fd *tests/MASS-Ex.Rout a309c9e287f72ae6fc239b0495d5685e *tests/confint.R e059b0451b55d1e8cb8ebf51ecce76ea *tests/cov.mcd.R a0854fee9b36691ca1d0880da31d6f79 *tests/fitdistr.R diff -Nru r-cran-mass-7.3-59/R/add.R r-cran-mass-7.3-60/R/add.R --- r-cran-mass-7.3-59/R/add.R 2014-08-21 07:57:21.000000000 +0000 +++ r-cran-mass-7.3-60/R/add.R 2023-05-02 16:24:05.000000000 +0000 @@ -1,5 +1,5 @@ # file MASS/R/add.R -# copyright (C) 1994-2008 W. N. Venables and B. D. Ripley +# copyright (C) 1994-2023 W. N. Venables and B. D. Ripley # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -176,7 +176,7 @@ fob <- list(call = oc, terms=Terms) class(fob) <- class(object) x <- model.matrix(Terms, model.frame(fob, xlev = object$xlevels), - contrasts = object$contrasts) + contrasts.arg = object$contrasts) n <- nrow(x) oldn <- length(object$residuals) y <- object$y diff -Nru r-cran-mass-7.3-59/R/lda.R r-cran-mass-7.3-60/R/lda.R --- r-cran-mass-7.3-59/R/lda.R 2015-02-01 14:56:42.000000000 +0000 +++ r-cran-mass-7.3-60/R/lda.R 2023-05-02 16:39:19.000000000 +0000 @@ -1,5 +1,5 @@ # file MASS/R/lda.R -# copyright (C) 1994-2013 W. N. Venables and B. D. Ripley +# copyright (C) 1994-2023 W. N. Venables and B. D. Ripley # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -223,7 +223,7 @@ if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, newdata) } - x <- model.matrix(Terms, newdata, contrasts = object$contrasts) + x <- model.matrix(Terms, newdata, contrasts.arg = object$contrasts) xint <- match("(Intercept)", colnames(x), nomatch = 0L) if(xint > 0L) x <- x[, -xint, drop = FALSE] } else { # matrix or data-frame fit diff -Nru r-cran-mass-7.3-59/R/lqs.R r-cran-mass-7.3-60/R/lqs.R --- r-cran-mass-7.3-59/R/lqs.R 2020-09-06 09:09:59.000000000 +0000 +++ r-cran-mass-7.3-60/R/lqs.R 2023-05-02 16:39:37.000000000 +0000 @@ -1,5 +1,5 @@ # file lqs/R/lqs.R -# copyright (C) 1998-2020 B. D. Ripley +# copyright (C) 1998-2023 B. D. Ripley # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -192,7 +192,7 @@ m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels) if(!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m) - X <- model.matrix(Terms, m, contrasts = object$contrasts) + X <- model.matrix(Terms, m, contrasts.arg = object$contrasts) drop(X %*% object$coefficients) } diff -Nru r-cran-mass-7.3-59/R/negbin.R r-cran-mass-7.3-60/R/negbin.R --- r-cran-mass-7.3-59/R/negbin.R 2015-02-01 11:18:38.000000000 +0000 +++ r-cran-mass-7.3-60/R/negbin.R 2023-05-02 16:23:56.000000000 +0000 @@ -1,5 +1,5 @@ # file MASS/R/negbin.R -# copyright (C) 1994-2014 W. N. Venables and B. D. Ripley +# copyright (C) 1994-2023 W. N. Venables and B. D. Ripley # anova.negbin <- function(object, ..., test = "Chisq") { @@ -119,7 +119,7 @@ glm.fitter <- stats::glm.fit } if(control$trace > 1) message("Initial fit:") - fit <- glm.fitter(x = X, y = Y, w = w, start = start, + fit <- glm.fitter(x = X, y = Y, weights = w, start = start, etastart = etastart, mustart = mustart, offset = offset, family = fam0, control = list(maxit=control$maxit, @@ -143,7 +143,7 @@ while((iter <- iter + 1) <= control$maxit && (abs(Lm0 - Lm)/d1 + abs(del)/d2) > control$epsilon) { eta <- g(mu) - fit <- glm.fitter(x = X, y = Y, w = w, etastart = + fit <- glm.fitter(x = X, y = Y, weights = w, etastart = eta, offset = offset, family = fam, control = list(maxit=control$maxit, epsilon = control$epsilon, diff -Nru r-cran-mass-7.3-59/R/polr.R r-cran-mass-7.3-60/R/polr.R --- r-cran-mass-7.3-59/R/polr.R 2015-08-30 07:29:59.000000000 +0000 +++ r-cran-mass-7.3-60/R/polr.R 2023-05-02 16:23:46.000000000 +0000 @@ -1,5 +1,5 @@ # file MASS/R/polr.R -# copyright (C) 1994-2013 W. N. Venables and B. D. Ripley +# copyright (C) 1994-2023 W. N. Venables and B. D. Ripley # Use of transformed intercepts contributed by David Firth # # This program is free software; you can redistribute it and/or modify @@ -224,7 +224,7 @@ xlev = object$xlevels) if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m) - X <- model.matrix(Terms, m, contrasts = object$contrasts) + X <- model.matrix(Terms, m, contrasts.arg = object$contrasts) xint <- match("(Intercept)", colnames(X), nomatch=0L) if(xint > 0L) X <- X[, -xint, drop=FALSE] n <- nrow(X) diff -Nru r-cran-mass-7.3-59/R/qda.R r-cran-mass-7.3-60/R/qda.R --- r-cran-mass-7.3-59/R/qda.R 2013-07-18 09:32:11.000000000 +0000 +++ r-cran-mass-7.3-60/R/qda.R 2023-05-02 16:39:11.000000000 +0000 @@ -1,5 +1,5 @@ # file MASS/R/qda.R -# copyright (C) 1994-2013 W. N. Venables and B. D. Ripley +# copyright (C) 1994-2023 W. N. Venables and B. D. Ripley # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -204,7 +204,7 @@ xlev = object$xlevels) } x <- model.matrix(delete.response(Terms), newdata, - contrasts = object$contrasts) + contrasts.arg = object$contrasts) xint <- match("(Intercept)", colnames(x), nomatch=0L) if(xint > 0) x <- x[, -xint, drop=FALSE] if(method == "looCV") g <- model.response(newdata) diff -Nru r-cran-mass-7.3-59/tests/Examples/MASS-Ex.Rout.save r-cran-mass-7.3-60/tests/Examples/MASS-Ex.Rout.save --- r-cran-mass-7.3-59/tests/Examples/MASS-Ex.Rout.save 2023-04-03 15:06:18.000000000 +0000 +++ r-cran-mass-7.3-60/tests/Examples/MASS-Ex.Rout.save 2023-05-02 16:41:59.000000000 +0000 @@ -1,5 +1,5 @@ -R version 4.3.0 alpha (2023-04-03 r84149) +R Under development (unstable) (2023-05-02 r84376) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: aarch64-apple-darwin22.4.0 (64-bit) @@ -203,7 +203,7 @@ > summary(lm(L75 ~ Noise/(Age + OME), data = OMEif, + subset = (Age >= 30 & Age <= 60), -+ na.action = na.omit), cor = FALSE) ++ na.action = na.omit), correlation = FALSE) Call: lm(formula = L75 ~ Noise/(Age + OME), data = OMEif, subset = (Age >= @@ -759,7 +759,7 @@ > stemp <- as.vector(temp - alpha*lag(temp, -1)) > sX <- X[-1, ] - alpha * X[-115,] > beav1.ls <- lm(stemp ~ -1 + sX, na.action = na.omit) -> summary(beav1.ls, cor = FALSE) +> summary(beav1.ls, correlation = FALSE) Call: lm(formula = stemp ~ -1 + sX, na.action = na.omit) @@ -857,7 +857,7 @@ > > ## IGNORE_RDIFF_BEGIN > library(nlme) # for gls and corAR1 -> beav2.gls <- gls(temp ~ activ, data = beav2, corr = corAR1(0.8), +> beav2.gls <- gls(temp ~ activ, data = beav2, correlation = corAR1(0.8), + method = "ML") > summary(beav2.gls) Generalized least squares fit by maximum likelihood @@ -978,10 +978,10 @@ > ### ** Examples > > boxcox(Volume ~ log(Height) + log(Girth), data = trees, -+ lambda = seq(-0.25, 0.25, length = 10)) ++ lambda = seq(-0.25, 0.25, length.out = 10)) > > boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, -+ lambda = seq(-0.05, 0.45, len = 20)) ++ lambda = seq(-0.05, 0.45, length.out = 20)) > > > @@ -1491,7 +1491,7 @@ > > ## IGNORE_RDIFF_BEGIN > summary(glm(y ~ lbase*trt + lage + V4, family = poisson, -+ data = epil), cor = FALSE) ++ data = epil), correlation = FALSE) Call: glm(formula = y ~ lbase * trt + lage + V4, family = poisson, @@ -1532,7 +1532,7 @@ + dimnames = list(NULL, c("placebo-base", "drug-placebo"))) > ## IGNORE_RDIFF_BEGIN > summary(glm(y ~ pred + factor(subject) + offset(log(time)), -+ family = poisson, data = epil3), cor = FALSE) ++ family = poisson, data = epil3), correlation = FALSE) Call: glm(formula = y ~ pred + factor(subject) + offset(log(time)), @@ -2341,7 +2341,7 @@ > house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family = poisson, + data = housing) > ## IGNORE_RDIFF_BEGIN -> summary(house.glm0, cor = FALSE) +> summary(house.glm0, correlation = FALSE) Call: glm(formula = Freq ~ Infl * Type * Cont + Sat, family = poisson, @@ -2403,7 +2403,7 @@ > > house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) > ## IGNORE_RDIFF_BEGIN -> summary(house.glm1, cor = FALSE) +> summary(house.glm1, correlation = FALSE) Call: glm(formula = Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + @@ -3154,7 +3154,7 @@ > ### ** Examples > > logtrans(Days ~ Age*Sex*Eth*Lrn, data = quine, -+ alpha = seq(0.75, 6.5, len=20)) ++ alpha = seq(0.75, 6.5, length.out = 20)) > > > @@ -3788,7 +3788,7 @@ > pet3.lme <- lme(Y ~ SG + VP + V10 + EP, + random = ~ 1 | No, data = Petrol) > pet3.lme <- update(pet3.lme, method = "ML") -> pet4.lme <- update(pet3.lme, fixed = Y ~ V10 + EP) +> pet4.lme <- update(pet3.lme, fixed. = Y ~ V10 + EP) > anova(pet4.lme, pet3.lme) Model df AIC BIC logLik Test L.Ratio p-value pet4.lme 1 5 149.6119 156.9406 -69.80594 @@ -3833,7 +3833,7 @@ > > ## a version of example(profile.nls) from R >= 2.8.0 > fm1 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD) -> pr1 <- profile(fm1, alpha = 0.1) +> pr1 <- profile(fm1, alphamax = 0.1) > MASS:::plot.profile(pr1) > pairs(pr1) # a little odd since the parameters are highly correlated > @@ -5124,7 +5124,7 @@ > cleanEx() > options(digits = 7L) > base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n") -Time elapsed: 2.206 0.12 2.438 0 0 +Time elapsed: 2.162 0.095 2.26 0 0 > grDevices::dev.off() null device 1 diff -Nru r-cran-mass-7.3-59/tests/MASS-Ex.Rout r-cran-mass-7.3-60/tests/MASS-Ex.Rout --- r-cran-mass-7.3-59/tests/MASS-Ex.Rout 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mass-7.3-60/tests/MASS-Ex.Rout 2023-05-02 16:40:18.000000000 +0000 @@ -0,0 +1,5150 @@ + +R Under development (unstable) (2023-05-02 r84376) -- "Unsuffered Consequences" +Copyright (C) 2023 The R Foundation for Statistical Computing +Platform: aarch64-apple-darwin22.4.0 (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. + + Natural language support but running in an English locale + +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. + +> pkgname <- "MASS" +> source(file.path(R.home("share"), "R", "examples-header.R")) +> options(warn = 1) +> library('MASS') +> +> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') +> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv') +> cleanEx() +> nameEx("Insurance") +> ### * Insurance +> +> flush(stderr()); flush(stdout()) +> +> ### Name: Insurance +> ### Title: Numbers of Car Insurance claims +> ### Aliases: Insurance +> ### Keywords: datasets +> +> ### ** Examples +> +> ## main-effects fit as Poisson GLM with offset +> glm(Claims ~ District + Group + Age + offset(log(Holders)), ++ data = Insurance, family = poisson) + +Call: glm(formula = Claims ~ District + Group + Age + offset(log(Holders)), + family = poisson, data = Insurance) + +Coefficients: +(Intercept) District2 District3 District4 Group.L Group.Q + -1.810508 0.025868 0.038524 0.234205 0.429708 0.004632 + Group.C Age.L Age.Q Age.C + -0.029294 -0.394432 -0.000355 -0.016737 + +Degrees of Freedom: 63 Total (i.e. Null); 54 Residual +Null Deviance: 236.3 +Residual Deviance: 51.42 AIC: 388.7 +> +> # same via loglm +> loglm(Claims ~ District + Group + Age + offset(log(Holders)), ++ data = Insurance) +Call: +loglm(formula = Claims ~ District + Group + Age + offset(log(Holders)), + data = Insurance) + +Statistics: + X^2 df P(> X^2) +Likelihood Ratio 51.42003 54 0.5745071 +Pearson 48.62933 54 0.6809086 +> +> +> +> cleanEx() +> nameEx("Null") +> ### * Null +> +> flush(stderr()); flush(stdout()) +> +> ### Name: Null +> ### Title: Null Spaces of Matrices +> ### Aliases: Null +> ### Keywords: algebra +> +> ### ** Examples +> +> # The function is currently defined as +> function(M) ++ { ++ tmp <- qr(M) ++ set <- if(tmp$rank == 0L) seq_len(ncol(M)) else -seq_len(tmp$rank) ++ qr.Q(tmp, complete = TRUE)[, set, drop = FALSE] ++ } +function (M) +{ + tmp <- qr(M) + set <- if (tmp$rank == 0L) + seq_len(ncol(M)) + else -seq_len(tmp$rank) + qr.Q(tmp, complete = TRUE)[, set, drop = FALSE] +} +> +> +> +> cleanEx() +> nameEx("OME") +> ### * OME +> +> flush(stderr()); flush(stdout()) +> +> ### Name: OME +> ### Title: Tests of Auditory Perception in Children with OME +> ### Aliases: OME +> ### Keywords: datasets +> +> ### ** Examples +> +> # Fit logistic curve from p = 0.5 to p = 1.0 +> fp1 <- deriv(~ 0.5 + 0.5/(1 + exp(-(x-L75)/scal)), ++ c("L75", "scal"), ++ function(x,L75,scal)NULL) +> nls(Correct/Trials ~ fp1(Loud, L75, scal), data = OME, ++ start = c(L75=45, scal=3)) +Nonlinear regression model + model: Correct/Trials ~ fp1(Loud, L75, scal) + data: OME + L75 scal +44.149 3.775 + residual sum-of-squares: 69.88 + +Number of iterations to convergence: 4 +Achieved convergence tolerance: 7.016e-06 +> nls(Correct/Trials ~ fp1(Loud, L75, scal), ++ data = OME[OME$Noise == "coherent",], ++ start=c(L75=45, scal=3)) +Nonlinear regression model + model: Correct/Trials ~ fp1(Loud, L75, scal) + data: OME[OME$Noise == "coherent", ] + L75 scal +47.993 1.259 + residual sum-of-squares: 30.35 + +Number of iterations to convergence: 5 +Achieved convergence tolerance: 4.895e-06 +> nls(Correct/Trials ~ fp1(Loud, L75, scal), ++ data = OME[OME$Noise == "incoherent",], ++ start = c(L75=45, scal=3)) +Nonlinear regression model + model: Correct/Trials ~ fp1(Loud, L75, scal) + data: OME[OME$Noise == "incoherent", ] + L75 scal +38.87 2.17 + residual sum-of-squares: 23.73 + +Number of iterations to convergence: 11 +Achieved convergence tolerance: 3.846e-06 +> +> # individual fits for each experiment +> +> aa <- factor(OME$Age) +> ab <- 10*OME$ID + unclass(aa) +> ac <- unclass(factor(ab)) +> OME$UID <- as.vector(ac) +> OME$UIDn <- OME$UID + 0.1*(OME$Noise == "incoherent") +> rm(aa, ab, ac) +> OMEi <- OME +> +> library(nlme) +> fp2 <- deriv(~ 0.5 + 0.5/(1 + exp(-(x-L75)/2)), ++ "L75", function(x,L75) NULL) +> dec <- getOption("OutDec") +> options(show.error.messages = FALSE, OutDec=".") +> OMEi.nls <- nlsList(Correct/Trials ~ fp2(Loud, L75) | UIDn, ++ data = OMEi, start = list(L75=45), control = list(maxiter=100)) +> options(show.error.messages = TRUE, OutDec=dec) +> tmp <- sapply(OMEi.nls, function(X) ++ {if(is.null(X)) NA else as.vector(coef(X))}) +> OMEif <- data.frame(UID = round(as.numeric((names(tmp)))), ++ Noise = rep(c("coherent", "incoherent"), 110), ++ L75 = as.vector(tmp), stringsAsFactors = TRUE) +> OMEif$Age <- OME$Age[match(OMEif$UID, OME$UID)] +> OMEif$OME <- OME$OME[match(OMEif$UID, OME$UID)] +> OMEif <- OMEif[OMEif$L75 > 30,] +> summary(lm(L75 ~ Noise/Age, data = OMEif, na.action = na.omit)) + +Call: +lm(formula = L75 ~ Noise/Age, data = OMEif, na.action = na.omit) + +Residuals: + Min 1Q Median 3Q Max +-13.0022 -1.9878 0.3346 2.0229 16.3260 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +(Intercept) 47.73580 0.76456 62.435 < 2e-16 *** +Noiseincoherent -4.87352 1.11247 -4.381 1.92e-05 *** +Noisecoherent:Age -0.02785 0.02349 -1.186 0.237 +Noiseincoherent:Age -0.12219 0.02589 -4.719 4.50e-06 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 3.774 on 196 degrees of freedom + (17 observations deleted due to missingness) +Multiple R-squared: 0.5246, Adjusted R-squared: 0.5173 +F-statistic: 72.09 on 3 and 196 DF, p-value: < 2.2e-16 + +> summary(lm(L75 ~ Noise/(Age + OME), data = OMEif, ++ subset = (Age >= 30 & Age <= 60), ++ na.action = na.omit), correlation = FALSE) + +Call: +lm(formula = L75 ~ Noise/(Age + OME), data = OMEif, subset = (Age >= + 30 & Age <= 60), na.action = na.omit) + +Residuals: + Min 1Q Median 3Q Max +-10.4514 -2.0588 0.0194 1.6827 15.9738 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +(Intercept) 50.21090 1.74482 28.777 < 2e-16 *** +Noiseincoherent -5.97491 2.70148 -2.212 0.02890 * +Noisecoherent:Age -0.09358 0.03586 -2.609 0.01023 * +Noiseincoherent:Age -0.15155 0.04151 -3.651 0.00039 *** +Noisecoherent:OMElow 0.45103 1.07594 0.419 0.67583 +Noiseincoherent:OMElow -0.14075 1.24537 -0.113 0.91021 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 3.7 on 119 degrees of freedom + (17 observations deleted due to missingness) +Multiple R-squared: 0.6073, Adjusted R-squared: 0.5908 +F-statistic: 36.81 on 5 and 119 DF, p-value: < 2.2e-16 + +> +> # Or fit by weighted least squares +> fpl75 <- deriv(~ sqrt(n)*(r/n - 0.5 - 0.5/(1 + exp(-(x-L75)/scal))), ++ c("L75", "scal"), ++ function(r,n,x,L75,scal) NULL) +> nls(0 ~ fpl75(Correct, Trials, Loud, L75, scal), ++ data = OME[OME$Noise == "coherent",], ++ start = c(L75=45, scal=3)) +Nonlinear regression model + model: 0 ~ fpl75(Correct, Trials, Loud, L75, scal) + data: OME[OME$Noise == "coherent", ] + L75 scal +47.798 1.296 + residual sum-of-squares: 91.72 + +Number of iterations to convergence: 5 +Achieved convergence tolerance: 9.302e-06 +> nls(0 ~ fpl75(Correct, Trials, Loud, L75, scal), ++ data = OME[OME$Noise == "incoherent",], ++ start = c(L75=45, scal=3)) +Nonlinear regression model + model: 0 ~ fpl75(Correct, Trials, Loud, L75, scal) + data: OME[OME$Noise == "incoherent", ] + L75 scal +38.553 2.078 + residual sum-of-squares: 60.19 + +Number of iterations to convergence: 8 +Achieved convergence tolerance: 4.55e-06 +> +> # Test to see if the curves shift with age +> fpl75age <- deriv(~sqrt(n)*(r/n - 0.5 - 0.5/(1 + ++ exp(-(x-L75-slope*age)/scal))), ++ c("L75", "slope", "scal"), ++ function(r,n,x,age,L75,slope,scal) NULL) +> OME.nls1 <- ++ nls(0 ~ fpl75age(Correct, Trials, Loud, Age, L75, slope, scal), ++ data = OME[OME$Noise == "coherent",], ++ start = c(L75=45, slope=0, scal=2)) +> sqrt(diag(vcov(OME.nls1))) + L75 slope scal +0.61091761 0.01665916 0.17566450 +> +> OME.nls2 <- ++ nls(0 ~ fpl75age(Correct, Trials, Loud, Age, L75, slope, scal), ++ data = OME[OME$Noise == "incoherent",], ++ start = c(L75=45, slope=0, scal=2)) +> sqrt(diag(vcov(OME.nls2))) + L75 slope scal +0.49553854 0.01348281 0.24453836 +> +> # Now allow random effects by using NLME +> OMEf <- OME[rep(1:nrow(OME), OME$Trials),] +> OMEf$Resp <- with(OME, rep(rep(c(1,0), length(Trials)), ++ t(cbind(Correct, Trials-Correct)))) +> OMEf <- OMEf[, -match(c("Correct", "Trials"), names(OMEf))] +> +> ## Not run: +> ##D ## these fail in R on most platforms +> ##D fp2 <- deriv(~ 0.5 + 0.5/(1 + exp(-(x-L75)/exp(lsc))), +> ##D c("L75", "lsc"), +> ##D function(x, L75, lsc) NULL) +> ##D try(summary(nlme(Resp ~ fp2(Loud, L75, lsc), +> ##D fixed = list(L75 ~ Age, lsc ~ 1), +> ##D random = L75 + lsc ~ 1 | UID, +> ##D data = OMEf[OMEf$Noise == "coherent",], method = "ML", +> ##D start = list(fixed=c(L75=c(48.7, -0.03), lsc=0.24)), verbose = TRUE))) +> ##D +> ##D try(summary(nlme(Resp ~ fp2(Loud, L75, lsc), +> ##D fixed = list(L75 ~ Age, lsc ~ 1), +> ##D random = L75 + lsc ~ 1 | UID, +> ##D data = OMEf[OMEf$Noise == "incoherent",], method = "ML", +> ##D start = list(fixed=c(L75=c(41.5, -0.1), lsc=0)), verbose = TRUE))) +> ## End(Not run) +> +> +> cleanEx() + +detaching ‘package:nlme’ + +> nameEx("Skye") +> ### * Skye +> +> flush(stderr()); flush(stdout()) +> +> ### Name: Skye +> ### Title: AFM Compositions of Aphyric Skye Lavas +> ### Aliases: Skye +> ### Keywords: datasets +> +> ### ** Examples +> +> # ternary() is from the on-line answers. +> ternary <- function(X, pch = par("pch"), lcex = 1, ++ add = FALSE, ord = 1:3, ...) ++ { ++ X <- as.matrix(X) ++ if(any(X < 0)) stop("X must be non-negative") ++ s <- drop(X %*% rep(1, ncol(X))) ++ if(any(s<=0)) stop("each row of X must have a positive sum") ++ if(max(abs(s-1)) > 1e-6) { ++ warning("row(s) of X will be rescaled") ++ X <- X / s ++ } ++ X <- X[, ord] ++ s3 <- sqrt(1/3) ++ if(!add) ++ { ++ oldpty <- par("pty") ++ on.exit(par(pty=oldpty)) ++ par(pty="s") ++ plot(c(-s3, s3), c(0.5-s3, 0.5+s3), type="n", axes=FALSE, ++ xlab="", ylab="") ++ polygon(c(0, -s3, s3), c(1, 0, 0), density=0) ++ lab <- NULL ++ if(!is.null(dn <- dimnames(X))) lab <- dn[[2]] ++ if(length(lab) < 3) lab <- as.character(1:3) ++ eps <- 0.05 * lcex ++ text(c(0, s3+eps*0.7, -s3-eps*0.7), ++ c(1+eps, -0.1*eps, -0.1*eps), lab, cex=lcex) ++ } ++ points((X[,2] - X[,3])*s3, X[,1], ...) ++ } +> +> ternary(Skye/100, ord=c(1,3,2)) +> +> +> +> graphics::par(get("par.postscript", pos = 'CheckExEnv')) +> cleanEx() +> nameEx("addterm") +> ### * addterm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: addterm +> ### Title: Try All One-Term Additions to a Model +> ### Aliases: addterm addterm.default addterm.glm addterm.lm +> ### Keywords: models +> +> ### ** Examples +> +> quine.hi <- aov(log(Days + 2.5) ~ .^4, quine) +> quine.lo <- aov(log(Days+2.5) ~ 1, quine) +> addterm(quine.lo, quine.hi, test="F") +Single term additions + +Model: +log(Days + 2.5) ~ 1 + Df Sum of Sq RSS AIC F Value Pr(F) + 106.787 -43.664 +Eth 1 10.6820 96.105 -57.052 16.0055 0.0001006 *** +Sex 1 0.5969 106.190 -42.483 0.8094 0.3698057 +Age 3 4.7469 102.040 -44.303 2.2019 0.0904804 . +Lrn 1 0.0043 106.783 -41.670 0.0058 0.9392083 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family=poisson, ++ data=housing) +> addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test="Chisq") +Single term additions + +Model: +Freq ~ Infl * Type * Cont + Sat + Df Deviance AIC LRT Pr(Chi) + 217.46 610.43 +Infl:Sat 4 111.08 512.05 106.371 < 2.2e-16 *** +Type:Sat 6 156.79 561.76 60.669 3.292e-11 *** +Cont:Sat 2 212.33 609.30 5.126 0.07708 . +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) +> addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test = "Chisq") +Single term additions + +Model: +Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont + + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont + Df Deviance AIC LRT Pr(Chi) + 38.662 455.63 +Infl:Type:Sat 12 16.107 457.08 22.5550 0.03175 * +Infl:Cont:Sat 4 37.472 462.44 1.1901 0.87973 +Type:Cont:Sat 6 28.256 457.23 10.4064 0.10855 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> +> +> cleanEx() +> nameEx("anova.negbin") +> ### * anova.negbin +> +> flush(stderr()); flush(stdout()) +> +> ### Name: anova.negbin +> ### Title: Likelihood Ratio Tests for Negative Binomial GLMs +> ### Aliases: anova.negbin +> ### Keywords: regression +> +> ### ** Examples +> +> m1 <- glm.nb(Days ~ Eth*Age*Lrn*Sex, quine, link = log) +> m2 <- update(m1, . ~ . - Eth:Age:Lrn:Sex) +> anova(m2, m1) +Likelihood ratio tests of Negative Binomial Models + +Response: Days + Model +1 Eth + Age + Lrn + Sex + Eth:Age + Eth:Lrn + Age:Lrn + Eth:Sex + Age:Sex + Lrn:Sex + Eth:Age:Lrn + Eth:Age:Sex + Eth:Lrn:Sex + Age:Lrn:Sex +2 Eth * Age * Lrn * Sex + theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi) +1 1.90799 120 -1040.728 +2 1.92836 118 -1039.324 1 vs 2 2 1.403843 0.4956319 +> anova(m2) +Warning in anova.negbin(m2) : tests made without re-estimating 'theta' +Calls: anova -> anova.negbin +Analysis of Deviance Table + +Model: Negative Binomial(1.908), link: log + +Response: Days + +Terms added sequentially (first to last) + + + Df Deviance Resid. Df Resid. Dev Pr(>Chi) +NULL 145 270.03 +Eth 1 19.0989 144 250.93 1.241e-05 *** +Age 3 16.3483 141 234.58 0.000962 *** +Lrn 1 3.5449 140 231.04 0.059730 . +Sex 1 0.3989 139 230.64 0.527666 +Eth:Age 3 14.6030 136 216.03 0.002189 ** +Eth:Lrn 1 0.0447 135 215.99 0.832601 +Age:Lrn 2 1.7482 133 214.24 0.417240 +Eth:Sex 1 1.1470 132 213.09 0.284183 +Age:Sex 3 21.9746 129 191.12 6.603e-05 *** +Lrn:Sex 1 0.0277 128 191.09 0.867712 +Eth:Age:Lrn 2 9.0099 126 182.08 0.011054 * +Eth:Age:Sex 3 4.8218 123 177.26 0.185319 +Eth:Lrn:Sex 1 3.3160 122 173.94 0.068608 . +Age:Lrn:Sex 2 6.3941 120 167.55 0.040882 * +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> +> +> cleanEx() +> nameEx("area") +> ### * area +> +> flush(stderr()); flush(stdout()) +> +> ### Name: area +> ### Title: Adaptive Numerical Integration +> ### Aliases: area +> ### Keywords: nonlinear +> +> ### ** Examples +> +> area(sin, 0, pi) # integrate the sin function from 0 to pi. +[1] 2 +> +> +> +> cleanEx() +> nameEx("bacteria") +> ### * bacteria +> +> flush(stderr()); flush(stdout()) +> +> ### Name: bacteria +> ### Title: Presence of Bacteria after Drug Treatments +> ### Aliases: bacteria +> ### Keywords: datasets +> +> ### ** Examples +> +> contrasts(bacteria$trt) <- structure(contr.sdif(3), ++ dimnames = list(NULL, c("drug", "encourage"))) +> ## fixed effects analyses +> ## IGNORE_RDIFF_BEGIN +> summary(glm(y ~ trt * week, binomial, data = bacteria)) + +Call: +glm(formula = y ~ trt * week, family = binomial, data = bacteria) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 1.97548 0.30053 6.573 4.92e-11 *** +trtdrug -0.99848 0.69490 -1.437 0.15075 +trtencourage 0.83865 0.73482 1.141 0.25374 +week -0.11814 0.04460 -2.649 0.00807 ** +trtdrug:week -0.01722 0.10570 -0.163 0.87061 +trtencourage:week -0.07043 0.10964 -0.642 0.52060 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for binomial family taken to be 1) + + Null deviance: 217.38 on 219 degrees of freedom +Residual deviance: 203.12 on 214 degrees of freedom +AIC: 215.12 + +Number of Fisher Scoring iterations: 4 + +> summary(glm(y ~ trt + week, binomial, data = bacteria)) + +Call: +glm(formula = y ~ trt + week, family = binomial, data = bacteria) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 1.96018 0.29705 6.599 4.15e-11 *** +trtdrug -1.10667 0.42519 -2.603 0.00925 ** +trtencourage 0.45502 0.42766 1.064 0.28735 +week -0.11577 0.04414 -2.623 0.00872 ** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for binomial family taken to be 1) + + Null deviance: 217.38 on 219 degrees of freedom +Residual deviance: 203.81 on 216 degrees of freedom +AIC: 211.81 + +Number of Fisher Scoring iterations: 4 + +> summary(glm(y ~ trt + I(week > 2), binomial, data = bacteria)) + +Call: +glm(formula = y ~ trt + I(week > 2), family = binomial, data = bacteria) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 2.2479 0.3560 6.315 2.71e-10 *** +trtdrug -1.1187 0.4288 -2.609 0.00909 ** +trtencourage 0.4815 0.4330 1.112 0.26614 +I(week > 2)TRUE -1.2949 0.4104 -3.155 0.00160 ** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for binomial family taken to be 1) + + Null deviance: 217.38 on 219 degrees of freedom +Residual deviance: 199.18 on 216 degrees of freedom +AIC: 207.18 + +Number of Fisher Scoring iterations: 5 + +> ## IGNORE_RDIFF_END +> +> # conditional random-effects analysis +> library(survival) +> bacteria$Time <- rep(1, nrow(bacteria)) +> coxph(Surv(Time, unclass(y)) ~ week + strata(ID), ++ data = bacteria, method = "exact") +Call: +coxph(formula = Surv(Time, unclass(y)) ~ week + strata(ID), data = bacteria, + method = "exact") + + coef exp(coef) se(coef) z p +week -0.16256 0.84996 0.05472 -2.971 0.00297 + +Likelihood ratio test=9.85 on 1 df, p=0.001696 +n= 220, number of events= 177 +> coxph(Surv(Time, unclass(y)) ~ factor(week) + strata(ID), ++ data = bacteria, method = "exact") +Call: +coxph(formula = Surv(Time, unclass(y)) ~ factor(week) + strata(ID), + data = bacteria, method = "exact") + + coef exp(coef) se(coef) z p +factor(week)2 0.1983 1.2193 0.7241 0.274 0.7842 +factor(week)4 -1.4206 0.2416 0.6665 -2.131 0.0331 +factor(week)6 -1.6615 0.1899 0.6825 -2.434 0.0149 +factor(week)11 -1.6752 0.1873 0.6780 -2.471 0.0135 + +Likelihood ratio test=15.45 on 4 df, p=0.003854 +n= 220, number of events= 177 +> coxph(Surv(Time, unclass(y)) ~ I(week > 2) + strata(ID), ++ data = bacteria, method = "exact") +Call: +coxph(formula = Surv(Time, unclass(y)) ~ I(week > 2) + strata(ID), + data = bacteria, method = "exact") + + coef exp(coef) se(coef) z p +I(week > 2)TRUE -1.6701 0.1882 0.4817 -3.467 0.000527 + +Likelihood ratio test=15.15 on 1 df, p=9.927e-05 +n= 220, number of events= 177 +> +> # PQL glmm analysis +> library(nlme) +> ## IGNORE_RDIFF_BEGIN +> summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID, ++ family = binomial, data = bacteria)) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +iteration 6 +Linear mixed-effects model fit by maximum likelihood + Data: bacteria + AIC BIC logLik + NA NA NA + +Random effects: + Formula: ~1 | ID + (Intercept) Residual +StdDev: 1.410637 0.7800511 + +Variance function: + Structure: fixed weights + Formula: ~invwt +Fixed effects: y ~ trt + I(week > 2) + Value Std.Error DF t-value p-value +(Intercept) 2.7447864 0.3784193 169 7.253294 0.0000 +trtdrug -1.2473553 0.6440635 47 -1.936696 0.0588 +trtencourage 0.4930279 0.6699339 47 0.735935 0.4654 +I(week > 2)TRUE -1.6072570 0.3583379 169 -4.485311 0.0000 + Correlation: + (Intr) trtdrg trtncr +trtdrug 0.009 +trtencourage 0.036 -0.518 +I(week > 2)TRUE -0.710 0.047 -0.046 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-5.1985361 0.1572336 0.3513075 0.4949482 1.7448845 + +Number of Observations: 220 +Number of Groups: 50 +> ## IGNORE_RDIFF_END +> +> +> +> cleanEx() + +detaching ‘package:nlme’, ‘package:survival’ + +> nameEx("bandwidth.nrd") +> ### * bandwidth.nrd +> +> flush(stderr()); flush(stdout()) +> +> ### Name: bandwidth.nrd +> ### Title: Bandwidth for density() via Normal Reference Distribution +> ### Aliases: bandwidth.nrd +> ### Keywords: dplot +> +> ### ** Examples +> +> # The function is currently defined as +> function(x) ++ { ++ r <- quantile(x, c(0.25, 0.75)) ++ h <- (r[2] - r[1])/1.34 ++ 4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5) ++ } +function (x) +{ + r <- quantile(x, c(0.25, 0.75)) + h <- (r[2] - r[1])/1.34 + 4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5) +} +> +> +> +> cleanEx() +> nameEx("bcv") +> ### * bcv +> +> flush(stderr()); flush(stdout()) +> +> ### Name: bcv +> ### Title: Biased Cross-Validation for Bandwidth Selection +> ### Aliases: bcv +> ### Keywords: dplot +> +> ### ** Examples +> +> bcv(geyser$duration) +[1] 0.8940809 +> +> +> +> cleanEx() +> nameEx("beav1") +> ### * beav1 +> +> flush(stderr()); flush(stdout()) +> +> ### Name: beav1 +> ### Title: Body Temperature Series of Beaver 1 +> ### Aliases: beav1 +> ### Keywords: datasets +> +> ### ** Examples +> +> beav1 <- within(beav1, ++ hours <- 24*(day-346) + trunc(time/100) + (time%%100)/60) +> plot(beav1$hours, beav1$temp, type="l", xlab="time", ++ ylab="temperature", main="Beaver 1") +> usr <- par("usr"); usr[3:4] <- c(-0.2, 8); par(usr=usr) +> lines(beav1$hours, beav1$activ, type="s", lty=2) +> temp <- ts(c(beav1$temp[1:82], NA, beav1$temp[83:114]), ++ start = 9.5, frequency = 6) +> activ <- ts(c(beav1$activ[1:82], NA, beav1$activ[83:114]), ++ start = 9.5, frequency = 6) +> +> acf(temp[1:53]) +> acf(temp[1:53], type = "partial") +> ar(temp[1:53]) + +Call: +ar(x = temp[1:53]) + +Coefficients: + 1 +0.8222 + +Order selected 1 sigma^2 estimated as 0.01011 +> act <- c(rep(0, 10), activ) +> X <- cbind(1, act = act[11:125], act1 = act[10:124], ++ act2 = act[9:123], act3 = act[8:122]) +> alpha <- 0.80 +> stemp <- as.vector(temp - alpha*lag(temp, -1)) +> sX <- X[-1, ] - alpha * X[-115,] +> beav1.ls <- lm(stemp ~ -1 + sX, na.action = na.omit) +> summary(beav1.ls, correlation = FALSE) + +Call: +lm(formula = stemp ~ -1 + sX, na.action = na.omit) + +Residuals: + Min 1Q Median 3Q Max +-0.21317 -0.04317 0.00683 0.05483 0.37683 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +sX 36.85587 0.03922 939.833 < 2e-16 *** +sXact 0.25400 0.03930 6.464 3.37e-09 *** +sXact1 0.17096 0.05100 3.352 0.00112 ** +sXact2 0.16202 0.05147 3.148 0.00215 ** +sXact3 0.10548 0.04310 2.448 0.01605 * +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 0.08096 on 104 degrees of freedom + (5 observations deleted due to missingness) +Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999 +F-statistic: 1.81e+05 on 5 and 104 DF, p-value: < 2.2e-16 + +> rm(temp, activ) +> +> +> +> graphics::par(get("par.postscript", pos = 'CheckExEnv')) +> cleanEx() +> nameEx("beav2") +> ### * beav2 +> +> flush(stderr()); flush(stdout()) +> +> ### Name: beav2 +> ### Title: Body Temperature Series of Beaver 2 +> ### Aliases: beav2 +> ### Keywords: datasets +> +> ### ** Examples +> +> attach(beav2) +> beav2$hours <- 24*(day-307) + trunc(time/100) + (time%%100)/60 +> plot(beav2$hours, beav2$temp, type = "l", xlab = "time", ++ ylab = "temperature", main = "Beaver 2") +> usr <- par("usr"); usr[3:4] <- c(-0.2, 8); par(usr = usr) +> lines(beav2$hours, beav2$activ, type = "s", lty = 2) +> +> temp <- ts(temp, start = 8+2/3, frequency = 6) +> activ <- ts(activ, start = 8+2/3, frequency = 6) +> acf(temp[activ == 0]); acf(temp[activ == 1]) # also look at PACFs +> ar(temp[activ == 0]); ar(temp[activ == 1]) + +Call: +ar(x = temp[activ == 0]) + +Coefficients: + 1 +0.7392 + +Order selected 1 sigma^2 estimated as 0.02011 + +Call: +ar(x = temp[activ == 1]) + +Coefficients: + 1 +0.7894 + +Order selected 1 sigma^2 estimated as 0.01792 +> +> arima(temp, order = c(1,0,0), xreg = activ) + +Call: +arima(x = temp, order = c(1, 0, 0), xreg = activ) + +Coefficients: + ar1 intercept activ + 0.8733 37.1920 0.6139 +s.e. 0.0684 0.1187 0.1381 + +sigma^2 estimated as 0.01518: log likelihood = 66.78, aic = -125.55 +> dreg <- cbind(sin = sin(2*pi*beav2$hours/24), cos = cos(2*pi*beav2$hours/24)) +> arima(temp, order = c(1,0,0), xreg = cbind(active=activ, dreg)) + +Call: +arima(x = temp, order = c(1, 0, 0), xreg = cbind(active = activ, dreg)) + +Coefficients: + ar1 intercept active dreg.sin dreg.cos + 0.7905 37.1674 0.5322 -0.282 0.1201 +s.e. 0.0681 0.0939 0.1282 0.105 0.0997 + +sigma^2 estimated as 0.01434: log likelihood = 69.83, aic = -127.67 +> +> ## IGNORE_RDIFF_BEGIN +> library(nlme) # for gls and corAR1 +> beav2.gls <- gls(temp ~ activ, data = beav2, correlation = corAR1(0.8), ++ method = "ML") +> summary(beav2.gls) +Generalized least squares fit by maximum likelihood + Model: temp ~ activ + Data: beav2 + AIC BIC logLik + -125.5505 -115.1298 66.77523 + +Correlation Structure: AR(1) + Formula: ~1 + Parameter estimate(s): + Phi +0.8731771 + +Coefficients: + Value Std.Error t-value p-value +(Intercept) 37.19195 0.1131328 328.7460 0 +activ 0.61418 0.1087286 5.6487 0 + + Correlation: + (Intr) +activ -0.582 + +Standardized residuals: + Min Q1 Med Q3 Max +-2.42080776 -0.61510519 -0.03573836 0.81641138 2.15153495 + +Residual standard error: 0.2527856 +Degrees of freedom: 100 total; 98 residual +> summary(update(beav2.gls, subset = 6:100)) +Generalized least squares fit by maximum likelihood + Model: temp ~ activ + Data: beav2 + Subset: 6:100 + AIC BIC logLik + -124.981 -114.7654 66.49048 + +Correlation Structure: AR(1) + Formula: ~1 + Parameter estimate(s): + Phi +0.8380448 + +Coefficients: + Value Std.Error t-value p-value +(Intercept) 37.25001 0.09634047 386.6496 0 +activ 0.60277 0.09931904 6.0690 0 + + Correlation: + (Intr) +activ -0.657 + +Standardized residuals: + Min Q1 Med Q3 Max +-2.0231494 -0.8910348 -0.1497564 0.7640939 2.2719468 + +Residual standard error: 0.2188542 +Degrees of freedom: 95 total; 93 residual +> detach("beav2"); rm(temp, activ) +> ## IGNORE_RDIFF_END +> +> +> +> graphics::par(get("par.postscript", pos = 'CheckExEnv')) +> cleanEx() + +detaching ‘package:nlme’ + +> nameEx("birthwt") +> ### * birthwt +> +> flush(stderr()); flush(stdout()) +> +> ### Name: birthwt +> ### Title: Risk Factors Associated with Low Infant Birth Weight +> ### Aliases: birthwt +> ### Keywords: datasets +> +> ### ** Examples +> +> bwt <- with(birthwt, { ++ race <- factor(race, labels = c("white", "black", "other")) ++ ptd <- factor(ptl > 0) ++ ftv <- factor(ftv) ++ levels(ftv)[-(1:2)] <- "2+" ++ data.frame(low = factor(low), age, lwt, race, smoke = (smoke > 0), ++ ptd, ht = (ht > 0), ui = (ui > 0), ftv) ++ }) +> options(contrasts = c("contr.treatment", "contr.poly")) +> glm(low ~ ., binomial, bwt) + +Call: glm(formula = low ~ ., family = binomial, data = bwt) + +Coefficients: +(Intercept) age lwt raceblack raceother smokeTRUE + 0.82302 -0.03723 -0.01565 1.19241 0.74068 0.75553 + ptdTRUE htTRUE uiTRUE ftv1 ftv2+ + 1.34376 1.91317 0.68020 -0.43638 0.17901 + +Degrees of Freedom: 188 Total (i.e. Null); 178 Residual +Null Deviance: 234.7 +Residual Deviance: 195.5 AIC: 217.5 +> +> +> +> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly")) +> cleanEx() +> nameEx("boxcox") +> ### * boxcox +> +> flush(stderr()); flush(stdout()) +> +> ### Name: boxcox +> ### Title: Box-Cox Transformations for Linear Models +> ### Aliases: boxcox boxcox.default boxcox.formula boxcox.lm +> ### Keywords: regression models hplot +> +> ### ** Examples +> +> boxcox(Volume ~ log(Height) + log(Girth), data = trees, ++ lambda = seq(-0.25, 0.25, length.out = 10)) +> +> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, ++ lambda = seq(-0.05, 0.45, length.out = 20)) +> +> +> +> cleanEx() +> nameEx("caith") +> ### * caith +> +> flush(stderr()); flush(stdout()) +> +> ### Name: caith +> ### Title: Colours of Eyes and Hair of People in Caithness +> ### Aliases: caith +> ### Keywords: datasets +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> ## The signs can vary by platform +> corresp(caith) +First canonical correlation(s): 0.4463684 + + Row scores: + blue light medium dark + 0.89679252 0.98731818 -0.07530627 -1.57434710 + + Column scores: + fair red medium dark black + 1.21871379 0.52257500 0.09414671 -1.31888486 -2.45176017 +> ## IGNORE_RDIFF_END +> dimnames(caith)[[2]] <- c("F", "R", "M", "D", "B") +> par(mfcol=c(1,3)) +> plot(corresp(caith, nf=2)); title("symmetric") +> plot(corresp(caith, nf=2), type="rows"); title("rows") +> plot(corresp(caith, nf=2), type="col"); title("columns") +> par(mfrow=c(1,1)) +> +> +> +> graphics::par(get("par.postscript", pos = 'CheckExEnv')) +> cleanEx() +> nameEx("cement") +> ### * cement +> +> flush(stderr()); flush(stdout()) +> +> ### Name: cement +> ### Title: Heat Evolved by Setting Cements +> ### Aliases: cement +> ### Keywords: datasets +> +> ### ** Examples +> +> lm(y ~ x1 + x2 + x3 + x4, cement) + +Call: +lm(formula = y ~ x1 + x2 + x3 + x4, data = cement) + +Coefficients: +(Intercept) x1 x2 x3 x4 + 62.4054 1.5511 0.5102 0.1019 -0.1441 + +> +> +> +> cleanEx() +> nameEx("confint") +> ### * confint +> +> flush(stderr()); flush(stdout()) +> +> ### Name: confint-MASS +> ### Title: Confidence Intervals for Model Parameters +> ### Aliases: confint.glm confint.nls confint.profile.glm +> ### confint.profile.nls +> ### Keywords: models +> +> ### ** Examples +> +> expn1 <- deriv(y ~ b0 + b1 * 2^(-x/th), c("b0", "b1", "th"), ++ function(b0, b1, th, x) {}) +> +> wtloss.gr <- nls(Weight ~ expn1(b0, b1, th, Days), ++ data = wtloss, start = c(b0=90, b1=95, th=120)) +> +> expn2 <- deriv(~b0 + b1*((w0 - b0)/b1)^(x/d0), ++ c("b0","b1","d0"), function(b0, b1, d0, x, w0) {}) +> +> wtloss.init <- function(obj, w0) { ++ p <- coef(obj) ++ d0 <- - log((w0 - p["b0"])/p["b1"])/log(2) * p["th"] ++ c(p[c("b0", "b1")], d0 = as.vector(d0)) ++ } +> +> out <- NULL +> w0s <- c(110, 100, 90) +> for(w0 in w0s) { ++ fm <- nls(Weight ~ expn2(b0, b1, d0, Days, w0), ++ wtloss, start = wtloss.init(wtloss.gr, w0)) ++ out <- rbind(out, c(coef(fm)["d0"], confint(fm, "d0"))) ++ } +Waiting for profiling to be done... +Waiting for profiling to be done... +Waiting for profiling to be done... +> dimnames(out) <- list(paste(w0s, "kg:"), c("d0", "low", "high")) +> out + d0 low high +110 kg: 261.5132 256.2303 267.5009 +100 kg: 349.4979 334.7293 368.0151 +90 kg: 507.0941 457.2637 594.8745 +> +> ldose <- rep(0:5, 2) +> numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16) +> sex <- factor(rep(c("M", "F"), c(6, 6))) +> SF <- cbind(numdead, numalive = 20 - numdead) +> budworm.lg0 <- glm(SF ~ sex + ldose - 1, family = binomial) +> confint(budworm.lg0) +Waiting for profiling to be done... + 2.5 % 97.5 % +sexF -4.4581438 -2.613610 +sexM -3.1728745 -1.655117 +ldose 0.8228708 1.339058 +> confint(budworm.lg0, "ldose") +Waiting for profiling to be done... + 2.5 % 97.5 % +0.8228708 1.3390581 +> +> +> +> cleanEx() +> nameEx("contr.sdif") +> ### * contr.sdif +> +> flush(stderr()); flush(stdout()) +> +> ### Name: contr.sdif +> ### Title: Successive Differences Contrast Coding +> ### Aliases: contr.sdif +> ### Keywords: models +> +> ### ** Examples +> +> (A <- contr.sdif(6)) + 2-1 3-2 4-3 5-4 6-5 +1 -0.8333333 -0.6666667 -0.5 -0.3333333 -0.1666667 +2 0.1666667 -0.6666667 -0.5 -0.3333333 -0.1666667 +3 0.1666667 0.3333333 -0.5 -0.3333333 -0.1666667 +4 0.1666667 0.3333333 0.5 -0.3333333 -0.1666667 +5 0.1666667 0.3333333 0.5 0.6666667 -0.1666667 +6 0.1666667 0.3333333 0.5 0.6666667 0.8333333 +> zapsmall(ginv(A)) + [,1] [,2] [,3] [,4] [,5] [,6] +[1,] -1 1 0 0 0 0 +[2,] 0 -1 1 0 0 0 +[3,] 0 0 -1 1 0 0 +[4,] 0 0 0 -1 1 0 +[5,] 0 0 0 0 -1 1 +> +> +> +> cleanEx() +> nameEx("corresp") +> ### * corresp +> +> flush(stderr()); flush(stdout()) +> +> ### Name: corresp +> ### Title: Simple Correspondence Analysis +> ### Aliases: corresp corresp.xtabs corresp.data.frame corresp.default +> ### corresp.factor corresp.formula corresp.matrix +> ### Keywords: category multivariate +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> ## The signs can vary by platform +> (ct <- corresp(~ Age + Eth, data = quine)) +First canonical correlation(s): 0.05317534 + + Age scores: + F0 F1 F2 F3 +-0.3344445 1.4246090 -1.0320002 -0.4612728 + + Eth scores: + A N +-1.0563816 0.9466276 +> plot(ct) +> +> corresp(caith) +First canonical correlation(s): 0.4463684 + + Row scores: + blue light medium dark + 0.89679252 0.98731818 -0.07530627 -1.57434710 + + Column scores: + fair red medium dark black + 1.21871379 0.52257500 0.09414671 -1.31888486 -2.45176017 +> biplot(corresp(caith, nf = 2)) +> ## IGNORE_RDIFF_END +> +> +> +> cleanEx() +> nameEx("cov.rob") +> ### * cov.rob +> +> flush(stderr()); flush(stdout()) +> +> ### Name: cov.rob +> ### Title: Resistant Estimation of Multivariate Location and Scatter +> ### Aliases: cov.rob cov.mve cov.mcd +> ### Keywords: robust multivariate +> +> ### ** Examples +> +> set.seed(123) +> cov.rob(stackloss) +$center + Air.Flow Water.Temp Acid.Conc. stack.loss + 56.3750 20.0000 85.4375 13.0625 + +$cov + Air.Flow Water.Temp Acid.Conc. stack.loss +Air.Flow 23.050000 6.666667 16.625000 19.308333 +Water.Temp 6.666667 5.733333 5.333333 7.733333 +Acid.Conc. 16.625000 5.333333 34.395833 13.837500 +stack.loss 19.308333 7.733333 13.837500 18.462500 + +$msg +[1] "20 singular samples of size 5 out of 2500" + +$crit +[1] 19.89056 + +$best + [1] 5 6 7 8 9 10 11 12 15 16 18 19 20 + +$n.obs +[1] 21 + +> cov.rob(stack.x, method = "mcd", nsamp = "exact") +$center + Air.Flow Water.Temp Acid.Conc. + 56.70588 20.23529 85.52941 + +$cov + Air.Flow Water.Temp Acid.Conc. +Air.Flow 23.470588 7.573529 16.102941 +Water.Temp 7.573529 6.316176 5.367647 +Acid.Conc. 16.102941 5.367647 32.389706 + +$msg +[1] "266 singular samples of size 4 out of 5985" + +$crit +[1] 5.472581 + +$best + [1] 4 5 6 7 8 9 10 11 12 13 14 20 + +$n.obs +[1] 21 + +> +> +> +> cleanEx() +> nameEx("cov.trob") +> ### * cov.trob +> +> flush(stderr()); flush(stdout()) +> +> ### Name: cov.trob +> ### Title: Covariance Estimation for Multivariate t Distribution +> ### Aliases: cov.trob +> ### Keywords: multivariate +> +> ### ** Examples +> +> cov.trob(stackloss) +$cov + Air.Flow Water.Temp Acid.Conc. stack.loss +Air.Flow 60.47035 17.027203 18.554452 62.28032 +Water.Temp 17.02720 8.085857 5.604132 20.50469 +Acid.Conc. 18.55445 5.604132 24.404633 16.91085 +stack.loss 62.28032 20.504687 16.910855 72.80743 + +$center + Air.Flow Water.Temp Acid.Conc. stack.loss + 58.96905 20.79263 86.05588 16.09028 + +$n.obs +[1] 21 + +$call +cov.trob(x = stackloss) + +$iter +[1] 5 + +> +> +> +> cleanEx() +> nameEx("denumerate") +> ### * denumerate +> +> flush(stderr()); flush(stdout()) +> +> ### Name: denumerate +> ### Title: Transform an Allowable Formula for 'loglm' into one for 'terms' +> ### Aliases: denumerate denumerate.formula +> ### Keywords: models +> +> ### ** Examples +> +> denumerate(~(1+2+3)^3 + a/b) +~(.v1 + .v2 + .v3)^3 + a/b +> ## which gives ~ (.v1 + .v2 + .v3)^3 + a/b +> +> +> +> cleanEx() +> nameEx("dose.p") +> ### * dose.p +> +> flush(stderr()); flush(stdout()) +> +> ### Name: dose.p +> ### Title: Predict Doses for Binomial Assay model +> ### Aliases: dose.p print.glm.dose +> ### Keywords: regression models +> +> ### ** Examples +> +> ldose <- rep(0:5, 2) +> numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16) +> sex <- factor(rep(c("M", "F"), c(6, 6))) +> SF <- cbind(numdead, numalive = 20 - numdead) +> budworm.lg0 <- glm(SF ~ sex + ldose - 1, family = binomial) +> +> dose.p(budworm.lg0, cf = c(1,3), p = 1:3/4) + Dose SE +p = 0.25: 2.231265 0.2499089 +p = 0.50: 3.263587 0.2297539 +p = 0.75: 4.295910 0.2746874 +> dose.p(update(budworm.lg0, family = binomial(link=probit)), ++ cf = c(1,3), p = 1:3/4) + Dose SE +p = 0.25: 2.191229 0.2384478 +p = 0.50: 3.257703 0.2240685 +p = 0.75: 4.324177 0.2668745 +> +> +> +> cleanEx() +> nameEx("dropterm") +> ### * dropterm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: dropterm +> ### Title: Try All One-Term Deletions from a Model +> ### Aliases: dropterm dropterm.default dropterm.glm dropterm.lm +> ### Keywords: models +> +> ### ** Examples +> +> quine.hi <- aov(log(Days + 2.5) ~ .^4, quine) +> quine.nxt <- update(quine.hi, . ~ . - Eth:Sex:Age:Lrn) +> dropterm(quine.nxt, test= "F") +Single term deletions + +Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Age + Eth:Sex:Lrn + + Eth:Age:Lrn + Sex:Age:Lrn + Df Sum of Sq RSS AIC F Value Pr(F) + 64.099 -68.184 +Eth:Sex:Age 3 0.97387 65.073 -71.982 0.60773 0.61125 +Eth:Sex:Lrn 1 1.57879 65.678 -66.631 2.95567 0.08816 . +Eth:Age:Lrn 2 2.12841 66.227 -67.415 1.99230 0.14087 +Sex:Age:Lrn 2 1.46623 65.565 -68.882 1.37247 0.25743 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> quine.stp <- stepAIC(quine.nxt, ++ scope = list(upper = ~Eth*Sex*Age*Lrn, lower = ~1), ++ trace = FALSE) +> dropterm(quine.stp, test = "F") +Single term deletions + +Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + Df Sum of Sq RSS AIC F Value Pr(F) + 66.600 -72.597 +Sex:Age 3 10.7959 77.396 -56.663 6.7542 0.0002933 *** +Eth:Sex:Lrn 1 3.0325 69.632 -68.096 5.6916 0.0185476 * +Eth:Age:Lrn 2 2.0960 68.696 -72.072 1.9670 0.1441822 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> quine.3 <- update(quine.stp, . ~ . - Eth:Age:Lrn) +> dropterm(quine.3, test = "F") +Single term deletions + +Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Df Sum of Sq RSS AIC F Value Pr(F) + 68.696 -72.072 +Eth:Age 3 3.0312 71.727 -71.768 1.8679 0.1383323 +Sex:Age 3 11.4272 80.123 -55.607 7.0419 0.0002037 *** +Age:Lrn 2 2.8149 71.511 -70.209 2.6020 0.0780701 . +Eth:Sex:Lrn 1 4.6956 73.391 -64.419 8.6809 0.0038268 ** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> quine.4 <- update(quine.3, . ~ . - Eth:Age) +> dropterm(quine.4, test = "F") +Single term deletions + +Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Lrn + + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Df Sum of Sq RSS AIC F Value Pr(F) + 71.727 -71.768 +Sex:Age 3 11.5656 83.292 -55.942 6.9873 0.0002147 *** +Age:Lrn 2 2.9118 74.639 -69.959 2.6387 0.0752793 . +Eth:Sex:Lrn 1 6.8181 78.545 -60.511 12.3574 0.0006052 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> quine.5 <- update(quine.4, . ~ . - Age:Lrn) +> dropterm(quine.5, test = "F") +Single term deletions + +Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Lrn + + Sex:Age + Sex:Lrn + Eth:Sex:Lrn + Df Sum of Sq RSS AIC F Value Pr(F) + 74.639 -69.959 +Sex:Age 3 9.9002 84.539 -57.774 5.8362 0.0008944 *** +Eth:Sex:Lrn 1 6.2988 80.937 -60.130 11.1396 0.0010982 ** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family=poisson, ++ data = housing) +> house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) +> dropterm(house.glm1, test = "Chisq") +Single term deletions + +Model: +Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont + + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont + Df Deviance AIC LRT Pr(Chi) + 38.662 455.63 +Infl:Sat 4 147.780 556.75 109.117 < 2.2e-16 *** +Type:Sat 6 100.889 505.86 62.227 1.586e-11 *** +Cont:Sat 2 54.722 467.69 16.060 0.0003256 *** +Infl:Type:Cont 6 43.952 448.92 5.290 0.5072454 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> +> +> cleanEx() +> nameEx("eagles") +> ### * eagles +> +> flush(stderr()); flush(stdout()) +> +> ### Name: eagles +> ### Title: Foraging Ecology of Bald Eagles +> ### Aliases: eagles +> ### Keywords: datasets +> +> ### ** Examples +> +> eagles.glm <- glm(cbind(y, n - y) ~ P*A + V, data = eagles, ++ family = binomial) +> dropterm(eagles.glm) +Single term deletions + +Model: +cbind(y, n - y) ~ P * A + V + Df Deviance AIC + 0.333 23.073 +V 1 53.737 74.478 +P:A 1 6.956 27.696 +> prof <- profile(eagles.glm) +> plot(prof) +> pairs(prof) +> +> +> +> cleanEx() +> nameEx("epil") +> ### * epil +> +> flush(stderr()); flush(stdout()) +> +> ### Name: epil +> ### Title: Seizure Counts for Epileptics +> ### Aliases: epil +> ### Keywords: datasets +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> summary(glm(y ~ lbase*trt + lage + V4, family = poisson, ++ data = epil), correlation = FALSE) + +Call: +glm(formula = y ~ lbase * trt + lage + V4, family = poisson, + data = epil) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 1.89791 0.04260 44.552 < 2e-16 *** +lbase 0.94862 0.04360 21.759 < 2e-16 *** +trtprogabide -0.34588 0.06100 -5.670 1.42e-08 *** +lage 0.88760 0.11650 7.619 2.56e-14 *** +V4 -0.15977 0.05458 -2.927 0.00342 ** +lbase:trtprogabide 0.56154 0.06352 8.841 < 2e-16 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for poisson family taken to be 1) + + Null deviance: 2517.83 on 235 degrees of freedom +Residual deviance: 869.07 on 230 degrees of freedom +AIC: 1647 + +Number of Fisher Scoring iterations: 5 + +> ## IGNORE_RDIFF_END +> epil2 <- epil[epil$period == 1, ] +> epil2["period"] <- rep(0, 59); epil2["y"] <- epil2["base"] +> epil["time"] <- 1; epil2["time"] <- 4 +> epil2 <- rbind(epil, epil2) +> epil2$pred <- unclass(epil2$trt) * (epil2$period > 0) +> epil2$subject <- factor(epil2$subject) +> epil3 <- aggregate(epil2, list(epil2$subject, epil2$period > 0), ++ function(x) if(is.numeric(x)) sum(x) else x[1]) +> epil3$pred <- factor(epil3$pred, ++ labels = c("base", "placebo", "drug")) +> +> contrasts(epil3$pred) <- structure(contr.sdif(3), ++ dimnames = list(NULL, c("placebo-base", "drug-placebo"))) +> ## IGNORE_RDIFF_BEGIN +> summary(glm(y ~ pred + factor(subject) + offset(log(time)), ++ family = poisson, data = epil3), correlation = FALSE) + +Call: +glm(formula = y ~ pred + factor(subject) + offset(log(time)), + family = poisson, data = epil3) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 1.122e+00 2.008e-01 5.590 2.28e-08 *** +predplacebo-base 1.087e-01 4.691e-02 2.318 0.020474 * +preddrug-placebo -1.016e-01 6.507e-02 -1.561 0.118431 +factor(subject)2 -9.105e-16 2.828e-01 0.000 1.000000 +factor(subject)3 -3.857e-01 3.144e-01 -1.227 0.219894 +factor(subject)4 -1.744e-01 2.960e-01 -0.589 0.555847 +factor(subject)5 1.577e+00 2.197e-01 7.178 7.08e-13 *** +factor(subject)6 6.729e-01 2.458e-01 2.738 0.006182 ** +factor(subject)7 -4.082e-02 2.858e-01 -0.143 0.886411 +factor(subject)8 1.758e+00 2.166e-01 8.117 4.77e-16 *** +factor(subject)9 5.878e-01 2.494e-01 2.356 0.018454 * +factor(subject)10 5.423e-01 2.515e-01 2.156 0.031060 * +factor(subject)11 1.552e+00 2.202e-01 7.048 1.81e-12 *** +factor(subject)12 9.243e-01 2.364e-01 3.910 9.22e-05 *** +factor(subject)13 3.075e-01 2.635e-01 1.167 0.243171 +factor(subject)14 1.212e+00 2.278e-01 5.320 1.04e-07 *** +factor(subject)15 1.765e+00 2.164e-01 8.153 3.54e-16 *** +factor(subject)16 9.708e-01 2.348e-01 4.134 3.57e-05 *** +factor(subject)17 -4.082e-02 2.858e-01 -0.143 0.886411 +factor(subject)18 2.236e+00 2.104e-01 10.629 < 2e-16 *** +factor(subject)19 2.776e-01 2.651e-01 1.047 0.295060 +factor(subject)20 3.646e-01 2.603e-01 1.401 0.161324 +factor(subject)21 3.922e-02 2.801e-01 0.140 0.888645 +factor(subject)22 -8.338e-02 2.889e-01 -0.289 0.772894 +factor(subject)23 1.823e-01 2.708e-01 0.673 0.500777 +factor(subject)24 8.416e-01 2.393e-01 3.517 0.000436 *** +factor(subject)25 2.069e+00 2.123e-01 9.750 < 2e-16 *** +factor(subject)26 -5.108e-01 3.266e-01 -1.564 0.117799 +factor(subject)27 -2.231e-01 3.000e-01 -0.744 0.456990 +factor(subject)28 1.386e+00 2.236e-01 6.200 5.66e-10 *** +factor(subject)29 1.604e+00 2.227e-01 7.203 5.90e-13 *** +factor(subject)30 1.023e+00 2.372e-01 4.313 1.61e-05 *** +factor(subject)31 9.149e-02 2.821e-01 0.324 0.745700 +factor(subject)32 -3.111e-02 2.909e-01 -0.107 0.914822 +factor(subject)33 4.710e-01 2.597e-01 1.814 0.069736 . +factor(subject)34 3.887e-01 2.640e-01 1.473 0.140879 +factor(subject)35 1.487e+00 2.250e-01 6.609 3.87e-11 *** +factor(subject)36 3.598e-01 2.656e-01 1.355 0.175551 +factor(subject)37 -1.221e-01 2.979e-01 -0.410 0.681943 +factor(subject)38 1.344e+00 2.283e-01 5.889 3.90e-09 *** +factor(subject)39 1.082e+00 2.354e-01 4.596 4.30e-06 *** +factor(subject)40 -7.687e-01 3.634e-01 -2.116 0.034384 * +factor(subject)41 1.656e-01 2.772e-01 0.597 0.550234 +factor(subject)42 5.227e-02 2.848e-01 0.184 0.854388 +factor(subject)43 1.543e+00 2.239e-01 6.891 5.54e-12 *** +factor(subject)44 9.605e-01 2.393e-01 4.014 5.96e-05 *** +factor(subject)45 1.177e+00 2.326e-01 5.061 4.18e-07 *** +factor(subject)46 -5.275e-01 3.355e-01 -1.572 0.115840 +factor(subject)47 1.053e+00 2.363e-01 4.456 8.35e-06 *** +factor(subject)48 -5.275e-01 3.355e-01 -1.572 0.115840 +factor(subject)49 2.949e+00 2.082e-01 14.168 < 2e-16 *** +factor(subject)50 3.887e-01 2.640e-01 1.473 0.140879 +factor(subject)51 1.038e+00 2.367e-01 4.385 1.16e-05 *** +factor(subject)52 5.711e-01 2.548e-01 2.241 0.025023 * +factor(subject)53 1.670e+00 2.215e-01 7.538 4.76e-14 *** +factor(subject)54 4.443e-01 2.611e-01 1.702 0.088759 . +factor(subject)55 2.674e-01 2.709e-01 0.987 0.323618 +factor(subject)56 1.124e+00 2.341e-01 4.800 1.59e-06 *** +factor(subject)57 2.674e-01 2.709e-01 0.987 0.323618 +factor(subject)58 -6.017e-01 3.436e-01 -1.751 0.079911 . +factor(subject)59 -7.556e-02 2.942e-01 -0.257 0.797331 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for poisson family taken to be 1) + + Null deviance: 3180.82 on 117 degrees of freedom +Residual deviance: 303.16 on 57 degrees of freedom +AIC: 1003.5 + +Number of Fisher Scoring iterations: 5 + +> ## IGNORE_RDIFF_END +> +> summary(glmmPQL(y ~ lbase*trt + lage + V4, ++ random = ~ 1 | subject, ++ family = poisson, data = epil)) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +Linear mixed-effects model fit by maximum likelihood + Data: epil + AIC BIC logLik + NA NA NA + +Random effects: + Formula: ~1 | subject + (Intercept) Residual +StdDev: 0.4442704 1.400807 + +Variance function: + Structure: fixed weights + Formula: ~invwt +Fixed effects: y ~ lbase * trt + lage + V4 + Value Std.Error DF t-value p-value +(Intercept) 1.8696677 0.1055620 176 17.711554 0.0000 +lbase 0.8818228 0.1292834 54 6.820849 0.0000 +trtprogabide -0.3095253 0.1490438 54 -2.076740 0.0426 +lage 0.5335460 0.3463119 54 1.540652 0.1292 +V4 -0.1597696 0.0774521 176 -2.062819 0.0406 +lbase:trtprogabide 0.3415425 0.2033325 54 1.679725 0.0988 + Correlation: + (Intr) lbase trtprg lage V4 +lbase -0.126 +trtprogabide -0.691 0.089 +lage -0.103 -0.038 0.088 +V4 -0.162 0.000 0.000 0.000 +lbase:trtprogabide 0.055 -0.645 -0.184 0.267 0.000 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-2.13240534 -0.63871136 -0.08486339 0.41960195 4.97872138 + +Number of Observations: 236 +Number of Groups: 59 +> summary(glmmPQL(y ~ pred, random = ~1 | subject, ++ family = poisson, data = epil3)) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +iteration 6 +iteration 7 +iteration 8 +Linear mixed-effects model fit by maximum likelihood + Data: epil3 + AIC BIC logLik + NA NA NA + +Random effects: + Formula: ~1 | subject + (Intercept) Residual +StdDev: 0.7257895 2.16629 + +Variance function: + Structure: fixed weights + Formula: ~invwt +Fixed effects: y ~ pred + Value Std.Error DF t-value p-value +(Intercept) 3.213631 0.10569117 58 30.405865 0.0000 +predplacebo-base 0.110855 0.09989089 57 1.109763 0.2718 +preddrug-placebo -0.105613 0.13480483 57 -0.783450 0.4366 + Correlation: + (Intr) prdpl- +predplacebo-base 0.081 +preddrug-placebo -0.010 -0.700 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-2.0446864 -0.4765135 -0.1975651 0.3145761 2.6532834 + +Number of Observations: 118 +Number of Groups: 59 +> +> +> +> cleanEx() +> nameEx("farms") +> ### * farms +> +> flush(stderr()); flush(stdout()) +> +> ### Name: farms +> ### Title: Ecological Factors in Farm Management +> ### Aliases: farms +> ### Keywords: datasets +> +> ### ** Examples +> +> farms.mca <- mca(farms, abbrev = TRUE) # Use levels as names +> eqscplot(farms.mca$cs, type = "n") +> text(farms.mca$rs, cex = 0.7) +> text(farms.mca$cs, labels = dimnames(farms.mca$cs)[[1]], cex = 0.7) +> +> +> +> cleanEx() +> nameEx("fitdistr") +> ### * fitdistr +> +> flush(stderr()); flush(stdout()) +> +> ### Name: fitdistr +> ### Title: Maximum-likelihood Fitting of Univariate Distributions +> ### Aliases: fitdistr +> ### Keywords: distribution htest +> +> ### ** Examples +> +> ## avoid spurious accuracy +> op <- options(digits = 3) +> set.seed(123) +> x <- rgamma(100, shape = 5, rate = 0.1) +> fitdistr(x, "gamma") + shape rate + 6.4870 0.1365 + (0.8946) (0.0196) +> ## now do this directly with more control. +> fitdistr(x, dgamma, list(shape = 1, rate = 0.1), lower = 0.001) + shape rate + 6.4869 0.1365 + (0.8944) (0.0196) +> +> set.seed(123) +> x2 <- rt(250, df = 9) +> fitdistr(x2, "t", df = 9) + m s + -0.0107 1.0441 + ( 0.0722) ( 0.0543) +> ## allow df to vary: not a very good idea! +> fitdistr(x2, "t") +Warning in dt((x - m)/s, df, log = TRUE) : NaNs produced +Calls: fitdistr ... -> -> fn -> dens -> densfun -> dt + m s df + -0.00965 1.00617 6.62729 + ( 0.07147) ( 0.07707) ( 2.71033) +> ## now do fixed-df fit directly with more control. +> mydt <- function(x, m, s, df) dt((x-m)/s, df)/s +> fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower = c(-Inf, 0)) + m s + -0.0107 1.0441 + ( 0.0722) ( 0.0543) +> +> set.seed(123) +> x3 <- rweibull(100, shape = 4, scale = 100) +> fitdistr(x3, "weibull") + shape scale + 4.080 99.984 + ( 0.313) ( 2.582) +> +> set.seed(123) +> x4 <- rnegbin(500, mu = 5, theta = 4) +> fitdistr(x4, "Negative Binomial") + size mu + 4.216 4.945 + (0.504) (0.147) +> options(op) +> +> +> +> cleanEx() +> nameEx("fractions") +> ### * fractions +> +> flush(stderr()); flush(stdout()) +> +> ### Name: fractions +> ### Title: Rational Approximation +> ### Aliases: fractions Math.fractions Ops.fractions Summary.fractions +> ### [.fractions [<-.fractions as.character.fractions as.fractions +> ### is.fractions print.fractions t.fractions +> ### Keywords: math +> +> ### ** Examples +> +> X <- matrix(runif(25), 5, 5) +> zapsmall(solve(X, X/5)) # print near-zeroes as zero + [,1] [,2] [,3] [,4] [,5] +[1,] 0.2 0.0 0.0 0.0 0.0 +[2,] 0.0 0.2 0.0 0.0 0.0 +[3,] 0.0 0.0 0.2 0.0 0.0 +[4,] 0.0 0.0 0.0 0.2 0.0 +[5,] 0.0 0.0 0.0 0.0 0.2 +> fractions(solve(X, X/5)) + [,1] [,2] [,3] [,4] [,5] +[1,] 1/5 0 0 0 0 +[2,] 0 1/5 0 0 0 +[3,] 0 0 1/5 0 0 +[4,] 0 0 0 1/5 0 +[5,] 0 0 0 0 1/5 +> fractions(solve(X, X/5)) + 1 + [,1] [,2] [,3] [,4] [,5] +[1,] 6/5 1 1 1 1 +[2,] 1 6/5 1 1 1 +[3,] 1 1 6/5 1 1 +[4,] 1 1 1 6/5 1 +[5,] 1 1 1 1 6/5 +> +> +> +> cleanEx() +> nameEx("galaxies") +> ### * galaxies +> +> flush(stderr()); flush(stdout()) +> +> ### Name: galaxies +> ### Title: Velocities for 82 Galaxies +> ### Aliases: galaxies +> ### Keywords: datasets +> +> ### ** Examples +> +> gal <- galaxies/1000 +> c(width.SJ(gal, method = "dpi"), width.SJ(gal)) +[1] 3.256151 2.566423 +> plot(x = c(0, 40), y = c(0, 0.3), type = "n", bty = "l", ++ xlab = "velocity of galaxy (1000km/s)", ylab = "density") +> rug(gal) +> lines(density(gal, width = 3.25, n = 200), lty = 1) +> lines(density(gal, width = 2.56, n = 200), lty = 3) +> +> +> +> cleanEx() +> nameEx("gamma.shape.glm") +> ### * gamma.shape.glm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: gamma.shape +> ### Title: Estimate the Shape Parameter of the Gamma Distribution in a GLM +> ### Fit +> ### Aliases: gamma.shape gamma.shape.glm print.gamma.shape +> ### Keywords: models +> +> ### ** Examples +> +> clotting <- data.frame( ++ u = c(5,10,15,20,30,40,60,80,100), ++ lot1 = c(118,58,42,35,27,25,21,19,18), ++ lot2 = c(69,35,26,21,18,16,13,12,12)) +> clot1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma) +> gamma.shape(clot1) + +Alpha: 538.1315 +SE: 253.5991 +> +> gm <- glm(Days + 0.1 ~ Age*Eth*Sex*Lrn, ++ quasi(link=log, variance="mu^2"), quine, ++ start = c(3, rep(0,31))) +> gamma.shape(gm, verbose = TRUE) +Initial estimate: 1.060344 +Iter. 1 Alpha: 1.238408 +Iter. 2 Alpha: 1.276997 +Iter. 3 Alpha: 1.278343 +Iter. 4 Alpha: 1.278345 + +Alpha: 1.2783449 +SE: 0.1345175 +> ## IGNORE_RDIFF_BEGIN +> summary(gm, dispersion = gamma.dispersion(gm)) # better summary + +Call: +glm(formula = Days + 0.1 ~ Age * Eth * Sex * Lrn, family = quasi(link = log, + variance = "mu^2"), data = quine, start = c(3, rep(0, 31))) + +Coefficients: (4 not defined because of singularities) + Estimate Std. Error z value Pr(>|z|) +(Intercept) 3.06105 0.44223 6.922 4.46e-12 *** +AgeF1 -0.61870 0.59331 -1.043 0.297041 +AgeF2 -2.31911 0.98885 -2.345 0.019014 * +AgeF3 -0.37623 0.53149 -0.708 0.479020 +EthN -0.13789 0.62540 -0.220 0.825496 +SexM -0.48844 0.59331 -0.823 0.410369 +LrnSL -1.92965 0.98885 -1.951 0.051009 . +AgeF1:EthN 0.10249 0.82338 0.124 0.900942 +AgeF2:EthN -0.50874 1.39845 -0.364 0.716017 +AgeF3:EthN 0.06314 0.74584 0.085 0.932534 +AgeF1:SexM 0.40695 0.94847 0.429 0.667884 +AgeF2:SexM 3.06173 1.11626 2.743 0.006091 ** +AgeF3:SexM 1.10841 0.74208 1.494 0.135267 +EthN:SexM -0.74217 0.82338 -0.901 0.367394 +AgeF1:LrnSL 2.60967 1.10114 2.370 0.017789 * +AgeF2:LrnSL 4.78434 1.36304 3.510 0.000448 *** +AgeF3:LrnSL NA NA NA NA +EthN:LrnSL 2.22936 1.39845 1.594 0.110899 +SexM:LrnSL 1.56531 1.18112 1.325 0.185077 +AgeF1:EthN:SexM -0.30235 1.32176 -0.229 0.819065 +AgeF2:EthN:SexM 0.29742 1.57035 0.189 0.849780 +AgeF3:EthN:SexM 0.82215 1.03277 0.796 0.425995 +AgeF1:EthN:LrnSL -3.50803 1.54655 -2.268 0.023311 * +AgeF2:EthN:LrnSL -3.33529 1.92481 -1.733 0.083133 . +AgeF3:EthN:LrnSL NA NA NA NA +AgeF1:SexM:LrnSL -2.39791 1.51050 -1.587 0.112400 +AgeF2:SexM:LrnSL -4.12161 1.60698 -2.565 0.010323 * +AgeF3:SexM:LrnSL NA NA NA NA +EthN:SexM:LrnSL -0.15305 1.66253 -0.092 0.926653 +AgeF1:EthN:SexM:LrnSL 2.13480 2.08685 1.023 0.306317 +AgeF2:EthN:SexM:LrnSL 2.11886 2.27882 0.930 0.352473 +AgeF3:EthN:SexM:LrnSL NA NA NA NA +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for quasi family taken to be 0.7822615) + + Null deviance: 190.40 on 145 degrees of freedom +Residual deviance: 128.36 on 118 degrees of freedom +AIC: NA + +Number of Fisher Scoring iterations: 7 + +> ## IGNORE_RDIFF_END +> +> +> +> cleanEx() +> nameEx("gehan") +> ### * gehan +> +> flush(stderr()); flush(stdout()) +> +> ### Name: gehan +> ### Title: Remission Times of Leukaemia Patients +> ### Aliases: gehan +> ### Keywords: datasets +> +> ### ** Examples +> +> library(survival) +> gehan.surv <- survfit(Surv(time, cens) ~ treat, data = gehan, ++ conf.type = "log-log") +> summary(gehan.surv) +Call: survfit(formula = Surv(time, cens) ~ treat, data = gehan, conf.type = "log-log") + + treat=6-MP + time n.risk n.event survival std.err lower 95% CI upper 95% CI + 6 21 3 0.857 0.0764 0.620 0.952 + 7 17 1 0.807 0.0869 0.563 0.923 + 10 15 1 0.753 0.0963 0.503 0.889 + 13 12 1 0.690 0.1068 0.432 0.849 + 16 11 1 0.627 0.1141 0.368 0.805 + 22 7 1 0.538 0.1282 0.268 0.747 + 23 6 1 0.448 0.1346 0.188 0.680 + + treat=control + time n.risk n.event survival std.err lower 95% CI upper 95% CI + 1 21 2 0.9048 0.0641 0.67005 0.975 + 2 19 2 0.8095 0.0857 0.56891 0.924 + 3 17 1 0.7619 0.0929 0.51939 0.893 + 4 16 2 0.6667 0.1029 0.42535 0.825 + 5 14 2 0.5714 0.1080 0.33798 0.749 + 8 12 4 0.3810 0.1060 0.18307 0.578 + 11 8 2 0.2857 0.0986 0.11656 0.482 + 12 6 2 0.1905 0.0857 0.05948 0.377 + 15 4 1 0.1429 0.0764 0.03566 0.321 + 17 3 1 0.0952 0.0641 0.01626 0.261 + 22 2 1 0.0476 0.0465 0.00332 0.197 + 23 1 1 0.0000 NaN NA NA + +> survreg(Surv(time, cens) ~ factor(pair) + treat, gehan, dist = "exponential") +Call: +survreg(formula = Surv(time, cens) ~ factor(pair) + treat, data = gehan, + dist = "exponential") +Warning in x$coef : partial match of 'coef' to 'coefficients' +Calls: -> print.survreg + +Coefficients: + (Intercept) factor(pair)2 factor(pair)3 factor(pair)4 factor(pair)5 + 2.0702861 2.1476909 1.8329493 1.7718527 1.4682566 + factor(pair)6 factor(pair)7 factor(pair)8 factor(pair)9 factor(pair)10 + 1.8954775 0.5583010 2.5187140 2.2970513 2.4862208 +factor(pair)11 factor(pair)12 factor(pair)13 factor(pair)14 factor(pair)15 + 1.0524472 1.8270477 1.6772567 1.7778672 2.0859913 +factor(pair)16 factor(pair)17 factor(pair)18 factor(pair)19 factor(pair)20 + 3.0634288 0.7996252 1.5855018 1.4083884 0.4023946 +factor(pair)21 treatcontrol + 1.9698390 -1.7671562 + +Scale fixed at 1 +Warning in x$linear : partial match of 'linear' to 'linear.predictors' +Calls: -> print.survreg + +Loglik(model)= -101.6 Loglik(intercept only)= -116.8 + Chisq= 30.27 on 21 degrees of freedom, p= 0.0866 +n= 42 +> summary(survreg(Surv(time, cens) ~ treat, gehan, dist = "exponential")) + +Call: +survreg(formula = Surv(time, cens) ~ treat, data = gehan, dist = "exponential") + Value Std. Error z p +(Intercept) 3.686 0.333 11.06 < 2e-16 +treatcontrol -1.527 0.398 -3.83 0.00013 + +Scale fixed at 1 + +Exponential distribution +Loglik(model)= -108.5 Loglik(intercept only)= -116.8 + Chisq= 16.49 on 1 degrees of freedom, p= 4.9e-05 +Number of Newton-Raphson Iterations: 4 +n= 42 + +> summary(survreg(Surv(time, cens) ~ treat, gehan)) + +Call: +survreg(formula = Surv(time, cens) ~ treat, data = gehan) + Value Std. Error z p +(Intercept) 3.516 0.252 13.96 < 2e-16 +treatcontrol -1.267 0.311 -4.08 4.5e-05 +Log(scale) -0.312 0.147 -2.12 0.034 + +Scale= 0.732 + +Weibull distribution +Loglik(model)= -106.6 Loglik(intercept only)= -116.4 + Chisq= 19.65 on 1 degrees of freedom, p= 9.3e-06 +Number of Newton-Raphson Iterations: 5 +n= 42 + +> gehan.cox <- coxph(Surv(time, cens) ~ treat, gehan) +> summary(gehan.cox) +Call: +coxph(formula = Surv(time, cens) ~ treat, data = gehan) + + n= 42, number of events= 30 +Warning in x$coef : partial match of 'coef' to 'coefficients' +Calls: -> print.summary.coxph -> nrow + + coef exp(coef) se(coef) z Pr(>|z|) +treatcontrol 1.5721 4.8169 0.4124 3.812 0.000138 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + + exp(coef) exp(-coef) lower .95 upper .95 +treatcontrol 4.817 0.2076 2.147 10.81 + +Concordance= 0.69 (se = 0.041 ) +Likelihood ratio test= 16.35 on 1 df, p=5e-05 +Wald test = 14.53 on 1 df, p=1e-04 +Score (logrank) test = 17.25 on 1 df, p=3e-05 + +> +> +> +> cleanEx() + +detaching ‘package:survival’ + +> nameEx("glm.convert") +> ### * glm.convert +> +> flush(stderr()); flush(stdout()) +> +> ### Name: glm.convert +> ### Title: Change a Negative Binomial fit to a GLM fit +> ### Aliases: glm.convert +> ### Keywords: regression models +> +> ### ** Examples +> +> quine.nb1 <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine) +> quine.nbA <- glm.convert(quine.nb1) +> quine.nbB <- update(quine.nb1, . ~ . + Sex:Age:Lrn) +> anova(quine.nbA, quine.nbB) +Analysis of Deviance Table + +Model 1: Days ~ Sex/(Age + Eth * Lrn) +Model 2: Days ~ Sex + Sex:Age + Sex:Eth + Sex:Lrn + Sex:Eth:Lrn + Sex:Age:Lrn + Resid. Df Resid. Dev Df Deviance +1 132 167.56 +2 128 166.83 4 0.723 +> +> +> +> cleanEx() +> nameEx("glm.nb") +> ### * glm.nb +> +> flush(stderr()); flush(stdout()) +> +> ### Name: glm.nb +> ### Title: Fit a Negative Binomial Generalized Linear Model +> ### Aliases: glm.nb family.negbin logLik.negbin +> ### Keywords: regression models +> +> ### ** Examples +> +> quine.nb1 <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine) +> quine.nb2 <- update(quine.nb1, . ~ . + Sex:Age:Lrn) +> quine.nb3 <- update(quine.nb2, Days ~ .^4) +> anova(quine.nb1, quine.nb2, quine.nb3) +Likelihood ratio tests of Negative Binomial Models + +Response: Days + Model +1 Sex/(Age + Eth * Lrn) +2 Sex + Sex:Age + Sex:Eth + Sex:Lrn + Sex:Eth:Lrn + Sex:Age:Lrn +3 Sex + Sex:Age + Sex:Eth + Sex:Lrn + Sex:Eth:Lrn + Sex:Age:Lrn + Sex:Age:Eth + Sex:Age:Eth:Lrn + theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi) +1 1.597991 132 -1063.025 +2 1.686899 128 -1055.398 1 vs 2 4 7.627279 0.10622602 +3 1.928360 118 -1039.324 2 vs 3 10 16.073723 0.09754136 +> ## Don't show: +> ## PR#1695 +> y <- c(7, 5, 4, 7, 5, 2, 11, 5, 5, 4, 2, 3, 4, 3, 5, 9, 6, 7, 10, 6, 12, ++ 6, 3, 5, 3, 9, 13, 0, 6, 1, 2, 0, 1, 0, 0, 4, 5, 1, 5, 3, 3, 4) +> +> lag1 <- c(0, 7, 5, 4, 7, 5, 2, 11, 5, 5, 4, 2, 3, 4, 3, 5, 9, 6, 7, 10, ++ 6, 12, 6, 3, 5, 3, 9, 13, 0, 6, 1, 2, 0, 1, 0, 0, 4, 5, 1, 5, 3, 3) +> +> lag2 <- c(0, 0, 7, 5, 4, 7, 5, 2, 11, 5, 5, 4, 2, 3, 4, 3, 5, 9, 6, 7, ++ 10, 6, 12, 6, 3, 5, 3, 9, 13, 0, 6, 1, 2, 0, 1, 0, 0, 4, 5, 1, 5, 3) +> +> lag3 <- c(0, 0, 0, 7, 5, 4, 7, 5, 2, 11, 5, 5, 4, 2, 3, 4, 3, 5, 9, 6, ++ 7, 10, 6, 12, 6, 3, 5, 3, 9, 13, 0, 6, 1, 2, 0, 1, 0, 0, 4, 5, 1, 5) +> +> (fit <- glm(y ~ lag1+lag2+lag3, family=poisson(link=identity), ++ start=c(2, 0.1, 0.1, 0.1))) + +Call: glm(formula = y ~ lag1 + lag2 + lag3, family = poisson(link = identity), + start = c(2, 0.1, 0.1, 0.1)) + +Coefficients: +(Intercept) lag1 lag2 lag3 + 2.6609 0.1573 0.1424 0.1458 + +Degrees of Freedom: 41 Total (i.e. Null); 38 Residual +Null Deviance: 100.2 +Residual Deviance: 90.34 AIC: 225.6 +> try(glm.nb(y ~ lag1+lag2+lag3, link=identity)) +Warning in log(y/mu) : NaNs produced +Calls: try ... tryCatchOne -> doTryCatch -> glm.nb -> glm.fitter -> dev.resids +Error : no valid set of coefficients has been found: please supply starting values +> glm.nb(y ~ lag1+lag2+lag3, link=identity, start=c(2, 0.1, 0.1, 0.1)) + +Call: glm.nb(formula = y ~ lag1 + lag2 + lag3, start = c(2, 0.1, 0.1, + 0.1), link = identity, init.theta = 4.406504429) + +Coefficients: +(Intercept) lag1 lag2 lag3 + 2.6298 0.1774 0.1407 0.1346 + +Degrees of Freedom: 41 Total (i.e. Null); 38 Residual +Null Deviance: 55.07 +Residual Deviance: 50.09 AIC: 215.9 +> glm.nb(y ~ lag1+lag2+lag3, link=identity, start=coef(fit)) + +Call: glm.nb(formula = y ~ lag1 + lag2 + lag3, start = coef(fit), link = identity, + init.theta = 4.406504429) + +Coefficients: +(Intercept) lag1 lag2 lag3 + 2.6298 0.1774 0.1407 0.1346 + +Degrees of Freedom: 41 Total (i.e. Null); 38 Residual +Null Deviance: 55.07 +Residual Deviance: 50.09 AIC: 215.9 +> glm.nb(y ~ lag1+lag2+lag3, link=identity, etastart=rep(5, 42)) + +Call: glm.nb(formula = y ~ lag1 + lag2 + lag3, etastart = rep(5, 42), + link = identity, init.theta = 4.406504429) + +Coefficients: +(Intercept) lag1 lag2 lag3 + 2.6298 0.1774 0.1407 0.1346 + +Degrees of Freedom: 41 Total (i.e. Null); 38 Residual +Null Deviance: 55.07 +Residual Deviance: 50.09 AIC: 215.9 +> ## End(Don't show) +> +> +> cleanEx() +> nameEx("glmmPQL") +> ### * glmmPQL +> +> flush(stderr()); flush(stdout()) +> +> ### Name: glmmPQL +> ### Title: Fit Generalized Linear Mixed Models via PQL +> ### Aliases: glmmPQL +> ### Keywords: models +> +> ### ** Examples +> +> summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID, ++ family = binomial, data = bacteria)) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +iteration 6 +Linear mixed-effects model fit by maximum likelihood + Data: bacteria + AIC BIC logLik + NA NA NA + +Random effects: + Formula: ~1 | ID + (Intercept) Residual +StdDev: 1.410637 0.7800511 + +Variance function: + Structure: fixed weights + Formula: ~invwt +Fixed effects: y ~ trt + I(week > 2) + Value Std.Error DF t-value p-value +(Intercept) 3.412014 0.5185033 169 6.580506 0.0000 +trtdrug -1.247355 0.6440635 47 -1.936696 0.0588 +trtdrug+ -0.754327 0.6453978 47 -1.168779 0.2484 +I(week > 2)TRUE -1.607257 0.3583379 169 -4.485311 0.0000 + Correlation: + (Intr) trtdrg trtdr+ +trtdrug -0.598 +trtdrug+ -0.571 0.460 +I(week > 2)TRUE -0.537 0.047 -0.001 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-5.1985361 0.1572336 0.3513075 0.4949482 1.7448845 + +Number of Observations: 220 +Number of Groups: 50 +> +> ## an example of an offset: the coefficient of 'week' changes by one. +> summary(glmmPQL(y ~ trt + week, random = ~ 1 | ID, ++ family = binomial, data = bacteria)) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +iteration 6 +Linear mixed-effects model fit by maximum likelihood + Data: bacteria + AIC BIC logLik + NA NA NA + +Random effects: + Formula: ~1 | ID + (Intercept) Residual +StdDev: 1.325243 0.7903088 + +Variance function: + Structure: fixed weights + Formula: ~invwt +Fixed effects: y ~ trt + week + Value Std.Error DF t-value p-value +(Intercept) 3.0302276 0.4791396 169 6.324310 0.0000 +trtdrug -1.2176812 0.6160113 47 -1.976719 0.0540 +trtdrug+ -0.7886376 0.6193895 47 -1.273250 0.2092 +week -0.1446463 0.0392343 169 -3.686730 0.0003 + Correlation: + (Intr) trtdrg trtdr+ +trtdrug -0.622 +trtdrug+ -0.609 0.464 +week -0.481 0.050 0.030 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-4.2868074 0.2039043 0.3140333 0.5440835 1.9754065 + +Number of Observations: 220 +Number of Groups: 50 +> summary(glmmPQL(y ~ trt + week + offset(week), random = ~ 1 | ID, ++ family = binomial, data = bacteria)) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +iteration 6 +Linear mixed-effects model fit by maximum likelihood + Data: bacteria + AIC BIC logLik + NA NA NA + +Random effects: + Formula: ~1 | ID + (Intercept) Residual +StdDev: 1.325243 0.7903088 + +Variance function: + Structure: fixed weights + Formula: ~invwt +Fixed effects: y ~ trt + week + offset(week) + Value Std.Error DF t-value p-value +(Intercept) 3.0302276 0.4791396 169 6.324310 0.0000 +trtdrug -1.2176812 0.6160113 47 -1.976719 0.0540 +trtdrug+ -0.7886376 0.6193895 47 -1.273250 0.2092 +week -1.1446463 0.0392343 169 -29.174622 0.0000 + Correlation: + (Intr) trtdrg trtdr+ +trtdrug -0.622 +trtdrug+ -0.609 0.464 +week -0.481 0.050 0.030 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-4.2868074 0.2039043 0.3140333 0.5440835 1.9754065 + +Number of Observations: 220 +Number of Groups: 50 +> +> +> +> cleanEx() +> nameEx("housing") +> ### * housing +> +> flush(stderr()); flush(stdout()) +> +> ### Name: housing +> ### Title: Frequency Table from a Copenhagen Housing Conditions Survey +> ### Aliases: housing +> ### Keywords: datasets +> +> ### ** Examples +> +> options(contrasts = c("contr.treatment", "contr.poly")) +> +> # Surrogate Poisson models +> house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family = poisson, ++ data = housing) +> ## IGNORE_RDIFF_BEGIN +> summary(house.glm0, correlation = FALSE) + +Call: +glm(formula = Freq ~ Infl * Type * Cont + Sat, family = poisson, + data = housing) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 3.136e+00 1.196e-01 26.225 < 2e-16 *** +InflMedium 2.733e-01 1.586e-01 1.723 0.084868 . +InflHigh -2.054e-01 1.784e-01 -1.152 0.249511 +TypeApartment 3.666e-01 1.555e-01 2.357 0.018403 * +TypeAtrium -7.828e-01 2.134e-01 -3.668 0.000244 *** +TypeTerrace -8.145e-01 2.157e-01 -3.775 0.000160 *** +ContHigh 1.402e-15 1.690e-01 0.000 1.000000 +Sat.L 1.159e-01 4.038e-02 2.871 0.004094 ** +Sat.Q 2.629e-01 4.515e-02 5.824 5.76e-09 *** +InflMedium:TypeApartment -1.177e-01 2.086e-01 -0.564 0.572571 +InflHigh:TypeApartment 1.753e-01 2.279e-01 0.769 0.441783 +InflMedium:TypeAtrium -4.068e-01 3.035e-01 -1.340 0.180118 +InflHigh:TypeAtrium -1.692e-01 3.294e-01 -0.514 0.607433 +InflMedium:TypeTerrace 6.292e-03 2.860e-01 0.022 0.982450 +InflHigh:TypeTerrace -9.305e-02 3.280e-01 -0.284 0.776633 +InflMedium:ContHigh -1.398e-01 2.279e-01 -0.613 0.539715 +InflHigh:ContHigh -6.091e-01 2.800e-01 -2.176 0.029585 * +TypeApartment:ContHigh 5.029e-01 2.109e-01 2.385 0.017083 * +TypeAtrium:ContHigh 6.774e-01 2.751e-01 2.462 0.013811 * +TypeTerrace:ContHigh 1.099e+00 2.675e-01 4.106 4.02e-05 *** +InflMedium:TypeApartment:ContHigh 5.359e-02 2.862e-01 0.187 0.851450 +InflHigh:TypeApartment:ContHigh 1.462e-01 3.380e-01 0.432 0.665390 +InflMedium:TypeAtrium:ContHigh 1.555e-01 3.907e-01 0.398 0.690597 +InflHigh:TypeAtrium:ContHigh 4.782e-01 4.441e-01 1.077 0.281619 +InflMedium:TypeTerrace:ContHigh -4.980e-01 3.671e-01 -1.357 0.174827 +InflHigh:TypeTerrace:ContHigh -4.470e-01 4.545e-01 -0.984 0.325326 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for poisson family taken to be 1) + + Null deviance: 833.66 on 71 degrees of freedom +Residual deviance: 217.46 on 46 degrees of freedom +AIC: 610.43 + +Number of Fisher Scoring iterations: 5 + +> ## IGNORE_RDIFF_END +> +> addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test = "Chisq") +Single term additions + +Model: +Freq ~ Infl * Type * Cont + Sat + Df Deviance AIC LRT Pr(Chi) + 217.46 610.43 +Infl:Sat 4 111.08 512.05 106.371 < 2.2e-16 *** +Type:Sat 6 156.79 561.76 60.669 3.292e-11 *** +Cont:Sat 2 212.33 609.30 5.126 0.07708 . +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) +> ## IGNORE_RDIFF_BEGIN +> summary(house.glm1, correlation = FALSE) + +Call: +glm(formula = Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + + Type:Cont + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont, + family = poisson, data = housing) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) 3.135074 0.120112 26.101 < 2e-16 *** +InflMedium 0.248327 0.159979 1.552 0.120602 +InflHigh -0.412645 0.184947 -2.231 0.025671 * +TypeApartment 0.292524 0.157477 1.858 0.063231 . +TypeAtrium -0.792847 0.214413 -3.698 0.000218 *** +TypeTerrace -1.018074 0.221263 -4.601 4.20e-06 *** +ContHigh -0.001407 0.169711 -0.008 0.993385 +Sat.L -0.098106 0.112592 -0.871 0.383570 +Sat.Q 0.285657 0.122283 2.336 0.019489 * +InflMedium:TypeApartment -0.017882 0.210496 -0.085 0.932302 +InflHigh:TypeApartment 0.386869 0.233297 1.658 0.097263 . +InflMedium:TypeAtrium -0.360311 0.304979 -1.181 0.237432 +InflHigh:TypeAtrium -0.036788 0.334793 -0.110 0.912503 +InflMedium:TypeTerrace 0.185154 0.288892 0.641 0.521580 +InflHigh:TypeTerrace 0.310749 0.334815 0.928 0.353345 +InflMedium:ContHigh -0.200060 0.228748 -0.875 0.381799 +InflHigh:ContHigh -0.725790 0.282352 -2.571 0.010155 * +TypeApartment:ContHigh 0.569691 0.212152 2.685 0.007247 ** +TypeAtrium:ContHigh 0.702115 0.276056 2.543 0.010979 * +TypeTerrace:ContHigh 1.215930 0.269968 4.504 6.67e-06 *** +InflMedium:Sat.L 0.519627 0.096830 5.366 8.03e-08 *** +InflHigh:Sat.L 1.140302 0.118180 9.649 < 2e-16 *** +InflMedium:Sat.Q -0.064474 0.102666 -0.628 0.530004 +InflHigh:Sat.Q 0.115436 0.127798 0.903 0.366380 +TypeApartment:Sat.L -0.520170 0.109793 -4.738 2.16e-06 *** +TypeAtrium:Sat.L -0.288484 0.149551 -1.929 0.053730 . +TypeTerrace:Sat.L -0.998666 0.141527 -7.056 1.71e-12 *** +TypeApartment:Sat.Q 0.055418 0.118515 0.468 0.640068 +TypeAtrium:Sat.Q -0.273820 0.149713 -1.829 0.067405 . +TypeTerrace:Sat.Q -0.032328 0.149251 -0.217 0.828520 +ContHigh:Sat.L 0.340703 0.087778 3.881 0.000104 *** +ContHigh:Sat.Q -0.097929 0.094068 -1.041 0.297851 +InflMedium:TypeApartment:ContHigh 0.046900 0.286212 0.164 0.869837 +InflHigh:TypeApartment:ContHigh 0.126229 0.338208 0.373 0.708979 +InflMedium:TypeAtrium:ContHigh 0.157239 0.390719 0.402 0.687364 +InflHigh:TypeAtrium:ContHigh 0.478611 0.444244 1.077 0.281320 +InflMedium:TypeTerrace:ContHigh -0.500162 0.367135 -1.362 0.173091 +InflHigh:TypeTerrace:ContHigh -0.463099 0.454713 -1.018 0.308467 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for poisson family taken to be 1) + + Null deviance: 833.657 on 71 degrees of freedom +Residual deviance: 38.662 on 34 degrees of freedom +AIC: 455.63 + +Number of Fisher Scoring iterations: 4 + +> ## IGNORE_RDIFF_END +> +> 1 - pchisq(deviance(house.glm1), house.glm1$df.residual) +[1] 0.2671363 +> +> dropterm(house.glm1, test = "Chisq") +Single term deletions + +Model: +Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont + + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont + Df Deviance AIC LRT Pr(Chi) + 38.662 455.63 +Infl:Sat 4 147.780 556.75 109.117 < 2.2e-16 *** +Type:Sat 6 100.889 505.86 62.227 1.586e-11 *** +Cont:Sat 2 54.722 467.69 16.060 0.0003256 *** +Infl:Type:Cont 6 43.952 448.92 5.290 0.5072454 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test = "Chisq") +Single term additions + +Model: +Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont + + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont + Df Deviance AIC LRT Pr(Chi) + 38.662 455.63 +Infl:Type:Sat 12 16.107 457.08 22.5550 0.03175 * +Infl:Cont:Sat 4 37.472 462.44 1.1901 0.87973 +Type:Cont:Sat 6 28.256 457.23 10.4064 0.10855 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> +> hnames <- lapply(housing[, -5], levels) # omit Freq +> newData <- expand.grid(hnames) +> newData$Sat <- ordered(newData$Sat) +> house.pm <- predict(house.glm1, newData, ++ type = "response") # poisson means +> house.pm <- matrix(house.pm, ncol = 3, byrow = TRUE, ++ dimnames = list(NULL, hnames[[1]])) +> house.pr <- house.pm/drop(house.pm %*% rep(1, 3)) +> cbind(expand.grid(hnames[-1]), round(house.pr, 2)) + Infl Type Cont Low Medium High +1 Low Tower Low 0.40 0.26 0.34 +2 Medium Tower Low 0.26 0.27 0.47 +3 High Tower Low 0.15 0.19 0.66 +4 Low Apartment Low 0.54 0.23 0.23 +5 Medium Apartment Low 0.39 0.26 0.34 +6 High Apartment Low 0.26 0.21 0.53 +7 Low Atrium Low 0.43 0.32 0.25 +8 Medium Atrium Low 0.30 0.35 0.36 +9 High Atrium Low 0.19 0.27 0.54 +10 Low Terrace Low 0.65 0.22 0.14 +11 Medium Terrace Low 0.51 0.27 0.22 +12 High Terrace Low 0.37 0.24 0.39 +13 Low Tower High 0.30 0.28 0.42 +14 Medium Tower High 0.18 0.27 0.54 +15 High Tower High 0.10 0.19 0.71 +16 Low Apartment High 0.44 0.27 0.30 +17 Medium Apartment High 0.30 0.28 0.42 +18 High Apartment High 0.18 0.21 0.61 +19 Low Atrium High 0.33 0.36 0.31 +20 Medium Atrium High 0.22 0.36 0.42 +21 High Atrium High 0.13 0.27 0.60 +22 Low Terrace High 0.55 0.27 0.19 +23 Medium Terrace High 0.40 0.31 0.29 +24 High Terrace High 0.27 0.26 0.47 +> +> # Iterative proportional scaling +> loglm(Freq ~ Infl*Type*Cont + Sat*(Infl+Type+Cont), data = housing) +Call: +loglm(formula = Freq ~ Infl * Type * Cont + Sat * (Infl + Type + + Cont), data = housing) + +Statistics: + X^2 df P(> X^2) +Likelihood Ratio 38.66222 34 0.2671359 +Pearson 38.90831 34 0.2582333 +> +> +> # multinomial model +> library(nnet) +> (house.mult<- multinom(Sat ~ Infl + Type + Cont, weights = Freq, ++ data = housing)) +# weights: 24 (14 variable) +initial value 1846.767257 +iter 10 value 1747.045232 +final value 1735.041933 +converged +Call: +multinom(formula = Sat ~ Infl + Type + Cont, data = housing, + weights = Freq) + +Coefficients: + (Intercept) InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace +Medium -0.4192316 0.4464003 0.6649367 -0.4356851 0.1313663 -0.6665728 +High -0.1387453 0.7348626 1.6126294 -0.7356261 -0.4079808 -1.4123333 + ContHigh +Medium 0.3608513 +High 0.4818236 + +Residual Deviance: 3470.084 +AIC: 3498.084 +> house.mult2 <- multinom(Sat ~ Infl*Type*Cont, weights = Freq, ++ data = housing) +# weights: 75 (48 variable) +initial value 1846.767257 +iter 10 value 1734.465581 +iter 20 value 1717.220153 +iter 30 value 1715.760679 +iter 40 value 1715.713306 +final value 1715.710836 +converged +> anova(house.mult, house.mult2) +Likelihood ratio tests of Multinomial Models + +Response: Sat + Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi) +1 Infl + Type + Cont 130 3470.084 +2 Infl * Type * Cont 96 3431.422 1 vs 2 34 38.66219 0.2671367 +> +> house.pm <- predict(house.mult, expand.grid(hnames[-1]), type = "probs") +> cbind(expand.grid(hnames[-1]), round(house.pm, 2)) + Infl Type Cont Low Medium High +1 Low Tower Low 0.40 0.26 0.34 +2 Medium Tower Low 0.26 0.27 0.47 +3 High Tower Low 0.15 0.19 0.66 +4 Low Apartment Low 0.54 0.23 0.23 +5 Medium Apartment Low 0.39 0.26 0.34 +6 High Apartment Low 0.26 0.21 0.53 +7 Low Atrium Low 0.43 0.32 0.25 +8 Medium Atrium Low 0.30 0.35 0.36 +9 High Atrium Low 0.19 0.27 0.54 +10 Low Terrace Low 0.65 0.22 0.14 +11 Medium Terrace Low 0.51 0.27 0.22 +12 High Terrace Low 0.37 0.24 0.39 +13 Low Tower High 0.30 0.28 0.42 +14 Medium Tower High 0.18 0.27 0.54 +15 High Tower High 0.10 0.19 0.71 +16 Low Apartment High 0.44 0.27 0.30 +17 Medium Apartment High 0.30 0.28 0.42 +18 High Apartment High 0.18 0.21 0.61 +19 Low Atrium High 0.33 0.36 0.31 +20 Medium Atrium High 0.22 0.36 0.42 +21 High Atrium High 0.13 0.27 0.60 +22 Low Terrace High 0.55 0.27 0.19 +23 Medium Terrace High 0.40 0.31 0.29 +24 High Terrace High 0.27 0.26 0.47 +> +> # proportional odds model +> house.cpr <- apply(house.pr, 1, cumsum) +> logit <- function(x) log(x/(1-x)) +> house.ld <- logit(house.cpr[2, ]) - logit(house.cpr[1, ]) +> (ratio <- sort(drop(house.ld))) + [1] 0.9357341 0.9854433 1.0573182 1.0680491 1.0772649 1.0803574 1.0824895 + [8] 1.0998759 1.1199975 1.1554228 1.1768138 1.1866427 1.2091541 1.2435026 +[15] 1.2724096 1.2750171 1.2849903 1.3062598 1.3123988 1.3904715 1.4540087 +[22] 1.4947753 1.4967585 1.6068789 +> mean(ratio) +[1] 1.223835 +> +> (house.plr <- polr(Sat ~ Infl + Type + Cont, ++ data = housing, weights = Freq)) +Call: +polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq) + +Coefficients: + InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace + 0.5663937 1.2888191 -0.5723501 -0.3661866 -1.0910149 + ContHigh + 0.3602841 + +Intercepts: + Low|Medium Medium|High + -0.4961353 0.6907083 + +Residual Deviance: 3479.149 +AIC: 3495.149 +> +> house.pr1 <- predict(house.plr, expand.grid(hnames[-1]), type = "probs") +> cbind(expand.grid(hnames[-1]), round(house.pr1, 2)) + Infl Type Cont Low Medium High +1 Low Tower Low 0.38 0.29 0.33 +2 Medium Tower Low 0.26 0.27 0.47 +3 High Tower Low 0.14 0.21 0.65 +4 Low Apartment Low 0.52 0.26 0.22 +5 Medium Apartment Low 0.38 0.29 0.33 +6 High Apartment Low 0.23 0.26 0.51 +7 Low Atrium Low 0.47 0.27 0.26 +8 Medium Atrium Low 0.33 0.29 0.38 +9 High Atrium Low 0.19 0.25 0.56 +10 Low Terrace Low 0.64 0.21 0.14 +11 Medium Terrace Low 0.51 0.26 0.23 +12 High Terrace Low 0.33 0.29 0.38 +13 Low Tower High 0.30 0.28 0.42 +14 Medium Tower High 0.19 0.25 0.56 +15 High Tower High 0.10 0.17 0.72 +16 Low Apartment High 0.43 0.28 0.29 +17 Medium Apartment High 0.30 0.28 0.42 +18 High Apartment High 0.17 0.23 0.60 +19 Low Atrium High 0.38 0.29 0.33 +20 Medium Atrium High 0.26 0.27 0.47 +21 High Atrium High 0.14 0.21 0.64 +22 Low Terrace High 0.56 0.25 0.19 +23 Medium Terrace High 0.42 0.28 0.30 +24 High Terrace High 0.26 0.27 0.47 +> +> Fr <- matrix(housing$Freq, ncol = 3, byrow = TRUE) +> 2*sum(Fr*log(house.pr/house.pr1)) +[1] 9.065433 +> +> house.plr2 <- stepAIC(house.plr, ~.^2) +Start: AIC=3495.15 +Sat ~ Infl + Type + Cont + + Df AIC ++ Infl:Type 6 3484.6 ++ Type:Cont 3 3492.5 + 3495.1 ++ Infl:Cont 2 3498.9 +- Cont 1 3507.5 +- Type 3 3545.1 +- Infl 2 3599.4 + +Step: AIC=3484.64 +Sat ~ Infl + Type + Cont + Infl:Type + + Df AIC ++ Type:Cont 3 3482.7 + 3484.6 ++ Infl:Cont 2 3488.5 +- Infl:Type 6 3495.1 +- Cont 1 3497.8 + +Step: AIC=3482.69 +Sat ~ Infl + Type + Cont + Infl:Type + Type:Cont + + Df AIC + 3482.7 +- Type:Cont 3 3484.6 ++ Infl:Cont 2 3486.6 +- Infl:Type 6 3492.5 +> house.plr2$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +Sat ~ Infl + Type + Cont + +Final Model: +Sat ~ Infl + Type + Cont + Infl:Type + Type:Cont + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 1673 3479.149 3495.149 +2 + Infl:Type 6 22.509347 1667 3456.640 3484.640 +3 + Type:Cont 3 7.945029 1664 3448.695 3482.695 +> +> +> +> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly")) +> cleanEx() + +detaching ‘package:nnet’ + +> nameEx("huber") +> ### * huber +> +> flush(stderr()); flush(stdout()) +> +> ### Name: huber +> ### Title: Huber M-estimator of Location with MAD Scale +> ### Aliases: huber +> ### Keywords: robust +> +> ### ** Examples +> +> huber(chem) +$mu +[1] 3.206724 + +$s +[1] 0.526323 + +> +> +> +> cleanEx() +> nameEx("hubers") +> ### * hubers +> +> flush(stderr()); flush(stdout()) +> +> ### Name: hubers +> ### Title: Huber Proposal 2 Robust Estimator of Location and/or Scale +> ### Aliases: hubers +> ### Keywords: robust +> +> ### ** Examples +> +> hubers(chem) +$mu +[1] 3.205498 + +$s +[1] 0.673652 + +> hubers(chem, mu=3.68) +$mu +[1] 3.68 + +$s +[1] 0.9409628 + +> +> +> +> cleanEx() +> nameEx("immer") +> ### * immer +> +> flush(stderr()); flush(stdout()) +> +> ### Name: immer +> ### Title: Yields from a Barley Field Trial +> ### Aliases: immer +> ### Keywords: datasets +> +> ### ** Examples +> +> immer.aov <- aov(cbind(Y1,Y2) ~ Loc + Var, data = immer) +> summary(immer.aov) + Response Y1 : + Df Sum Sq Mean Sq F value Pr(>F) +Loc 5 17829.8 3566.0 21.8923 1.751e-07 *** +Var 4 2756.6 689.2 4.2309 0.01214 * +Residuals 20 3257.7 162.9 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + + Response Y2 : + Df Sum Sq Mean Sq F value Pr(>F) +Loc 5 10285.0 2056.99 10.3901 5.049e-05 *** +Var 4 2845.2 711.29 3.5928 0.02306 * +Residuals 20 3959.5 197.98 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +> +> immer.aov <- aov((Y1+Y2)/2 ~ Var + Loc, data = immer) +> summary(immer.aov) + Df Sum Sq Mean Sq F value Pr(>F) +Var 4 2655 663.7 5.989 0.00245 ** +Loc 5 10610 2122.1 19.148 5.21e-07 *** +Residuals 20 2217 110.8 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> model.tables(immer.aov, type = "means", se = TRUE, cterms = "Var") +Tables of means +Grand mean + +101.09 + + Var +Var + M P S T V + 94.39 102.54 91.13 118.20 99.18 + +Standard errors for differences of means + Var + 6.078 +replic. 6 +> +> +> +> cleanEx() +> nameEx("isoMDS") +> ### * isoMDS +> +> flush(stderr()); flush(stdout()) +> +> ### Name: isoMDS +> ### Title: Kruskal's Non-metric Multidimensional Scaling +> ### Aliases: isoMDS Shepard +> ### Keywords: multivariate +> +> ### ** Examples +> +> swiss.x <- as.matrix(swiss[, -1]) +> swiss.dist <- dist(swiss.x) +> swiss.mds <- isoMDS(swiss.dist) +initial value 2.979731 +iter 5 value 2.431486 +iter 10 value 2.343353 +final value 2.338839 +converged +> plot(swiss.mds$points, type = "n") +> text(swiss.mds$points, labels = as.character(1:nrow(swiss.x))) +> swiss.sh <- Shepard(swiss.dist, swiss.mds$points) +> plot(swiss.sh, pch = ".") +> lines(swiss.sh$x, swiss.sh$yf, type = "S") +> +> +> +> cleanEx() +> nameEx("kde2d") +> ### * kde2d +> +> flush(stderr()); flush(stdout()) +> +> ### Name: kde2d +> ### Title: Two-Dimensional Kernel Density Estimation +> ### Aliases: kde2d +> ### Keywords: dplot +> +> ### ** Examples +> +> attach(geyser) +> plot(duration, waiting, xlim = c(0.5,6), ylim = c(40,100)) +> f1 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6, 40, 100)) +> image(f1, zlim = c(0, 0.05)) +> f2 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6, 40, 100), ++ h = c(width.SJ(duration), width.SJ(waiting)) ) +> image(f2, zlim = c(0, 0.05)) +> persp(f2, phi = 30, theta = 20, d = 5) +> +> plot(duration[-272], duration[-1], xlim = c(0.5, 6), ++ ylim = c(1, 6),xlab = "previous duration", ylab = "duration") +> f1 <- kde2d(duration[-272], duration[-1], ++ h = rep(1.5, 2), n = 50, lims = c(0.5, 6, 0.5, 6)) +> contour(f1, xlab = "previous duration", ++ ylab = "duration", levels = c(0.05, 0.1, 0.2, 0.4) ) +> f1 <- kde2d(duration[-272], duration[-1], ++ h = rep(0.6, 2), n = 50, lims = c(0.5, 6, 0.5, 6)) +> contour(f1, xlab = "previous duration", ++ ylab = "duration", levels = c(0.05, 0.1, 0.2, 0.4) ) +> f1 <- kde2d(duration[-272], duration[-1], ++ h = rep(0.4, 2), n = 50, lims = c(0.5, 6, 0.5, 6)) +> contour(f1, xlab = "previous duration", ++ ylab = "duration", levels = c(0.05, 0.1, 0.2, 0.4) ) +> detach("geyser") +> +> +> +> cleanEx() +> nameEx("lda") +> ### * lda +> +> flush(stderr()); flush(stdout()) +> +> ### Name: lda +> ### Title: Linear Discriminant Analysis +> ### Aliases: lda lda.default lda.data.frame lda.formula lda.matrix +> ### model.frame.lda print.lda coef.lda +> ### Keywords: multivariate +> +> ### ** Examples +> +> Iris <- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]), ++ Sp = rep(c("s","c","v"), rep(50,3))) +> train <- sample(1:150, 75) +> table(Iris$Sp[train]) + + c s v +20 28 27 +> ## your answer may differ +> ## c s v +> ## 22 23 30 +> z <- lda(Sp ~ ., Iris, prior = c(1,1,1)/3, subset = train) +> predict(z, Iris[-train, ])$class + [1] s s s s s s s s s s s s s s s s s s s s s s c c c c c c c c c c c c c c c v +[39] c c c c c c c c c c c c c c v v v v v v v v v v v v v v c v v v v v v v v +Levels: c s v +> ## [1] s s s s s s s s s s s s s s s s s s s s s s s s s s s c c c +> ## [31] c c c c c c c v c c c c v c c c c c c c c c c c c v v v v v +> ## [61] v v v v v v v v v v v v v v v +> (z1 <- update(z, . ~ . - Petal.W.)) +Call: +lda(Sp ~ Sepal.L. + Sepal.W. + Petal.L., data = Iris, prior = c(1, + 1, 1)/3, subset = train) + +Prior probabilities of groups: + c s v +0.3333333 0.3333333 0.3333333 + +Group means: + Sepal.L. Sepal.W. Petal.L. +c 5.975000 2.810000 4.395000 +s 4.978571 3.432143 1.460714 +v 6.748148 2.988889 5.637037 + +Coefficients of linear discriminants: + LD1 LD2 +Sepal.L. 1.1643015 0.68235619 +Sepal.W. 0.7945307 2.23093702 +Petal.L. -3.0421425 0.01236265 + +Proportion of trace: + LD1 LD2 +0.9929 0.0071 +> +> +> +> cleanEx() +> nameEx("leuk") +> ### * leuk +> +> flush(stderr()); flush(stdout()) +> +> ### Name: leuk +> ### Title: Survival Times and White Blood Counts for Leukaemia Patients +> ### Aliases: leuk +> ### Keywords: datasets +> +> ### ** Examples +> +> library(survival) +> plot(survfit(Surv(time) ~ ag, data = leuk), lty = 2:3, col = 2:3) +> +> # now Cox models +> leuk.cox <- coxph(Surv(time) ~ ag + log(wbc), leuk) +> summary(leuk.cox) +Call: +coxph(formula = Surv(time) ~ ag + log(wbc), data = leuk) + + n= 33, number of events= 33 +Warning in x$coef : partial match of 'coef' to 'coefficients' +Calls: -> print.summary.coxph -> nrow + + coef exp(coef) se(coef) z Pr(>|z|) +agpresent -1.0691 0.3433 0.4293 -2.490 0.01276 * +log(wbc) 0.3677 1.4444 0.1360 2.703 0.00687 ** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + + exp(coef) exp(-coef) lower .95 upper .95 +agpresent 0.3433 2.9126 0.148 0.7964 +log(wbc) 1.4444 0.6923 1.106 1.8857 + +Concordance= 0.726 (se = 0.047 ) +Likelihood ratio test= 15.64 on 2 df, p=4e-04 +Wald test = 15.06 on 2 df, p=5e-04 +Score (logrank) test = 16.49 on 2 df, p=3e-04 + +> +> +> +> cleanEx() + +detaching ‘package:survival’ + +> nameEx("lm.ridge") +> ### * lm.ridge +> +> flush(stderr()); flush(stdout()) +> +> ### Name: lm.ridge +> ### Title: Ridge Regression +> ### Aliases: lm.ridge plot.ridgelm print.ridgelm select select.ridgelm +> ### Keywords: models +> +> ### ** Examples +> +> longley # not the same as the S-PLUS dataset + GNP.deflator GNP Unemployed Armed.Forces Population Year Employed +1947 83.0 234.289 235.6 159.0 107.608 1947 60.323 +1948 88.5 259.426 232.5 145.6 108.632 1948 61.122 +1949 88.2 258.054 368.2 161.6 109.773 1949 60.171 +1950 89.5 284.599 335.1 165.0 110.929 1950 61.187 +1951 96.2 328.975 209.9 309.9 112.075 1951 63.221 +1952 98.1 346.999 193.2 359.4 113.270 1952 63.639 +1953 99.0 365.385 187.0 354.7 115.094 1953 64.989 +1954 100.0 363.112 357.8 335.0 116.219 1954 63.761 +1955 101.2 397.469 290.4 304.8 117.388 1955 66.019 +1956 104.6 419.180 282.2 285.7 118.734 1956 67.857 +1957 108.4 442.769 293.6 279.8 120.445 1957 68.169 +1958 110.8 444.546 468.1 263.7 121.950 1958 66.513 +1959 112.6 482.704 381.3 255.2 123.366 1959 68.655 +1960 114.2 502.601 393.1 251.4 125.368 1960 69.564 +1961 115.7 518.173 480.6 257.2 127.852 1961 69.331 +1962 116.9 554.894 400.7 282.7 130.081 1962 70.551 +> names(longley)[1] <- "y" +> lm.ridge(y ~ ., longley) + GNP Unemployed Armed.Forces Population +2946.85636017 0.26352725 0.03648291 0.01116105 -1.73702984 + Year Employed + -1.41879853 0.23128785 +> plot(lm.ridge(y ~ ., longley, ++ lambda = seq(0,0.1,0.001))) +> select(lm.ridge(y ~ ., longley, ++ lambda = seq(0,0.1,0.0001))) +modified HKB estimator is 0.006836982 +modified L-W estimator is 0.05267247 +smallest value of GCV at 0.0057 +> +> +> +> cleanEx() +> nameEx("loglm") +> ### * loglm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: loglm +> ### Title: Fit Log-Linear Models by Iterative Proportional Scaling +> ### Aliases: loglm +> ### Keywords: category models +> +> ### ** Examples +> +> # The data frames Cars93, minn38 and quine are available +> # in the MASS package. +> +> # Case 1: frequencies specified as an array. +> sapply(minn38, function(x) length(levels(x))) + hs phs fol sex f + 3 4 7 2 0 +> ## hs phs fol sex f +> ## 3 4 7 2 0 +> ##minn38a <- array(0, c(3,4,7,2), lapply(minn38[, -5], levels)) +> ##minn38a[data.matrix(minn38[,-5])] <- minn38$f +> +> ## or more simply +> minn38a <- xtabs(f ~ ., minn38) +> +> fm <- loglm(~ 1 + 2 + 3 + 4, minn38a) # numerals as names. +> deviance(fm) +[1] 3711.914 +> ## [1] 3711.9 +> fm1 <- update(fm, .~.^2) +> fm2 <- update(fm, .~.^3, print = TRUE) +5 iterations: deviation 0.07512432 +> ## 5 iterations: deviation 0.075 +> anova(fm, fm1, fm2) +LR tests for hierarchical log-linear models + +Model 1: + ~1 + 2 + 3 + 4 +Model 2: + . ~ `1` + `2` + `3` + `4` + `1`:`2` + `1`:`3` + `1`:`4` + `2`:`3` + `2`:`4` + `3`:`4` +Model 3: + . ~ `1` + `2` + `3` + `4` + `1`:`2` + `1`:`3` + `1`:`4` + `2`:`3` + `2`:`4` + `3`:`4` + `1`:`2`:`3` + `1`:`2`:`4` + `1`:`3`:`4` + `2`:`3`:`4` + + Deviance df Delta(Dev) Delta(df) P(> Delta(Dev) +Model 1 3711.91367 155 +Model 2 220.04285 108 3491.87082 47 0.00000 +Model 3 47.74492 36 172.29794 72 0.00000 +Saturated 0.00000 0 47.74492 36 0.09114 +> +> # Case 1. An array generated with xtabs. +> +> loglm(~ Type + Origin, xtabs(~ Type + Origin, Cars93)) +Call: +loglm(formula = ~Type + Origin, data = xtabs(~Type + Origin, + Cars93)) + +Statistics: + X^2 df P(> X^2) +Likelihood Ratio 18.36179 5 0.00252554 +Pearson 14.07985 5 0.01511005 +> +> # Case 2. Frequencies given as a vector in a data frame +> names(quine) +[1] "Eth" "Sex" "Age" "Lrn" "Days" +> ## [1] "Eth" "Sex" "Age" "Lrn" "Days" +> fm <- loglm(Days ~ .^2, quine) +> gm <- glm(Days ~ .^2, poisson, quine) # check glm. +> c(deviance(fm), deviance(gm)) # deviances agree +[1] 1368.669 1368.669 +> ## [1] 1368.7 1368.7 +> c(fm$df, gm$df) # resid df do not! +[1] 127 +> c(fm$df, gm$df.residual) # resid df do not! +[1] 127 128 +> ## [1] 127 128 +> # The loglm residual degrees of freedom is wrong because of +> # a non-detectable redundancy in the model matrix. +> +> +> +> cleanEx() +> nameEx("logtrans") +> ### * logtrans +> +> flush(stderr()); flush(stdout()) +> +> ### Name: logtrans +> ### Title: Estimate log Transformation Parameter +> ### Aliases: logtrans logtrans.formula logtrans.lm logtrans.default +> ### Keywords: regression models hplot +> +> ### ** Examples +> +> logtrans(Days ~ Age*Sex*Eth*Lrn, data = quine, ++ alpha = seq(0.75, 6.5, length.out = 20)) +> +> +> +> cleanEx() +> nameEx("lqs") +> ### * lqs +> +> flush(stderr()); flush(stdout()) +> +> ### Name: lqs +> ### Title: Resistant Regression +> ### Aliases: lqs lqs.formula lqs.default lmsreg ltsreg +> ### Keywords: models robust +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> set.seed(123) # make reproducible +> lqs(stack.loss ~ ., data = stackloss) +Call: +lqs.formula(formula = stack.loss ~ ., data = stackloss) + +Coefficients: +(Intercept) Air.Flow Water.Temp Acid.Conc. + -3.631e+01 7.292e-01 4.167e-01 -1.659e-16 + +Scale estimates 0.9149 1.0148 + +> lqs(stack.loss ~ ., data = stackloss, method = "S", nsamp = "exact") +Call: +lqs.formula(formula = stack.loss ~ ., data = stackloss, nsamp = "exact", + method = "S") + +Coefficients: +(Intercept) Air.Flow Water.Temp Acid.Conc. + -35.37611 0.82522 0.44248 -0.07965 + +Scale estimates 1.912 + +> ## IGNORE_RDIFF_END +> +> +> +> cleanEx() +> nameEx("mca") +> ### * mca +> +> flush(stderr()); flush(stdout()) +> +> ### Name: mca +> ### Title: Multiple Correspondence Analysis +> ### Aliases: mca print.mca +> ### Keywords: category multivariate +> +> ### ** Examples +> +> farms.mca <- mca(farms, abbrev=TRUE) +> farms.mca +Call: +mca(df = farms, abbrev = TRUE) + +Multiple correspondence analysis of 20 cases of 4 factors + +Correlations 0.806 0.745 cumulative % explained 26.87 51.71 +> plot(farms.mca) +> +> +> +> cleanEx() +> nameEx("menarche") +> ### * menarche +> +> flush(stderr()); flush(stdout()) +> +> ### Name: menarche +> ### Title: Age of Menarche in Warsaw +> ### Aliases: menarche +> ### Keywords: datasets +> +> ### ** Examples +> +> mprob <- glm(cbind(Menarche, Total - Menarche) ~ Age, ++ binomial(link = probit), data = menarche) +> +> +> +> cleanEx() +> nameEx("motors") +> ### * motors +> +> flush(stderr()); flush(stdout()) +> +> ### Name: motors +> ### Title: Accelerated Life Testing of Motorettes +> ### Aliases: motors +> ### Keywords: datasets +> +> ### ** Examples +> +> library(survival) +> plot(survfit(Surv(time, cens) ~ factor(temp), motors), conf.int = FALSE) +> # fit Weibull model +> motor.wei <- survreg(Surv(time, cens) ~ temp, motors) +> summary(motor.wei) + +Call: +survreg(formula = Surv(time, cens) ~ temp, data = motors) + Value Std. Error z p +(Intercept) 16.31852 0.62296 26.2 < 2e-16 +temp -0.04531 0.00319 -14.2 < 2e-16 +Log(scale) -1.09564 0.21480 -5.1 3.4e-07 + +Scale= 0.334 + +Weibull distribution +Loglik(model)= -147.4 Loglik(intercept only)= -169.5 + Chisq= 44.32 on 1 degrees of freedom, p= 2.8e-11 +Number of Newton-Raphson Iterations: 7 +n= 40 + +> # and predict at 130C +> unlist(predict(motor.wei, data.frame(temp=130), se.fit = TRUE)) + fit.1 se.fit.1 +33813.06 7506.36 +> +> motor.cox <- coxph(Surv(time, cens) ~ temp, motors) +> summary(motor.cox) +Call: +coxph(formula = Surv(time, cens) ~ temp, data = motors) + + n= 40, number of events= 17 +Warning in x$coef : partial match of 'coef' to 'coefficients' +Calls: -> print.summary.coxph -> nrow + + coef exp(coef) se(coef) z Pr(>|z|) +temp 0.09185 1.09620 0.02736 3.358 0.000786 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + + exp(coef) exp(-coef) lower .95 upper .95 +temp 1.096 0.9122 1.039 1.157 + +Concordance= 0.84 (se = 0.035 ) +Likelihood ratio test= 25.56 on 1 df, p=4e-07 +Wald test = 11.27 on 1 df, p=8e-04 +Score (logrank) test = 22.73 on 1 df, p=2e-06 + +> # predict at temperature 200 +> plot(survfit(motor.cox, newdata = data.frame(temp=200), ++ conf.type = "log-log")) +> summary( survfit(motor.cox, newdata = data.frame(temp=130)) ) +Call: survfit(formula = motor.cox, newdata = data.frame(temp = 130)) + + time n.risk n.event survival std.err lower 95% CI upper 95% CI + 408 40 4 1.000 0.000254 0.999 1 + 504 36 3 1.000 0.000498 0.999 1 + 1344 28 2 0.999 0.001910 0.995 1 + 1440 26 1 0.998 0.002697 0.993 1 + 1764 20 1 0.996 0.005325 0.986 1 + 2772 19 1 0.994 0.007920 0.978 1 + 3444 18 1 0.991 0.010673 0.971 1 + 3542 17 1 0.988 0.013667 0.962 1 + 3780 16 1 0.985 0.016976 0.952 1 + 4860 15 1 0.981 0.020692 0.941 1 + 5196 14 1 0.977 0.024941 0.929 1 +> +> +> +> cleanEx() + +detaching ‘package:survival’ + +> nameEx("muscle") +> ### * muscle +> +> flush(stderr()); flush(stdout()) +> +> ### Name: muscle +> ### Title: Effect of Calcium Chloride on Muscle Contraction in Rat Hearts +> ### Aliases: muscle +> ### Keywords: datasets +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> A <- model.matrix(~ Strip - 1, data=muscle) +> rats.nls1 <- nls(log(Length) ~ cbind(A, rho^Conc), ++ data = muscle, start = c(rho=0.1), algorithm="plinear") +> (B <- coef(rats.nls1)) + rho .lin.StripS01 .lin.StripS02 .lin.StripS03 .lin.StripS04 + 0.07776401 3.08304824 3.30137838 3.44562531 2.80464434 +.lin.StripS05 .lin.StripS06 .lin.StripS07 .lin.StripS08 .lin.StripS09 + 2.60835015 3.03357725 3.52301734 3.38711844 3.46709396 +.lin.StripS10 .lin.StripS11 .lin.StripS12 .lin.StripS13 .lin.StripS14 + 3.81438456 3.73878664 3.51332581 3.39741115 3.47088608 +.lin.StripS15 .lin.StripS16 .lin.StripS17 .lin.StripS18 .lin.StripS19 + 3.72895847 3.31863862 3.37938673 2.96452195 3.58468686 +.lin.StripS20 .lin.StripS21 .lin22 + 3.39628029 3.36998872 -2.96015460 +> +> st <- list(alpha = B[2:22], beta = B[23], rho = B[1]) +> (rats.nls2 <- nls(log(Length) ~ alpha[Strip] + beta*rho^Conc, ++ data = muscle, start = st)) +Nonlinear regression model + model: log(Length) ~ alpha[Strip] + beta * rho^Conc + data: muscle +alpha..lin.StripS01 alpha..lin.StripS02 alpha..lin.StripS03 alpha..lin.StripS04 + 3.08305 3.30138 3.44563 2.80464 +alpha..lin.StripS05 alpha..lin.StripS06 alpha..lin.StripS07 alpha..lin.StripS08 + 2.60835 3.03358 3.52302 3.38712 +alpha..lin.StripS09 alpha..lin.StripS10 alpha..lin.StripS11 alpha..lin.StripS12 + 3.46709 3.81438 3.73879 3.51333 +alpha..lin.StripS13 alpha..lin.StripS14 alpha..lin.StripS15 alpha..lin.StripS16 + 3.39741 3.47089 3.72896 3.31864 +alpha..lin.StripS17 alpha..lin.StripS18 alpha..lin.StripS19 alpha..lin.StripS20 + 3.37939 2.96452 3.58469 3.39628 +alpha..lin.StripS21 beta..lin22 rho.rho + 3.36999 -2.96015 0.07776 + residual sum-of-squares: 1.045 + +Number of iterations to convergence: 0 +Achieved convergence tolerance: 4.918e-06 +> ## IGNORE_RDIFF_END +> +> Muscle <- with(muscle, { ++ Muscle <- expand.grid(Conc = sort(unique(Conc)), Strip = levels(Strip)) ++ Muscle$Yhat <- predict(rats.nls2, Muscle) ++ Muscle <- cbind(Muscle, logLength = rep(as.numeric(NA), 126)) ++ ind <- match(paste(Strip, Conc), ++ paste(Muscle$Strip, Muscle$Conc)) ++ Muscle$logLength[ind] <- log(Length) ++ Muscle}) +> +> lattice::xyplot(Yhat ~ Conc | Strip, Muscle, as.table = TRUE, ++ ylim = range(c(Muscle$Yhat, Muscle$logLength), na.rm = TRUE), ++ subscripts = TRUE, xlab = "Calcium Chloride concentration (mM)", ++ ylab = "log(Length in mm)", panel = ++ function(x, y, subscripts, ...) { ++ panel.xyplot(x, Muscle$logLength[subscripts], ...) ++ llines(spline(x, y)) ++ }) +> +> +> +> cleanEx() +> nameEx("mvrnorm") +> ### * mvrnorm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: mvrnorm +> ### Title: Simulate from a Multivariate Normal Distribution +> ### Aliases: mvrnorm +> ### Keywords: distribution multivariate +> +> ### ** Examples +> +> Sigma <- matrix(c(10,3,3,2),2,2) +> Sigma + [,1] [,2] +[1,] 10 3 +[2,] 3 2 +> var(mvrnorm(n = 1000, rep(0, 2), Sigma)) + [,1] [,2] +[1,] 10.697849 3.228279 +[2,] 3.228279 2.165271 +> var(mvrnorm(n = 1000, rep(0, 2), Sigma, empirical = TRUE)) + [,1] [,2] +[1,] 10 3 +[2,] 3 2 +> +> +> +> cleanEx() +> nameEx("negative.binomial") +> ### * negative.binomial +> +> flush(stderr()); flush(stdout()) +> +> ### Name: negative.binomial +> ### Title: Family function for Negative Binomial GLMs +> ### Aliases: negative.binomial +> ### Keywords: regression models +> +> ### ** Examples +> +> # Fitting a Negative Binomial model to the quine data +> # with theta = 2 assumed known. +> # +> glm(Days ~ .^4, family = negative.binomial(2), data = quine) + +Call: glm(formula = Days ~ .^4, family = negative.binomial(2), data = quine) + +Coefficients: + (Intercept) EthN SexM + 3.0564 -0.1386 -0.4914 + AgeF1 AgeF2 AgeF3 + -0.6227 -2.3632 -0.3784 + LrnSL EthN:SexM EthN:AgeF1 + -1.9577 -0.7524 0.1029 + EthN:AgeF2 EthN:AgeF3 EthN:LrnSL + -0.5546 0.0633 2.2588 + SexM:AgeF1 SexM:AgeF2 SexM:AgeF3 + 0.4092 3.1098 1.1145 + SexM:LrnSL AgeF1:LrnSL AgeF2:LrnSL + 1.5900 2.6421 4.8585 + AgeF3:LrnSL EthN:SexM:AgeF1 EthN:SexM:AgeF2 + NA -0.3105 0.3469 + EthN:SexM:AgeF3 EthN:SexM:LrnSL EthN:AgeF1:LrnSL + 0.8329 -0.1639 -3.5493 + EthN:AgeF2:LrnSL EthN:AgeF3:LrnSL SexM:AgeF1:LrnSL + -3.3315 NA -2.4285 + SexM:AgeF2:LrnSL SexM:AgeF3:LrnSL EthN:SexM:AgeF1:LrnSL + -4.1914 NA 2.1711 +EthN:SexM:AgeF2:LrnSL EthN:SexM:AgeF3:LrnSL + 2.1029 NA + +Degrees of Freedom: 145 Total (i.e. Null); 118 Residual +Null Deviance: 280.2 +Residual Deviance: 172 AIC: 1095 +> +> +> +> cleanEx() +> nameEx("nlschools") +> ### * nlschools +> +> flush(stderr()); flush(stdout()) +> +> ### Name: nlschools +> ### Title: Eighth-Grade Pupils in the Netherlands +> ### Aliases: nlschools +> ### Keywords: datasets +> +> ### ** Examples +> +> ## Don't show: +> op <- options(digits=5) +> ## End(Don't show) +> nl1 <- within(nlschools, { ++ IQave <- tapply(IQ, class, mean)[as.character(class)] ++ IQ <- IQ - IQave ++ }) +> cen <- c("IQ", "IQave", "SES") +> nl1[cen] <- scale(nl1[cen], center = TRUE, scale = FALSE) +> +> nl.lme <- nlme::lme(lang ~ IQ*COMB + IQave + SES, ++ random = ~ IQ | class, data = nl1) +> ## IGNORE_RDIFF_BEGIN +> summary(nl.lme) +Linear mixed-effects model fit by REML + Data: nl1 + AIC BIC logLik + 15120 15178 -7550.2 + +Random effects: + Formula: ~IQ | class + Structure: General positive-definite, Log-Cholesky parametrization + StdDev Corr +(Intercept) 2.78707 (Intr) +IQ 0.48424 -0.516 +Residual 6.24839 + +Fixed effects: lang ~ IQ * COMB + IQave + SES + Value Std.Error DF t-value p-value +(Intercept) 41.370 0.35364 2151 116.985 0.0000 +IQ 2.124 0.10070 2151 21.088 0.0000 +COMB1 -1.672 0.58719 130 -2.847 0.0051 +IQave 3.248 0.30021 130 10.818 0.0000 +SES 0.157 0.01465 2151 10.697 0.0000 +IQ:COMB1 0.431 0.18594 2151 2.317 0.0206 + Correlation: + (Intr) IQ COMB1 IQave SES +IQ -0.257 +COMB1 -0.609 0.155 +IQave -0.049 0.041 0.171 +SES 0.010 -0.190 -0.001 -0.168 +IQ:COMB1 0.139 -0.522 -0.206 -0.016 -0.003 + +Standardized Within-Group Residuals: + Min Q1 Med Q3 Max +-4.059387 -0.631084 0.065519 0.717864 2.794540 + +Number of Observations: 2287 +Number of Groups: 133 +> ## IGNORE_RDIFF_END +> ## Don't show: +> options(op) +> ## End(Don't show) +> +> +> +> cleanEx() +> nameEx("npk") +> ### * npk +> +> flush(stderr()); flush(stdout()) +> +> ### Name: npk +> ### Title: Classical N, P, K Factorial Experiment +> ### Aliases: npk +> ### Keywords: datasets +> +> ### ** Examples +> +> options(contrasts = c("contr.sum", "contr.poly")) +> npk.aov <- aov(yield ~ block + N*P*K, npk) +> ## IGNORE_RDIFF_BEGIN +> npk.aov +Call: + aov(formula = yield ~ block + N * P * K, data = npk) + +Terms: + block N P K N:P N:K P:K +Sum of Squares 343.2950 189.2817 8.4017 95.2017 21.2817 33.1350 0.4817 +Deg. of Freedom 5 1 1 1 1 1 1 + Residuals +Sum of Squares 185.2867 +Deg. of Freedom 12 + +Residual standard error: 3.929447 +1 out of 13 effects not estimable +Estimated effects may be unbalanced +> summary(npk.aov) + Df Sum Sq Mean Sq F value Pr(>F) +block 5 343.3 68.66 4.447 0.01594 * +N 1 189.3 189.28 12.259 0.00437 ** +P 1 8.4 8.40 0.544 0.47490 +K 1 95.2 95.20 6.166 0.02880 * +N:P 1 21.3 21.28 1.378 0.26317 +N:K 1 33.1 33.14 2.146 0.16865 +P:K 1 0.5 0.48 0.031 0.86275 +Residuals 12 185.3 15.44 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> alias(npk.aov) +Model : +yield ~ block + N * P * K + +Complete : + (Intercept) block1 block2 block3 block4 block5 N1 P1 K1 N1:P1 N1:K1 +N1:P1:K1 0 1 -1 -1 -1 1 0 0 0 0 0 + P1:K1 +N1:P1:K1 0 + +> coef(npk.aov) +(Intercept) block1 block2 block3 block4 block5 + 54.8750000 -0.8500000 2.5750000 5.9000000 -4.7500000 -4.3500000 + N1 P1 K1 N1:P1 N1:K1 P1:K1 + -2.8083333 0.5916667 1.9916667 -0.9416667 -1.1750000 0.1416667 +> options(contrasts = c("contr.treatment", "contr.poly")) +> npk.aov1 <- aov(yield ~ block + N + K, data = npk) +> summary.lm(npk.aov1) + +Call: +aov(formula = yield ~ block + N + K, data = npk) + +Residuals: + Min 1Q Median 3Q Max +-6.4083 -2.1438 0.2042 2.3292 7.0750 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +(Intercept) 53.208 2.276 23.381 8.5e-14 *** +block2 3.425 2.787 1.229 0.23690 +block3 6.750 2.787 2.422 0.02769 * +block4 -3.900 2.787 -1.399 0.18082 +block5 -3.500 2.787 -1.256 0.22723 +block6 2.325 2.787 0.834 0.41646 +N1 5.617 1.609 3.490 0.00302 ** +K1 -3.983 1.609 -2.475 0.02487 * +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 3.942 on 16 degrees of freedom +Multiple R-squared: 0.7163, Adjusted R-squared: 0.5922 +F-statistic: 5.772 on 7 and 16 DF, p-value: 0.001805 + +> se.contrast(npk.aov1, list(N=="0", N=="1"), data = npk) +[1] 1.609175 +> model.tables(npk.aov1, type = "means", se = TRUE) +Tables of means +Grand mean + +54.875 + + block +block + 1 2 3 4 5 6 +54.03 57.45 60.78 50.12 50.52 56.35 + + N +N + 0 1 +52.07 57.68 + + K +K + 0 1 +56.87 52.88 + +Standard errors for differences of means + block N K + 2.787 1.609 1.609 +replic. 4 12 12 +> ## IGNORE_RDIFF_END +> +> +> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly")) +> cleanEx() +> nameEx("oats") +> ### * oats +> +> flush(stderr()); flush(stdout()) +> +> ### Name: oats +> ### Title: Data from an Oats Field Trial +> ### Aliases: oats +> ### Keywords: datasets +> +> ### ** Examples +> +> oats$Nf <- ordered(oats$N, levels = sort(levels(oats$N))) +> oats.aov <- aov(Y ~ Nf*V + Error(B/V), data = oats, qr = TRUE) +> ## IGNORE_RDIFF_BEGIN +> summary(oats.aov) + +Error: B + Df Sum Sq Mean Sq F value Pr(>F) +Residuals 5 15875 3175 + +Error: B:V + Df Sum Sq Mean Sq F value Pr(>F) +V 2 1786 893.2 1.485 0.272 +Residuals 10 6013 601.3 + +Error: Within + Df Sum Sq Mean Sq F value Pr(>F) +Nf 3 20021 6674 37.686 2.46e-12 *** +Nf:V 6 322 54 0.303 0.932 +Residuals 45 7969 177 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> summary(oats.aov, split = list(Nf=list(L=1, Dev=2:3))) + +Error: B + Df Sum Sq Mean Sq F value Pr(>F) +Residuals 5 15875 3175 + +Error: B:V + Df Sum Sq Mean Sq F value Pr(>F) +V 2 1786 893.2 1.485 0.272 +Residuals 10 6013 601.3 + +Error: Within + Df Sum Sq Mean Sq F value Pr(>F) +Nf 3 20021 6674 37.686 2.46e-12 *** + Nf: L 1 19536 19536 110.323 1.09e-13 *** + Nf: Dev 2 484 242 1.367 0.265 +Nf:V 6 322 54 0.303 0.932 + Nf:V: L 2 168 84 0.475 0.625 + Nf:V: Dev 4 153 38 0.217 0.928 +Residuals 45 7969 177 +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> ## IGNORE_RDIFF_END +> par(mfrow = c(1,2), pty = "s") +> plot(fitted(oats.aov[[4]]), studres(oats.aov[[4]])) +> abline(h = 0, lty = 2) +> oats.pr <- proj(oats.aov) +> qqnorm(oats.pr[[4]][,"Residuals"], ylab = "Stratum 4 residuals") +> qqline(oats.pr[[4]][,"Residuals"]) +> +> par(mfrow = c(1,1), pty = "m") +> oats.aov2 <- aov(Y ~ N + V + Error(B/V), data = oats, qr = TRUE) +> model.tables(oats.aov2, type = "means", se = TRUE) +Warning in model.tables.aovlist(oats.aov2, type = "means", se = TRUE) : + SEs for type 'means' are not yet implemented +Calls: model.tables -> model.tables.aovlist +Tables of means +Grand mean + +103.9722 + + N +N +0.0cwt 0.2cwt 0.4cwt 0.6cwt + 79.39 98.89 114.22 123.39 + + V +V +Golden.rain Marvellous Victory + 104.50 109.79 97.63 +> +> +> +> graphics::par(get("par.postscript", pos = 'CheckExEnv')) +> cleanEx() +> nameEx("parcoord") +> ### * parcoord +> +> flush(stderr()); flush(stdout()) +> +> ### Name: parcoord +> ### Title: Parallel Coordinates Plot +> ### Aliases: parcoord +> ### Keywords: hplot +> +> ### ** Examples +> +> parcoord(state.x77[, c(7, 4, 6, 2, 5, 3)]) +> +> ir <- rbind(iris3[,,1], iris3[,,2], iris3[,,3]) +> parcoord(log(ir)[, c(3, 4, 2, 1)], col = 1 + (0:149)%/%50) +> +> +> +> cleanEx() +> nameEx("petrol") +> ### * petrol +> +> flush(stderr()); flush(stdout()) +> +> ### Name: petrol +> ### Title: N. L. Prater's Petrol Refinery Data +> ### Aliases: petrol +> ### Keywords: datasets +> +> ### ** Examples +> +> library(nlme) +> Petrol <- petrol +> Petrol[, 2:5] <- scale(as.matrix(Petrol[, 2:5]), scale = FALSE) +> pet3.lme <- lme(Y ~ SG + VP + V10 + EP, ++ random = ~ 1 | No, data = Petrol) +> pet3.lme <- update(pet3.lme, method = "ML") +> pet4.lme <- update(pet3.lme, fixed. = Y ~ V10 + EP) +> anova(pet4.lme, pet3.lme) + Model df AIC BIC logLik Test L.Ratio p-value +pet4.lme 1 5 149.6119 156.9406 -69.80594 +pet3.lme 2 7 149.3833 159.6435 -67.69166 1 vs 2 4.22855 0.1207 +> +> +> +> cleanEx() + +detaching ‘package:nlme’ + +> nameEx("plot.mca") +> ### * plot.mca +> +> flush(stderr()); flush(stdout()) +> +> ### Name: plot.mca +> ### Title: Plot Method for Objects of Class 'mca' +> ### Aliases: plot.mca +> ### Keywords: hplot multivariate +> +> ### ** Examples +> +> plot(mca(farms, abbrev = TRUE)) +> +> +> +> cleanEx() +> nameEx("plot.profile") +> ### * plot.profile +> +> flush(stderr()); flush(stdout()) +> +> ### Name: plot.profile +> ### Title: Plotting Functions for 'profile' Objects +> ### Aliases: plot.profile pairs.profile +> ### Keywords: models hplot +> +> ### ** Examples +> +> ## see ?profile.glm for an example using glm fits. +> +> ## a version of example(profile.nls) from R >= 2.8.0 +> fm1 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD) +> pr1 <- profile(fm1, alphamax = 0.1) +> MASS:::plot.profile(pr1) +> pairs(pr1) # a little odd since the parameters are highly correlated +> +> ## an example from ?nls +> x <- -(1:100)/10 +> y <- 100 + 10 * exp(x / 2) + rnorm(x)/10 +> nlmod <- nls(y ~ Const + A * exp(B * x), start=list(Const=100, A=10, B=1)) +> pairs(profile(nlmod)) +> +> +> +> cleanEx() +> nameEx("polr") +> ### * polr +> +> flush(stderr()); flush(stdout()) +> +> ### Name: polr +> ### Title: Ordered Logistic or Probit Regression +> ### Aliases: polr +> ### Keywords: models +> +> ### ** Examples +> +> options(contrasts = c("contr.treatment", "contr.poly")) +> house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) +> house.plr +Call: +polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq) + +Coefficients: + InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace + 0.5663937 1.2888191 -0.5723501 -0.3661866 -1.0910149 + ContHigh + 0.3602841 + +Intercepts: + Low|Medium Medium|High + -0.4961353 0.6907083 + +Residual Deviance: 3479.149 +AIC: 3495.149 +> summary(house.plr, digits = 3) + +Re-fitting to get Hessian + +Call: +polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq) + +Coefficients: + Value Std. Error t value +InflMedium 0.566 0.1047 5.41 +InflHigh 1.289 0.1272 10.14 +TypeApartment -0.572 0.1192 -4.80 +TypeAtrium -0.366 0.1552 -2.36 +TypeTerrace -1.091 0.1515 -7.20 +ContHigh 0.360 0.0955 3.77 + +Intercepts: + Value Std. Error t value +Low|Medium -0.496 0.125 -3.974 +Medium|High 0.691 0.125 5.505 + +Residual Deviance: 3479.149 +AIC: 3495.149 +> ## slightly worse fit from +> summary(update(house.plr, method = "probit", Hess = TRUE), digits = 3) +Call: +polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq, + Hess = TRUE, method = "probit") + +Coefficients: + Value Std. Error t value +InflMedium 0.346 0.0641 5.40 +InflHigh 0.783 0.0764 10.24 +TypeApartment -0.348 0.0723 -4.81 +TypeAtrium -0.218 0.0948 -2.30 +TypeTerrace -0.664 0.0918 -7.24 +ContHigh 0.222 0.0581 3.83 + +Intercepts: + Value Std. Error t value +Low|Medium -0.300 0.076 -3.937 +Medium|High 0.427 0.076 5.585 + +Residual Deviance: 3479.689 +AIC: 3495.689 +> ## although it is not really appropriate, can fit +> summary(update(house.plr, method = "loglog", Hess = TRUE), digits = 3) +Call: +polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq, + Hess = TRUE, method = "loglog") + +Coefficients: + Value Std. Error t value +InflMedium 0.367 0.0727 5.05 +InflHigh 0.790 0.0806 9.81 +TypeApartment -0.349 0.0757 -4.61 +TypeAtrium -0.196 0.0988 -1.98 +TypeTerrace -0.698 0.1043 -6.69 +ContHigh 0.268 0.0636 4.21 + +Intercepts: + Value Std. Error t value +Low|Medium 0.086 0.083 1.038 +Medium|High 0.892 0.087 10.223 + +Residual Deviance: 3491.41 +AIC: 3507.41 +> summary(update(house.plr, method = "cloglog", Hess = TRUE), digits = 3) +Call: +polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq, + Hess = TRUE, method = "cloglog") + +Coefficients: + Value Std. Error t value +InflMedium 0.382 0.0703 5.44 +InflHigh 0.915 0.0926 9.89 +TypeApartment -0.407 0.0861 -4.73 +TypeAtrium -0.281 0.1111 -2.52 +TypeTerrace -0.742 0.1013 -7.33 +ContHigh 0.209 0.0651 3.21 + +Intercepts: + Value Std. Error t value +Low|Medium -0.796 0.090 -8.881 +Medium|High 0.055 0.086 0.647 + +Residual Deviance: 3484.053 +AIC: 3500.053 +> +> predict(house.plr, housing, type = "p") + Low Medium High +1 0.3784493 0.2876752 0.3338755 +2 0.3784493 0.2876752 0.3338755 +3 0.3784493 0.2876752 0.3338755 +4 0.2568264 0.2742122 0.4689613 +5 0.2568264 0.2742122 0.4689613 +6 0.2568264 0.2742122 0.4689613 +7 0.1436924 0.2110836 0.6452240 +8 0.1436924 0.2110836 0.6452240 +9 0.1436924 0.2110836 0.6452240 +10 0.5190445 0.2605077 0.2204478 +11 0.5190445 0.2605077 0.2204478 +12 0.5190445 0.2605077 0.2204478 +13 0.3798514 0.2875965 0.3325521 +14 0.3798514 0.2875965 0.3325521 +15 0.3798514 0.2875965 0.3325521 +16 0.2292406 0.2643196 0.5064398 +17 0.2292406 0.2643196 0.5064398 +18 0.2292406 0.2643196 0.5064398 +19 0.4675584 0.2745383 0.2579033 +20 0.4675584 0.2745383 0.2579033 +21 0.4675584 0.2745383 0.2579033 +22 0.3326236 0.2876008 0.3797755 +23 0.3326236 0.2876008 0.3797755 +24 0.3326236 0.2876008 0.3797755 +25 0.1948548 0.2474226 0.5577225 +26 0.1948548 0.2474226 0.5577225 +27 0.1948548 0.2474226 0.5577225 +28 0.6444840 0.2114256 0.1440905 +29 0.6444840 0.2114256 0.1440905 +30 0.6444840 0.2114256 0.1440905 +31 0.5071210 0.2641196 0.2287594 +32 0.5071210 0.2641196 0.2287594 +33 0.5071210 0.2641196 0.2287594 +34 0.3331573 0.2876330 0.3792097 +35 0.3331573 0.2876330 0.3792097 +36 0.3331573 0.2876330 0.3792097 +37 0.2980880 0.2837746 0.4181374 +38 0.2980880 0.2837746 0.4181374 +39 0.2980880 0.2837746 0.4181374 +40 0.1942209 0.2470589 0.5587202 +41 0.1942209 0.2470589 0.5587202 +42 0.1942209 0.2470589 0.5587202 +43 0.1047770 0.1724227 0.7228003 +44 0.1047770 0.1724227 0.7228003 +45 0.1047770 0.1724227 0.7228003 +46 0.4294564 0.2820629 0.2884807 +47 0.4294564 0.2820629 0.2884807 +48 0.4294564 0.2820629 0.2884807 +49 0.2993357 0.2839753 0.4166890 +50 0.2993357 0.2839753 0.4166890 +51 0.2993357 0.2839753 0.4166890 +52 0.1718050 0.2328648 0.5953302 +53 0.1718050 0.2328648 0.5953302 +54 0.1718050 0.2328648 0.5953302 +55 0.3798387 0.2875972 0.3325641 +56 0.3798387 0.2875972 0.3325641 +57 0.3798387 0.2875972 0.3325641 +58 0.2579546 0.2745537 0.4674917 +59 0.2579546 0.2745537 0.4674917 +60 0.2579546 0.2745537 0.4674917 +61 0.1444202 0.2117081 0.6438717 +62 0.1444202 0.2117081 0.6438717 +63 0.1444202 0.2117081 0.6438717 +64 0.5583813 0.2471826 0.1944361 +65 0.5583813 0.2471826 0.1944361 +66 0.5583813 0.2471826 0.1944361 +67 0.4178031 0.2838213 0.2983756 +68 0.4178031 0.2838213 0.2983756 +69 0.4178031 0.2838213 0.2983756 +70 0.2584149 0.2746916 0.4668935 +71 0.2584149 0.2746916 0.4668935 +72 0.2584149 0.2746916 0.4668935 +> addterm(house.plr, ~.^2, test = "Chisq") +Single term additions + +Model: +Sat ~ Infl + Type + Cont + Df AIC LRT Pr(Chi) + 3495.1 +Infl:Type 6 3484.6 22.5093 0.0009786 *** +Infl:Cont 2 3498.9 0.2090 0.9007957 +Type:Cont 3 3492.5 8.6662 0.0340752 * +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> house.plr2 <- stepAIC(house.plr, ~.^2) +Start: AIC=3495.15 +Sat ~ Infl + Type + Cont + + Df AIC ++ Infl:Type 6 3484.6 ++ Type:Cont 3 3492.5 + 3495.1 ++ Infl:Cont 2 3498.9 +- Cont 1 3507.5 +- Type 3 3545.1 +- Infl 2 3599.4 + +Step: AIC=3484.64 +Sat ~ Infl + Type + Cont + Infl:Type + + Df AIC ++ Type:Cont 3 3482.7 + 3484.6 ++ Infl:Cont 2 3488.5 +- Infl:Type 6 3495.1 +- Cont 1 3497.8 + +Step: AIC=3482.69 +Sat ~ Infl + Type + Cont + Infl:Type + Type:Cont + + Df AIC + 3482.7 +- Type:Cont 3 3484.6 ++ Infl:Cont 2 3486.6 +- Infl:Type 6 3492.5 +> house.plr2$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +Sat ~ Infl + Type + Cont + +Final Model: +Sat ~ Infl + Type + Cont + Infl:Type + Type:Cont + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 1673 3479.149 3495.149 +2 + Infl:Type 6 22.509347 1667 3456.640 3484.640 +3 + Type:Cont 3 7.945029 1664 3448.695 3482.695 +> anova(house.plr, house.plr2) +Likelihood ratio tests of ordinal regression models + +Response: Sat + Model Resid. df Resid. Dev Test Df +1 Infl + Type + Cont 1673 3479.149 +2 Infl + Type + Cont + Infl:Type + Type:Cont 1664 3448.695 1 vs 2 9 + LR stat. Pr(Chi) +1 +2 30.45438 0.0003670555 +> +> house.plr <- update(house.plr, Hess=TRUE) +> pr <- profile(house.plr) +> confint(pr) + 2.5 % 97.5 % +InflMedium 0.3616415 0.77195375 +InflHigh 1.0409701 1.53958138 +TypeApartment -0.8069590 -0.33940432 +TypeAtrium -0.6705862 -0.06204495 +TypeTerrace -1.3893863 -0.79533958 +ContHigh 0.1733589 0.54792854 +> plot(pr) +> pairs(pr) +> +> +> +> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly")) +> cleanEx() +> nameEx("predict.glmmPQL") +> ### * predict.glmmPQL +> +> flush(stderr()); flush(stdout()) +> +> ### Name: predict.glmmPQL +> ### Title: Predict Method for glmmPQL Fits +> ### Aliases: predict.glmmPQL +> ### Keywords: models +> +> ### ** Examples +> +> fit <- glmmPQL(y ~ trt + I(week > 2), random = ~1 | ID, ++ family = binomial, data = bacteria) +iteration 1 +iteration 2 +iteration 3 +iteration 4 +iteration 5 +iteration 6 +> predict(fit, bacteria, level = 0, type="response") + [1] 0.9680779 0.9680779 0.8587270 0.8587270 0.9344832 0.9344832 0.7408574 + [8] 0.7408574 0.8970307 0.8970307 0.6358511 0.6358511 0.6358511 0.9680779 + [15] 0.9680779 0.8587270 0.8587270 0.8587270 0.9680779 0.9680779 0.8587270 + [22] 0.8587270 0.8587270 0.8970307 0.8970307 0.6358511 0.6358511 0.9344832 + [29] 0.9344832 0.7408574 0.7408574 0.7408574 0.9680779 0.9680779 0.8587270 + [36] 0.8587270 0.8587270 0.9680779 0.9680779 0.8587270 0.8587270 0.8587270 + [43] 0.9344832 0.7408574 0.9680779 0.9680779 0.8587270 0.8587270 0.8587270 + [50] 0.8970307 0.8970307 0.6358511 0.6358511 0.6358511 0.9680779 0.9680779 + [57] 0.8587270 0.8587270 0.8587270 0.9680779 0.9680779 0.8587270 0.8970307 + [64] 0.8970307 0.6358511 0.6358511 0.6358511 0.9344832 0.9344832 0.7408574 + [71] 0.7408574 0.7408574 0.9680779 0.9680779 0.8587270 0.8587270 0.8587270 + [78] 0.8970307 0.8970307 0.6358511 0.6358511 0.6358511 0.9680779 0.9680779 + [85] 0.8587270 0.8587270 0.8587270 0.9344832 0.9344832 0.7408574 0.7408574 + [92] 0.9680779 0.9680779 0.8587270 0.8587270 0.8587270 0.9680779 0.9680779 + [99] 0.8587270 0.8587270 0.8587270 0.9680779 0.9680779 0.8587270 0.8587270 +[106] 0.8587270 0.9344832 0.9344832 0.7408574 0.7408574 0.7408574 0.8970307 +[113] 0.8970307 0.6358511 0.6358511 0.9680779 0.9680779 0.8587270 0.9680779 +[120] 0.9680779 0.8587270 0.8587270 0.8970307 0.8970307 0.6358511 0.6358511 +[127] 0.6358511 0.9344832 0.7408574 0.7408574 0.7408574 0.9680779 0.8587270 +[134] 0.8587270 0.8587270 0.8970307 0.8970307 0.6358511 0.6358511 0.6358511 +[141] 0.9680779 0.9680779 0.8587270 0.8587270 0.8587270 0.9344832 0.7408574 +[148] 0.8970307 0.8970307 0.6358511 0.6358511 0.9680779 0.9680779 0.8587270 +[155] 0.8970307 0.8970307 0.6358511 0.9680779 0.9680779 0.8587270 0.8587270 +[162] 0.8587270 0.9344832 0.9344832 0.7408574 0.7408574 0.7408574 0.9680779 +[169] 0.9680779 0.8587270 0.8587270 0.8587270 0.9344832 0.7408574 0.8970307 +[176] 0.8970307 0.6358511 0.6358511 0.6358511 0.9344832 0.9344832 0.7408574 +[183] 0.7408574 0.9680779 0.9680779 0.8587270 0.8587270 0.8587270 0.8970307 +[190] 0.8970307 0.6358511 0.6358511 0.6358511 0.9344832 0.9344832 0.7408574 +[197] 0.7408574 0.7408574 0.8970307 0.6358511 0.6358511 0.9344832 0.9344832 +[204] 0.7408574 0.7408574 0.7408574 0.8970307 0.8970307 0.6358511 0.6358511 +[211] 0.9344832 0.9344832 0.7408574 0.7408574 0.7408574 0.9344832 0.9344832 +[218] 0.7408574 0.7408574 0.7408574 +attr(,"label") +[1] "Predicted values" +> predict(fit, bacteria, level = 1, type="response") + X01 X01 X01 X01 X02 X02 X02 X02 +0.9828449 0.9828449 0.9198935 0.9198935 0.9050782 0.9050782 0.6564944 0.6564944 + X03 X03 X03 X03 X03 X04 X04 X04 +0.9724022 0.9724022 0.8759665 0.8759665 0.8759665 0.9851548 0.9851548 0.9300763 + X04 X04 X05 X05 X05 X05 X05 X06 +0.9300763 0.9300763 0.9851548 0.9851548 0.9300763 0.9300763 0.9300763 0.9662755 + X06 X06 X06 X07 X07 X07 X07 X07 +0.9662755 0.8516962 0.8516962 0.7291679 0.7291679 0.3504978 0.3504978 0.3504978 + X08 X08 X08 X08 X08 X09 X09 X09 +0.9426815 0.9426815 0.7672499 0.7672499 0.7672499 0.9851548 0.9851548 0.9300763 + X09 X09 X10 X10 X11 X11 X11 X11 +0.9300763 0.9300763 0.9640326 0.8430706 0.9851548 0.9851548 0.9300763 0.9300763 + X11 X12 X12 X12 X12 X12 X13 X13 +0.9300763 0.8334870 0.8334870 0.5008219 0.5008219 0.5008219 0.9851548 0.9851548 + X13 X13 X13 X14 X14 X14 X15 X15 +0.9300763 0.9300763 0.9300763 0.8907227 0.8907227 0.6203155 0.9724022 0.9724022 + X15 X15 X15 X16 X16 X16 X16 X16 +0.8759665 0.8759665 0.8759665 0.9287777 0.9287777 0.7232833 0.7232833 0.7232833 + X17 X17 X17 X17 X17 X18 X18 X18 +0.9426815 0.9426815 0.7672499 0.7672499 0.7672499 0.7070916 0.7070916 0.3260827 + X18 X18 X19 X19 X19 X19 X19 X20 +0.3260827 0.3260827 0.8702991 0.8702991 0.5735499 0.5735499 0.5735499 0.9736293 + X20 X20 X20 X21 X21 X21 X21 X21 +0.9736293 0.8809564 0.8809564 0.9851548 0.9851548 0.9300763 0.9300763 0.9300763 + Y01 Y01 Y01 Y01 Y01 Y02 Y02 Y02 +0.9851548 0.9851548 0.9300763 0.9300763 0.9300763 0.7607971 0.7607971 0.3893126 + Y02 Y02 Y03 Y03 Y03 Y03 Y03 Y04 +0.3893126 0.3893126 0.8487181 0.8487181 0.5292976 0.5292976 0.5292976 0.5734482 + Y04 Y04 Y04 Y05 Y05 Y05 Y06 Y06 +0.5734482 0.2122655 0.2122655 0.7144523 0.7144523 0.3339997 0.9828449 0.9828449 + Y06 Y06 Y07 Y07 Y07 Y07 Y07 Y08 +0.9198935 0.9198935 0.8334870 0.8334870 0.5008219 0.5008219 0.5008219 0.9238389 + Y08 Y08 Y08 Y09 Y09 Y09 Y09 Y10 +0.7085660 0.7085660 0.7085660 0.9847299 0.9281899 0.9281899 0.9281899 0.9188296 + Y10 Y10 Y10 Y10 Y11 Y11 Y11 Y11 +0.9188296 0.6940862 0.6940862 0.6940862 0.9851548 0.9851548 0.9300763 0.9300763 + Y11 Y12 Y12 Y13 Y13 Y13 Y13 Y14 +0.9300763 0.9640326 0.8430706 0.5734482 0.5734482 0.2122655 0.2122655 0.9793383 + Y14 Y14 Z01 Z01 Z01 Z02 Z02 Z02 +0.9793383 0.9047659 0.9556329 0.9556329 0.8119328 0.9851548 0.9851548 0.9300763 + Z02 Z02 Z03 Z03 Z03 Z03 Z03 Z05 +0.9300763 0.9300763 0.9779690 0.9779690 0.8989642 0.8989642 0.8989642 0.8702991 + Z05 Z05 Z05 Z05 Z06 Z06 Z07 Z07 +0.8702991 0.5735499 0.5735499 0.5735499 0.8306525 0.4957505 0.8334870 0.8334870 + Z07 Z07 Z07 Z09 Z09 Z09 Z09 Z10 +0.5008219 0.5008219 0.5008219 0.9736293 0.9736293 0.8809564 0.8809564 0.9851548 + Z10 Z10 Z10 Z10 Z11 Z11 Z11 Z11 +0.9851548 0.9300763 0.9300763 0.9300763 0.9724022 0.9724022 0.8759665 0.8759665 + Z11 Z14 Z14 Z14 Z14 Z14 Z15 Z15 +0.8759665 0.9287777 0.9287777 0.7232833 0.7232833 0.7232833 0.9643851 0.8444172 + Z15 Z19 Z19 Z19 Z19 Z19 Z20 Z20 +0.8444172 0.9779690 0.9779690 0.8989642 0.8989642 0.8989642 0.7620490 0.7620490 + Z20 Z20 Z24 Z24 Z24 Z24 Z24 Z26 +0.3909523 0.3909523 0.8487181 0.8487181 0.5292976 0.5292976 0.5292976 0.9287777 + Z26 Z26 Z26 Z26 +0.9287777 0.7232833 0.7232833 0.7232833 +attr(,"label") +[1] "Predicted values" +> +> +> +> cleanEx() +> nameEx("predict.lda") +> ### * predict.lda +> +> flush(stderr()); flush(stdout()) +> +> ### Name: predict.lda +> ### Title: Classify Multivariate Observations by Linear Discrimination +> ### Aliases: predict.lda +> ### Keywords: multivariate +> +> ### ** Examples +> +> tr <- sample(1:50, 25) +> train <- rbind(iris3[tr,,1], iris3[tr,,2], iris3[tr,,3]) +> test <- rbind(iris3[-tr,,1], iris3[-tr,,2], iris3[-tr,,3]) +> cl <- factor(c(rep("s",25), rep("c",25), rep("v",25))) +> z <- lda(train, cl) +> predict(z, test)$class + [1] s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c c c c c c c c c +[39] c c c c c c c c c c c c v v v v v v v v v v v v v v v v v c v v v v v v v +Levels: c s v +> +> +> +> cleanEx() +> nameEx("predict.lqs") +> ### * predict.lqs +> +> flush(stderr()); flush(stdout()) +> +> ### Name: predict.lqs +> ### Title: Predict from an lqs Fit +> ### Aliases: predict.lqs +> ### Keywords: models +> +> ### ** Examples +> +> set.seed(123) +> fm <- lqs(stack.loss ~ ., data = stackloss, method = "S", nsamp = "exact") +> predict(fm, stackloss) + 1 2 3 4 5 6 7 8 +35.500000 35.579646 30.409292 19.477876 18.592920 19.035398 19.000000 19.000000 + 9 10 11 12 13 14 15 16 +15.734513 14.079646 13.362832 13.000000 13.920354 13.486726 6.761062 7.000000 + 17 18 19 20 21 + 8.557522 8.000000 8.362832 13.154867 23.991150 +> +> +> +> cleanEx() +> nameEx("predict.qda") +> ### * predict.qda +> +> flush(stderr()); flush(stdout()) +> +> ### Name: predict.qda +> ### Title: Classify from Quadratic Discriminant Analysis +> ### Aliases: predict.qda +> ### Keywords: multivariate +> +> ### ** Examples +> +> tr <- sample(1:50, 25) +> train <- rbind(iris3[tr,,1], iris3[tr,,2], iris3[tr,,3]) +> test <- rbind(iris3[-tr,,1], iris3[-tr,,2], iris3[-tr,,3]) +> cl <- factor(c(rep("s",25), rep("c",25), rep("v",25))) +> zq <- qda(train, cl) +> predict(zq, test)$class + [1] s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c c c v c c c c c +[39] c c c c c c c c c c c c v v v v v v v v v v v v v v v v v v v v v v v v v +Levels: c s v +> +> +> +> cleanEx() +> nameEx("profile.glm") +> ### * profile.glm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: profile.glm +> ### Title: Method for Profiling glm Objects +> ### Aliases: profile.glm +> ### Keywords: regression models +> +> ### ** Examples +> +> options(contrasts = c("contr.treatment", "contr.poly")) +> ldose <- rep(0:5, 2) +> numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16) +> sex <- factor(rep(c("M", "F"), c(6, 6))) +> SF <- cbind(numdead, numalive = 20 - numdead) +> budworm.lg <- glm(SF ~ sex*ldose, family = binomial) +> pr1 <- profile(budworm.lg) +> plot(pr1) +> pairs(pr1) +> +> +> +> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly")) +> cleanEx() +> nameEx("qda") +> ### * qda +> +> flush(stderr()); flush(stdout()) +> +> ### Name: qda +> ### Title: Quadratic Discriminant Analysis +> ### Aliases: qda qda.data.frame qda.default qda.formula qda.matrix +> ### model.frame.qda print.qda +> ### Keywords: multivariate +> +> ### ** Examples +> +> tr <- sample(1:50, 25) +> train <- rbind(iris3[tr,,1], iris3[tr,,2], iris3[tr,,3]) +> test <- rbind(iris3[-tr,,1], iris3[-tr,,2], iris3[-tr,,3]) +> cl <- factor(c(rep("s",25), rep("c",25), rep("v",25))) +> z <- qda(train, cl) +> predict(z,test)$class + [1] s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c c c v c c c c c +[39] c c c c c c c c c c c c v v v v v v v v v v v v v v v v v v v v v v v v v +Levels: c s v +> +> +> +> cleanEx() +> nameEx("rational") +> ### * rational +> +> flush(stderr()); flush(stdout()) +> +> ### Name: rational +> ### Title: Rational Approximation +> ### Aliases: rational .rat +> ### Keywords: math +> +> ### ** Examples +> +> X <- matrix(runif(25), 5, 5) +> zapsmall(solve(X, X/5)) # print near-zeroes as zero + [,1] [,2] [,3] [,4] [,5] +[1,] 0.2 0.0 0.0 0.0 0.0 +[2,] 0.0 0.2 0.0 0.0 0.0 +[3,] 0.0 0.0 0.2 0.0 0.0 +[4,] 0.0 0.0 0.0 0.2 0.0 +[5,] 0.0 0.0 0.0 0.0 0.2 +> rational(solve(X, X/5)) + [,1] [,2] [,3] [,4] [,5] +[1,] 0.2 0.0 0.0 0.0 0.0 +[2,] 0.0 0.2 0.0 0.0 0.0 +[3,] 0.0 0.0 0.2 0.0 0.0 +[4,] 0.0 0.0 0.0 0.2 0.0 +[5,] 0.0 0.0 0.0 0.0 0.2 +> +> +> +> cleanEx() +> nameEx("renumerate") +> ### * renumerate +> +> flush(stderr()); flush(stdout()) +> +> ### Name: renumerate +> ### Title: Convert a Formula Transformed by 'denumerate' +> ### Aliases: renumerate renumerate.formula +> ### Keywords: models +> +> ### ** Examples +> +> denumerate(~(1+2+3)^3 + a/b) +~(.v1 + .v2 + .v3)^3 + a/b +> ## ~ (.v1 + .v2 + .v3)^3 + a/b +> renumerate(.Last.value) +~(`1` + `2` + `3`)^3 + a/b +> ## ~ (1 + 2 + 3)^3 + a/b +> +> +> +> cleanEx() +> nameEx("rlm") +> ### * rlm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: rlm +> ### Title: Robust Fitting of Linear Models +> ### Aliases: rlm rlm.default rlm.formula print.rlm predict.rlm psi.bisquare +> ### psi.hampel psi.huber +> ### Keywords: models robust +> +> ### ** Examples +> +> summary(rlm(stack.loss ~ ., stackloss)) + +Call: rlm(formula = stack.loss ~ ., data = stackloss) +Residuals: + Min 1Q Median 3Q Max +-8.91753 -1.73127 0.06187 1.54306 6.50163 + +Coefficients: + Value Std. Error t value +(Intercept) -41.0265 9.8073 -4.1832 +Air.Flow 0.8294 0.1112 7.4597 +Water.Temp 0.9261 0.3034 3.0524 +Acid.Conc. -0.1278 0.1289 -0.9922 + +Residual standard error: 2.441 on 17 degrees of freedom +> rlm(stack.loss ~ ., stackloss, psi = psi.hampel, init = "lts") +Call: +rlm(formula = stack.loss ~ ., data = stackloss, psi = psi.hampel, + init = "lts") +Converged in 9 iterations + +Coefficients: +(Intercept) Air.Flow Water.Temp Acid.Conc. +-40.4747826 0.7410853 1.2250730 -0.1455245 + +Degrees of freedom: 21 total; 17 residual +Scale estimate: 3.09 +> rlm(stack.loss ~ ., stackloss, psi = psi.bisquare) +Call: +rlm(formula = stack.loss ~ ., data = stackloss, psi = psi.bisquare) +Converged in 11 iterations + +Coefficients: +(Intercept) Air.Flow Water.Temp Acid.Conc. +-42.2852537 0.9275471 0.6507322 -0.1123310 + +Degrees of freedom: 21 total; 17 residual +Scale estimate: 2.28 +> +> +> +> cleanEx() +> nameEx("rms.curv") +> ### * rms.curv +> +> flush(stderr()); flush(stdout()) +> +> ### Name: rms.curv +> ### Title: Relative Curvature Measures for Non-Linear Regression +> ### Aliases: rms.curv print.rms.curv +> ### Keywords: nonlinear +> +> ### ** Examples +> +> # The treated sample from the Puromycin data +> mmcurve <- deriv3(~ Vm * conc/(K + conc), c("Vm", "K"), ++ function(Vm, K, conc) NULL) +> Treated <- Puromycin[Puromycin$state == "treated", ] +> (Purfit1 <- nls(rate ~ mmcurve(Vm, K, conc), data = Treated, ++ start = list(Vm=200, K=0.1))) +Nonlinear regression model + model: rate ~ mmcurve(Vm, K, conc) + data: Treated + Vm K +212.68363 0.06412 + residual sum-of-squares: 1195 + +Number of iterations to convergence: 6 +Achieved convergence tolerance: 6.096e-06 +> rms.curv(Purfit1) +Parameter effects: c^theta x sqrt(F) = 0.2121 + Intrinsic: c^iota x sqrt(F) = 0.092 +> ##Parameter effects: c^theta x sqrt(F) = 0.2121 +> ## Intrinsic: c^iota x sqrt(F) = 0.092 +> +> +> +> cleanEx() +> nameEx("rnegbin") +> ### * rnegbin +> +> flush(stderr()); flush(stdout()) +> +> ### Name: rnegbin +> ### Title: Simulate Negative Binomial Variates +> ### Aliases: rnegbin +> ### Keywords: distribution +> +> ### ** Examples +> +> # Negative Binomials with means fitted(fm) and theta = 4.5 +> fm <- glm.nb(Days ~ ., data = quine) +> dummy <- rnegbin(fitted(fm), theta = 4.5) +> +> +> +> cleanEx() +> nameEx("sammon") +> ### * sammon +> +> flush(stderr()); flush(stdout()) +> +> ### Name: sammon +> ### Title: Sammon's Non-Linear Mapping +> ### Aliases: sammon +> ### Keywords: multivariate +> +> ### ** Examples +> +> swiss.x <- as.matrix(swiss[, -1]) +> swiss.sam <- sammon(dist(swiss.x)) +Initial stress : 0.00824 +stress after 10 iters: 0.00439, magic = 0.338 +stress after 20 iters: 0.00383, magic = 0.500 +stress after 30 iters: 0.00383, magic = 0.500 +> plot(swiss.sam$points, type = "n") +> text(swiss.sam$points, labels = as.character(1:nrow(swiss.x))) +> +> +> +> cleanEx() +> nameEx("stepAIC") +> ### * stepAIC +> +> flush(stderr()); flush(stdout()) +> +> ### Name: stepAIC +> ### Title: Choose a model by AIC in a Stepwise Algorithm +> ### Aliases: stepAIC extractAIC.gls terms.gls extractAIC.lme terms.lme +> ### Keywords: models +> +> ### ** Examples +> +> quine.hi <- aov(log(Days + 2.5) ~ .^4, quine) +> quine.nxt <- update(quine.hi, . ~ . - Eth:Sex:Age:Lrn) +> quine.stp <- stepAIC(quine.nxt, ++ scope = list(upper = ~Eth*Sex*Age*Lrn, lower = ~1), ++ trace = FALSE) +> quine.stp$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Age + Eth:Sex:Lrn + + Eth:Age:Lrn + Sex:Age:Lrn + +Final Model: +log(Days + 2.5) ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 120 64.09900 -68.18396 +2 - Eth:Sex:Age 3 0.973869 123 65.07287 -71.98244 +3 - Sex:Age:Lrn 2 1.526754 125 66.59962 -72.59652 +> +> cpus1 <- cpus +> for(v in names(cpus)[2:7]) ++ cpus1[[v]] <- cut(cpus[[v]], unique(quantile(cpus[[v]])), ++ include.lowest = TRUE) +> cpus0 <- cpus1[, 2:8] # excludes names, authors' predictions +> cpus.samp <- sample(1:209, 100) +> cpus.lm <- lm(log10(perf) ~ ., data = cpus1[cpus.samp,2:8]) +> cpus.lm2 <- stepAIC(cpus.lm, trace = FALSE) +> cpus.lm2$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +log10(perf) ~ syct + mmin + mmax + cach + chmin + chmax + +Final Model: +log10(perf) ~ syct + mmax + cach + chmax + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 82 3.458189 -300.4425 +2 - chmin 3 0.02548983 85 3.483679 -305.7081 +3 - mmin 3 0.12039102 88 3.604070 -308.3106 +> +> example(birthwt) + +brthwt> bwt <- with(birthwt, { +brthwt+ race <- factor(race, labels = c("white", "black", "other")) +brthwt+ ptd <- factor(ptl > 0) +brthwt+ ftv <- factor(ftv) +brthwt+ levels(ftv)[-(1:2)] <- "2+" +brthwt+ data.frame(low = factor(low), age, lwt, race, smoke = (smoke > 0), +brthwt+ ptd, ht = (ht > 0), ui = (ui > 0), ftv) +brthwt+ }) + +brthwt> options(contrasts = c("contr.treatment", "contr.poly")) + +brthwt> glm(low ~ ., binomial, bwt) + +Call: glm(formula = low ~ ., family = binomial, data = bwt) + +Coefficients: +(Intercept) age lwt raceblack raceother smokeTRUE + 0.82302 -0.03723 -0.01565 1.19241 0.74068 0.75553 + ptdTRUE htTRUE uiTRUE ftv1 ftv2+ + 1.34376 1.91317 0.68020 -0.43638 0.17901 + +Degrees of Freedom: 188 Total (i.e. Null); 178 Residual +Null Deviance: 234.7 +Residual Deviance: 195.5 AIC: 217.5 +> birthwt.glm <- glm(low ~ ., family = binomial, data = bwt) +> birthwt.step <- stepAIC(birthwt.glm, trace = FALSE) +> birthwt.step$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +low ~ age + lwt + race + smoke + ptd + ht + ui + ftv + +Final Model: +low ~ lwt + race + smoke + ptd + ht + ui + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 178 195.4755 217.4755 +2 - ftv 2 1.358185 180 196.8337 214.8337 +3 - age 1 1.017866 181 197.8516 213.8516 +> birthwt.step2 <- stepAIC(birthwt.glm, ~ .^2 + I(scale(age)^2) ++ + I(scale(lwt)^2), trace = FALSE) +> birthwt.step2$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +low ~ age + lwt + race + smoke + ptd + ht + ui + ftv + +Final Model: +low ~ age + lwt + smoke + ptd + ht + ui + ftv + age:ftv + smoke:ui + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 178 195.4755 217.4755 +2 + age:ftv 2 12.474896 176 183.0006 209.0006 +3 + smoke:ui 1 3.056805 175 179.9438 207.9438 +4 - race 2 3.129586 177 183.0734 207.0734 +> +> quine.nb <- glm.nb(Days ~ .^4, data = quine) +> quine.nb2 <- stepAIC(quine.nb) +Start: AIC=1095.32 +Days ~ (Eth + Sex + Age + Lrn)^4 + + Df AIC +- Eth:Sex:Age:Lrn 2 1092.7 + 1095.3 + +Step: AIC=1092.73 +Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Age + Eth:Sex:Lrn + + Eth:Age:Lrn + Sex:Age:Lrn + + Df AIC +- Eth:Sex:Age 3 1089.4 + 1092.7 +- Eth:Sex:Lrn 1 1093.3 +- Eth:Age:Lrn 2 1094.7 +- Sex:Age:Lrn 2 1095.0 + +Step: AIC=1089.41 +Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + + Sex:Age:Lrn + + Df AIC + 1089.4 +- Sex:Age:Lrn 2 1091.1 +- Eth:Age:Lrn 2 1091.2 +- Eth:Sex:Lrn 1 1092.5 +> quine.nb2$anova +Stepwise Model Path +Analysis of Deviance Table + +Initial Model: +Days ~ (Eth + Sex + Age + Lrn)^4 + +Final Model: +Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + + Sex:Age:Lrn + + + Step Df Deviance Resid. Df Resid. Dev AIC +1 118 167.4535 1095.324 +2 - Eth:Sex:Age:Lrn 2 0.09746244 120 167.5509 1092.728 +3 - Eth:Sex:Age 3 0.11060087 123 167.4403 1089.409 +> +> +> +> cleanEx() +> nameEx("summary.negbin") +> ### * summary.negbin +> +> flush(stderr()); flush(stdout()) +> +> ### Name: summary.negbin +> ### Title: Summary Method Function for Objects of Class 'negbin' +> ### Aliases: summary.negbin print.summary.negbin +> ### Keywords: models +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> summary(glm.nb(Days ~ Eth*Age*Lrn*Sex, quine, link = log)) + +Call: +glm.nb(formula = Days ~ Eth * Age * Lrn * Sex, data = quine, + link = log, init.theta = 1.928360145) + +Coefficients: (4 not defined because of singularities) + Estimate Std. Error z value Pr(>|z|) +(Intercept) 3.0564 0.3760 8.128 4.38e-16 *** +EthN -0.1386 0.5334 -0.260 0.795023 +AgeF1 -0.6227 0.5125 -1.215 0.224334 +AgeF2 -2.3632 1.0770 -2.194 0.028221 * +AgeF3 -0.3784 0.4546 -0.832 0.405215 +LrnSL -1.9577 0.9967 -1.964 0.049493 * +SexM -0.4914 0.5104 -0.963 0.335653 +EthN:AgeF1 0.1029 0.7123 0.144 0.885175 +EthN:AgeF2 -0.5546 1.6798 -0.330 0.741297 +EthN:AgeF3 0.0633 0.6396 0.099 0.921159 +EthN:LrnSL 2.2588 1.3019 1.735 0.082743 . +AgeF1:LrnSL 2.6421 1.0821 2.442 0.014618 * +AgeF2:LrnSL 4.8585 1.4423 3.369 0.000755 *** +AgeF3:LrnSL NA NA NA NA +EthN:SexM -0.7524 0.7220 -1.042 0.297400 +AgeF1:SexM 0.4092 0.8299 0.493 0.621973 +AgeF2:SexM 3.1098 1.1655 2.668 0.007624 ** +AgeF3:SexM 1.1145 0.6365 1.751 0.079926 . +LrnSL:SexM 1.5900 1.1499 1.383 0.166750 +EthN:AgeF1:LrnSL -3.5493 1.4270 -2.487 0.012876 * +EthN:AgeF2:LrnSL -3.3315 2.0919 -1.593 0.111256 +EthN:AgeF3:LrnSL NA NA NA NA +EthN:AgeF1:SexM -0.3105 1.2055 -0.258 0.796735 +EthN:AgeF2:SexM 0.3469 1.7965 0.193 0.846875 +EthN:AgeF3:SexM 0.8329 0.8970 0.929 0.353092 +EthN:LrnSL:SexM -0.1639 1.5250 -0.107 0.914411 +AgeF1:LrnSL:SexM -2.4285 1.4201 -1.710 0.087246 . +AgeF2:LrnSL:SexM -4.1914 1.6201 -2.587 0.009679 ** +AgeF3:LrnSL:SexM NA NA NA NA +EthN:AgeF1:LrnSL:SexM 2.1711 1.9192 1.131 0.257963 +EthN:AgeF2:LrnSL:SexM 2.1029 2.3444 0.897 0.369718 +EthN:AgeF3:LrnSL:SexM NA NA NA NA +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for Negative Binomial(1.9284) family taken to be 1) + + Null deviance: 272.29 on 145 degrees of freedom +Residual deviance: 167.45 on 118 degrees of freedom +AIC: 1097.3 + +Number of Fisher Scoring iterations: 1 + + + Theta: 1.928 + Std. Err.: 0.269 + + 2 x log-likelihood: -1039.324 +> ## IGNORE_RDIFF_END +> +> +> +> cleanEx() +> nameEx("summary.rlm") +> ### * summary.rlm +> +> flush(stderr()); flush(stdout()) +> +> ### Name: summary.rlm +> ### Title: Summary Method for Robust Linear Models +> ### Aliases: summary.rlm print.summary.rlm +> ### Keywords: robust +> +> ### ** Examples +> +> summary(rlm(calls ~ year, data = phones, maxit = 50)) + +Call: rlm(formula = calls ~ year, data = phones, maxit = 50) +Residuals: + Min 1Q Median 3Q Max +-18.314 -5.953 -1.681 26.460 173.769 + +Coefficients: + Value Std. Error t value +(Intercept) -102.6222 26.6082 -3.8568 +year 2.0414 0.4299 4.7480 + +Residual standard error: 9.032 on 22 degrees of freedom +> +> +> +> cleanEx() +> nameEx("theta.md") +> ### * theta.md +> +> flush(stderr()); flush(stdout()) +> +> ### Name: theta.md +> ### Title: Estimate theta of the Negative Binomial +> ### Aliases: theta.md theta.ml theta.mm +> ### Keywords: models +> +> ### ** Examples +> +> quine.nb <- glm.nb(Days ~ .^2, data = quine) +> theta.md(quine$Days, fitted(quine.nb), dfr = df.residual(quine.nb)) +[1] 1.135441 +> theta.ml(quine$Days, fitted(quine.nb)) +[1] 1.603641 +attr(,"SE") +[1] 0.2138379 +> theta.mm(quine$Days, fitted(quine.nb), dfr = df.residual(quine.nb)) +[1] 1.562879 +> +> ## weighted example +> yeast <- data.frame(cbind(numbers = 0:5, fr = c(213, 128, 37, 18, 3, 1))) +> fit <- glm.nb(numbers ~ 1, weights = fr, data = yeast) +> ## IGNORE_RDIFF_BEGIN +> summary(fit) + +Call: +glm.nb(formula = numbers ~ 1, data = yeast, weights = fr, init.theta = 3.586087428, + link = log) + +Coefficients: + Estimate Std. Error z value Pr(>|z|) +(Intercept) -0.38199 0.06603 -5.785 7.25e-09 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +(Dispersion parameter for Negative Binomial(3.5861) family taken to be 1) + + Null deviance: 408.9 on 5 degrees of freedom +Residual deviance: 408.9 on 5 degrees of freedom +AIC: 897.06 + +Number of Fisher Scoring iterations: 1 + + + Theta: 3.59 + Std. Err.: 1.75 + + 2 x log-likelihood: -893.063 +> ## IGNORE_RDIFF_END +> mu <- fitted(fit) +> theta.md(yeast$numbers, mu, dfr = 399, weights = yeast$fr) +[1] 3.027079 +> theta.ml(yeast$numbers, mu, limit = 15, weights = yeast$fr) +[1] 3.586087 +attr(,"SE") +[1] 1.749609 +> theta.mm(yeast$numbers, mu, dfr = 399, weights = yeast$fr) +[1] 3.549593 +> +> +> +> cleanEx() +> nameEx("ucv") +> ### * ucv +> +> flush(stderr()); flush(stdout()) +> +> ### Name: ucv +> ### Title: Unbiased Cross-Validation for Bandwidth Selection +> ### Aliases: ucv +> ### Keywords: dplot +> +> ### ** Examples +> +> ucv(geyser$duration) +Warning in ucv(geyser$duration) : + minimum occurred at one end of the range +[1] 0.1746726 +> +> +> +> cleanEx() +> nameEx("waders") +> ### * waders +> +> flush(stderr()); flush(stdout()) +> +> ### Name: waders +> ### Title: Counts of Waders at 15 Sites in South Africa +> ### Aliases: waders +> ### Keywords: datasets +> +> ### ** Examples +> +> plot(corresp(waders, nf=2)) +> +> +> +> cleanEx() +> nameEx("whiteside") +> ### * whiteside +> +> flush(stderr()); flush(stdout()) +> +> ### Name: whiteside +> ### Title: House Insulation: Whiteside's Data +> ### Aliases: whiteside +> ### Keywords: datasets +> +> ### ** Examples +> +> require(lattice) +Loading required package: lattice +> xyplot(Gas ~ Temp | Insul, whiteside, panel = ++ function(x, y, ...) { ++ panel.xyplot(x, y, ...) ++ panel.lmline(x, y, ...) ++ }, xlab = "Average external temperature (deg. C)", ++ ylab = "Gas consumption (1000 cubic feet)", aspect = "xy", ++ strip = function(...) strip.default(..., style = 1)) +> +> gasB <- lm(Gas ~ Temp, whiteside, subset = Insul=="Before") +> gasA <- update(gasB, subset = Insul=="After") +> summary(gasB) + +Call: +lm(formula = Gas ~ Temp, data = whiteside, subset = Insul == + "Before") + +Residuals: + Min 1Q Median 3Q Max +-0.62020 -0.19947 0.06068 0.16770 0.59778 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +(Intercept) 6.85383 0.11842 57.88 <2e-16 *** +Temp -0.39324 0.01959 -20.08 <2e-16 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 0.2813 on 24 degrees of freedom +Multiple R-squared: 0.9438, Adjusted R-squared: 0.9415 +F-statistic: 403.1 on 1 and 24 DF, p-value: < 2.2e-16 + +> summary(gasA) + +Call: +lm(formula = Gas ~ Temp, data = whiteside, subset = Insul == + "After") + +Residuals: + Min 1Q Median 3Q Max +-0.97802 -0.11082 0.02672 0.25294 0.63803 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +(Intercept) 4.72385 0.12974 36.41 < 2e-16 *** +Temp -0.27793 0.02518 -11.04 1.05e-11 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 0.3548 on 28 degrees of freedom +Multiple R-squared: 0.8131, Adjusted R-squared: 0.8064 +F-statistic: 121.8 on 1 and 28 DF, p-value: 1.046e-11 + +> gasBA <- lm(Gas ~ Insul/Temp - 1, whiteside) +> summary(gasBA) + +Call: +lm(formula = Gas ~ Insul/Temp - 1, data = whiteside) + +Residuals: + Min 1Q Median 3Q Max +-0.97802 -0.18011 0.03757 0.20930 0.63803 + +Coefficients: + Estimate Std. Error t value Pr(>|t|) +InsulBefore 6.85383 0.13596 50.41 <2e-16 *** +InsulAfter 4.72385 0.11810 40.00 <2e-16 *** +InsulBefore:Temp -0.39324 0.02249 -17.49 <2e-16 *** +InsulAfter:Temp -0.27793 0.02292 -12.12 <2e-16 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 + +Residual standard error: 0.323 on 52 degrees of freedom +Multiple R-squared: 0.9946, Adjusted R-squared: 0.9942 +F-statistic: 2391 on 4 and 52 DF, p-value: < 2.2e-16 + +> +> gasQ <- lm(Gas ~ Insul/(Temp + I(Temp^2)) - 1, whiteside) +> coef(summary(gasQ)) + Estimate Std. Error t value Pr(>|t|) +InsulBefore 6.759215179 0.150786777 44.826312 4.854615e-42 +InsulAfter 4.496373920 0.160667904 27.985514 3.302572e-32 +InsulBefore:Temp -0.317658735 0.062965170 -5.044991 6.362323e-06 +InsulAfter:Temp -0.137901603 0.073058019 -1.887563 6.489554e-02 +InsulBefore:I(Temp^2) -0.008472572 0.006624737 -1.278930 2.068259e-01 +InsulAfter:I(Temp^2) -0.014979455 0.007447107 -2.011446 4.968398e-02 +> +> gasPR <- lm(Gas ~ Insul + Temp, whiteside) +> anova(gasPR, gasBA) +Analysis of Variance Table + +Model 1: Gas ~ Insul + Temp +Model 2: Gas ~ Insul/Temp - 1 + Res.Df RSS Df Sum of Sq F Pr(>F) +1 53 6.7704 +2 52 5.4252 1 1.3451 12.893 0.0007307 *** +--- +Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +> options(contrasts = c("contr.treatment", "contr.poly")) +> gasBA1 <- lm(Gas ~ Insul*Temp, whiteside) +> coef(summary(gasBA1)) + Estimate Std. Error t value Pr(>|t|) +(Intercept) 6.8538277 0.13596397 50.409146 7.997414e-46 +InsulAfter -2.1299780 0.18009172 -11.827185 2.315921e-16 +Temp -0.3932388 0.02248703 -17.487358 1.976009e-23 +InsulAfter:Temp 0.1153039 0.03211212 3.590665 7.306852e-04 +> +> +> +> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly")) +> cleanEx() + +detaching ‘package:lattice’ + +> nameEx("width.SJ") +> ### * width.SJ +> +> flush(stderr()); flush(stdout()) +> +> ### Name: width.SJ +> ### Title: Bandwidth Selection by Pilot Estimation of Derivatives +> ### Aliases: width.SJ +> ### Keywords: dplot +> +> ### ** Examples +> +> width.SJ(geyser$duration, method = "dpi") +[1] 0.5747852 +> width.SJ(geyser$duration) +[1] 0.360518 +> +> width.SJ(galaxies, method = "dpi") +[1] 3256.151 +> width.SJ(galaxies) +[1] 2566.423 +> +> +> +> cleanEx() +> nameEx("wtloss") +> ### * wtloss +> +> flush(stderr()); flush(stdout()) +> +> ### Name: wtloss +> ### Title: Weight Loss Data from an Obese Patient +> ### Aliases: wtloss +> ### Keywords: datasets +> +> ### ** Examples +> +> ## IGNORE_RDIFF_BEGIN +> wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th), ++ data = wtloss, start = list(b0=90, b1=95, th=120)) +> wtloss.fm +Nonlinear regression model + model: Weight ~ b0 + b1 * 2^(-Days/th) + data: wtloss + b0 b1 th + 81.37 102.68 141.91 + residual sum-of-squares: 39.24 + +Number of iterations to convergence: 3 +Achieved convergence tolerance: 4.389e-06 +> ## IGNORE_RDIFF_END +> plot(wtloss) +> with(wtloss, lines(Days, fitted(wtloss.fm))) +> +> +> +> ### *