File: C:/Program Files/R/R-4.4.0/library/MASS/tests/Examples/MASS-Ex.Rout.save
R Under development (unstable) (2023-07-13 r84685) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin22.5.0
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)
<none> 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)
<none> 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)
<none> 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'
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("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)
<none> 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)
<none> 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)
<none> 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)
<none> 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)
<none> 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)
<none> 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
<none> 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
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")
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
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
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
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)
<none> 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)
<none> 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)
<none> 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
<none> 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
<none> 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
<none> 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
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
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
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("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)
<none> 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
<none> 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
<none> 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
<none> 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("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
<none> 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
<none> 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
<none> 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)))
>
>
>
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 2.902 0.152 3.684 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')