library(MASS);set.seed(1) ### 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) # same via loglm loglm(Claims ~ District + Group + Age + offset(log(Holders)), data = Insurance) ### 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 == 0) 1:ncol(M) else - (1:tmp$rank) qr.Q(tmp, complete = TRUE)[, set, drop = FALSE] } ### 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)) nls(Correct/Trials ~ fp1(Loud, L75, scal), data = OME[OME$Noise == "coherent",], start=c(L75=45, scal=3)) nls(Correct/Trials ~ fp1(Loud, L75, scal), data = OME[OME$Noise == "incoherent",], start = c(L75=45, scal=3)) # 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)) summary(lm(L75 ~ Noise/(Age + OME), data = OMEif, subset = (Age >= 30 & Age <= 60), na.action = na.omit), cor = FALSE) # 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)) nls(0 ~ fpl75(Correct, Trials, Loud, L75, scal), data = OME[OME$Noise == "incoherent",], start = c(L75=45, scal=3)) # 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))) 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))) # Now allow random effects by using NLME OMEf <- OME[rep(1:nrow(OME), OME$Trials),] attach(OME) OMEf$Resp <- rep(rep(c(1,0), length(Trials)), t(cbind(Correct, Trials-Correct))) OMEf <- OMEf[, -match(c("Correct", "Trials"), names(OMEf))] detach("OME") ## Not run: ## this fails in R on some 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 G1.nlme <- 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 summary(G1.nlme) ##D ##D G2.nlme <- 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) ##D summary(G2.nlme) ## End(Not run) ### 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)) ### 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") house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family=poisson, data=housing) addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test="Chisq") house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test = "Chisq") ### 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) anova(m2) ### Name: area ### Title: Adaptive Numerical Integration ### Aliases: area ### Keywords: nonlinear ### ** Examples area(sin, 0, pi) # integrate the sin function from 0 to pi. ### 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 summary(glm(y ~ trt * week, binomial, data = bacteria)) summary(glm(y ~ trt + week, binomial, data = bacteria)) summary(glm(y ~ trt + I(week > 2), binomial, data = bacteria)) # conditional random-effects analysis library(survival) bacteria$Time <- rep(1, nrow(bacteria)) coxph(Surv(Time, unclass(y)) ~ week + strata(ID), data = bacteria, method = "exact") coxph(Surv(Time, unclass(y)) ~ factor(week) + strata(ID), data = bacteria, method = "exact") coxph(Surv(Time, unclass(y)) ~ I(week > 2) + strata(ID), data = bacteria, method = "exact") # PQL glmm analysis library(nlme) summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID, family = binomial, data = bacteria)) ### 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) } ### Name: bcv ### Title: Biased Cross-Validation for Bandwidth Selection ### Aliases: bcv ### Keywords: dplot ### ** Examples bcv(geyser$duration) ### Name: beav1 ### Title: Body Temperature Series of Beaver 1 ### Aliases: beav1 ### Keywords: datasets ### ** Examples attach(beav1) 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]) 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, cor = FALSE) detach("beav1"); rm(temp, activ) ### 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]) arima(temp, order = c(1,0,0), xreg = activ) 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)) library(nlme) beav2.gls <- gls(temp ~ activ, data = beav2, corr = corAR1(0.8), method = "ML") summary(beav2.gls) summary(update(beav2.gls, subset = 6:100)) detach("beav2"); rm(temp, activ) ### Name: birthwt ### Title: Risk Factors Associated with Low Infant Birth Weight ### Aliases: birthwt ### Keywords: datasets ### ** Examples attach(birthwt) race <- factor(race, labels = c("white", "black", "other")) ptd <- factor(ptl > 0) ftv <- factor(ftv) levels(ftv)[-(1:2)] <- "2+" bwt <- data.frame(low = factor(low), age, lwt, race, smoke = (smoke > 0), ptd, ht = (ht > 0), ui = (ui > 0), ftv) detach("birthwt") options(contrasts = c("contr.treatment", "contr.poly")) glm(low ~ ., binomial, bwt) ### Name: boxcox ### Title: Box-Cox Transformations for Linear Models ### Aliases: boxcox boxcox.default boxcox.formula boxcox.lm ### Keywords: regression models hplot ### ** Examples data(trees) boxcox(Volume ~ log(Height) + log(Girth), data = trees, lambda = seq(-0.25, 0.25, length = 10)) boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = seq(-0.05, 0.45, len = 20)) ### Name: caith ### Title: Colours of Eyes and Hair of People in Caithness ### Aliases: caith ### Keywords: datasets ### ** Examples corresp(caith) 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)) ### Name: cement ### Title: Heat Evolved by Setting Cements ### Aliases: cement ### Keywords: datasets ### ** Examples lm(y ~ x1 + x2 + x3 + x4, cement) ### Name: confint-MASS ### Title: Confidence Intervals for Model Parameters ### Aliases: confint.glm confint.nls confint.profile.glm ### confint.profile.nls profile.glm ### Keywords: models ### ** Examples expn1 <- deriv(y ~ b0 + b1 * 2^(-x/th), c("b0", "b1", "th"), function(b0, b1, th, x) {}) wtloss.gr <- nls(Weight ~ expn1(b0, b1, th, Days), data = wtloss, start = c(b0=90, b1=95, th=120)) expn2 <- deriv(~b0 + b1*((w0 - b0)/b1)^(x/d0), c("b0","b1","d0"), function(b0, b1, d0, x, w0) {}) wtloss.init <- function(obj, w0) { p <- coef(obj) d0 <- - log((w0 - p["b0"])/p["b1"])/log(2) * p["th"] c(p[c("b0", "b1")], d0 = as.vector(d0)) } out <- NULL w0s <- c(110, 100, 90) for(w0 in w0s) { fm <- nls(Weight ~ expn2(b0, b1, d0, Days, w0), wtloss, start = wtloss.init(wtloss.gr, w0)) out <- rbind(out, c(coef(fm)["d0"], confint(fm, "d0"))) } dimnames(out) <- list(paste(w0s, "kg:"), c("d0", "low", "high")) out ldose <- rep(0:5, 2) numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16) sex <- factor(rep(c("M", "F"), c(6, 6))) SF <- cbind(numdead, numalive = 20 - numdead) budworm.lg0 <- glm(SF ~ sex + ldose - 1, family = binomial) confint(budworm.lg0) confint(budworm.lg0, "ldose") ### Name: contr.sdif ### Title: Successive Differences contrast coding ### Aliases: contr.sdif ### Keywords: models ### ** Examples contr.sdif(6) ### 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 (ct <- corresp(~ Age + Eth, data = quine)) ## Not run: plot(ct) corresp(caith) biplot(corresp(caith, nf = 2)) ### Name: cov.rob ### Title: Resistant Estimation of Multivariate Location and Scatter ### Aliases: cov.rob cov.mve cov.mcd ### Keywords: robust multivariate ### ** Examples data(stackloss) set.seed(123) cov.rob(stackloss) cov.rob(stack.x, method = "mcd", nsamp = "exact") ### Name: cov.trob ### Title: Covariance Estimation for Multivariate t Distribution ### Aliases: cov.trob ### Keywords: multivariate ### ** Examples data(stackloss) cov.trob(stackloss) ### 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) ## Not run: ~ (.v1 + .v2 + .v3)^3 + a/b ### 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.p(update(budworm.lg0, family = binomial(link=probit)), cf = c(1,3), p = 1:3/4) ### 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") quine.stp <- stepAIC(quine.nxt, scope = list(upper = ~Eth*Sex*Age*Lrn, lower = ~1), trace = FALSE) dropterm(quine.stp, test = "F") quine.3 <- update(quine.stp, . ~ . - Eth:Age:Lrn) dropterm(quine.3, test = "F") quine.4 <- update(quine.3, . ~ . - Eth:Age) dropterm(quine.4, test = "F") quine.5 <- update(quine.4, . ~ . - Age:Lrn) dropterm(quine.5, test = "F") 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") ### 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) prof <- profile(eagles.glm) plot(prof) pairs(prof) ### Name: epil ### Title: Seizure Counts for Epileptics ### Aliases: epil ### Keywords: datasets ### ** Examples summary(glm(y ~ lbase*trt + lage + V4, family = poisson, data = epil), cor = FALSE) 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"))) summary(glm(y ~ pred + factor(subject) + offset(log(time)), family = poisson, data = epil3), cor = FALSE) summary(glmmPQL(y ~ lbase*trt + lage + V4, random = ~ 1 | subject, family = poisson, data = epil)) summary(glmmPQL(y ~ pred, random = ~1 | subject, family = poisson, data = epil3)) ### 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) ### Name: fitdistr ### Title: Maximum-likelihood Fitting of Univariate Distributions ### Aliases: fitdistr print.fitdistr coef.fitdistr logLik.fitdistr ### Keywords: distribution htest ### ** Examples set.seed(123) x <- rgamma(100, shape = 5, rate = 0.1) fitdistr(x, "gamma") ## now do this directly with more control. fitdistr(x, dgamma, list(shape = 1, rate = 0.1), lower = 0.01) set.seed(123) x2 <- rt(250, df = 9) fitdistr(x2, "t", df = 9) ## allow df to vary: not a very good idea! fitdistr(x2, "t") ## 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)) set.seed(123) x3 <- rweibull(100, shape = 4, scale = 100) fitdistr(x3, "weibull") set.seed(123) x4 <- rnegbin(500, mu = 5, theta = 4) fitdistr(x4, "Negative Binomial") # R only ### 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) solve(X, X/5) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 2.0000e-01 3.7199e-17 1.2214e-16 5.7887e-17 -8.7841e-17 ## [2,] -1.1473e-16 2.0000e-01 7.0955e-17 2.0300e-17 -1.0566e-16 ## [3,] 2.7975e-16 1.3653e-17 2.0000e-01 -1.3397e-16 1.5577e-16 ## [4,] -2.9196e-16 2.0412e-17 1.5618e-16 2.0000e-01 -2.1921e-16 ## [5,] -3.6476e-17 -3.6430e-17 3.6432e-17 4.7690e-17 2.0000e-01 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 ### Name: galaxies ### Title: Velocities for 82 Galaxies ### Aliases: galaxies ### Keywords: datasets ### ** Examples gal <- galaxies/1000 c(width.SJ(gal, method = "dpi"), width.SJ(gal)) 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) ### 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) ## Not run: ##D Alpha: 538.13 ##D SE: 253.60 ## End(Not run) 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) ## Not run: ##D Initial estimate: 1.0603 ##D Iter. 1 Alpha: 1.23840774338543 ##D Iter. 2 Alpha: 1.27699745778205 ##D Iter. 3 Alpha: 1.27834332265501 ##D Iter. 4 Alpha: 1.27834485787226 ##D ##D Alpha: 1.27834 ##D SE: 0.13452 ## End(Not run) summary(gm, dispersion = gamma.dispersion(gm)) # better summary ### 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) survreg(Surv(time, cens) ~ factor(pair) + treat, gehan, dist = "exponential") summary(survreg(Surv(time, cens) ~ treat, gehan, dist = "exponential")) summary(survreg(Surv(time, cens) ~ treat, gehan)) gehan.cox <- coxph(Surv(time, cens) ~ treat, gehan) summary(gehan.cox) ### Name: ginv ### Title: Generalized Inverse of a Matrix ### Aliases: ginv ### Keywords: algebra ### ** Examples ## Not run: ##D # The function is currently defined as ##D function(X, tol = sqrt(.Machine$double.eps)) ##D { ##D ## Generalized Inverse of a Matrix ##D dnx <- dimnames(X) ##D if(is.null(dnx)) dnx <- vector("list", 2) ##D s <- svd(X) ##D nz <- s$d > tol * s$d[1] ##D structure( ##D if(any(nz)) s$v[, nz] %*% (t(s$u[, nz])/s$d[nz]) else X, ##D dimnames = dnx[2:1]) ##D } ## End(Not run) ### 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) ### 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) ## 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))) try(glm.nb(y ~ lag1+lag2+lag3, link=identity)) glm.nb(y ~ lag1+lag2+lag3, link=identity, start=c(2, 0.1, 0.1, 0.1)) glm.nb(y ~ lag1+lag2+lag3, link=identity, start=coef(fit)) glm.nb(y ~ lag1+lag2+lag3, link=identity, etastart=rep(5, 42)) ## End Don't show ### Name: glmmPQL ### Title: Fit Generalized Linear Mixed Models via PQL ### Aliases: glmmPQL ### Keywords: models ### ** Examples library(nlme) # will be loaded automatically if omitted summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID, family = binomial, data = bacteria)) ## Don't show: # an example of offset summary(glmmPQL(y ~ trt + week, random = ~ 1 | ID, family = binomial, data = bacteria)) summary(glmmPQL(y ~ trt + week + offset(week), random = ~ 1 | ID, family = binomial, data = bacteria)) ## End Don't show ### 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) summary(house.glm0, cor = FALSE) addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test = "Chisq") house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont)) summary(house.glm1, cor = FALSE) 1 - pchisq(deviance(house.glm1), house.glm1$df.residual) dropterm(house.glm1, test = "Chisq") addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test = "Chisq") 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)) # Iterative proportional scaling loglm(Freq ~ Infl*Type*Cont + Sat*(Infl+Type+Cont), data = housing) # multinomial model library(nnet) (house.mult<- multinom(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)) house.mult2 <- multinom(Sat ~ Infl*Type*Cont, weights = Freq, data = housing) anova(house.mult, house.mult2) house.pm <- predict(house.mult, expand.grid(hnames[-1]), type = "probs") cbind(expand.grid(hnames[-1]), round(house.pm, 2)) # 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))) mean(ratio) (house.plr <- polr(Sat ~ Infl + Type + Cont, data = housing, weights = Freq)) house.pr1 <- predict(house.plr, expand.grid(hnames[-1]), type = "probs") cbind(expand.grid(hnames[-1]), round(house.pr1, 2)) Fr <- matrix(housing$Freq, ncol = 3, byrow = TRUE) 2*sum(Fr*log(house.pr/house.pr1)) house.plr2 <- stepAIC(house.plr, ~.^2) house.plr2$anova ### Name: huber ### Title: Huber M-estimator of Location with MAD Scale ### Aliases: huber ### Keywords: robust ### ** Examples huber(chem) ### Name: hubers ### Title: Huber Proposal 2 Robust Estimator of Location and/or Scale ### Aliases: hubers ### Keywords: robust ### ** Examples hubers(chem) hubers(chem, mu=3.68) ### 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) immer.aov <- aov((Y1+Y2)/2 ~ Var + Loc, data = immer) summary(immer.aov) model.tables(immer.aov, type = "means", se = TRUE, cterms = "Var") ### Name: isoMDS ### Title: Kruskal's Non-metric Multidimensional Scaling ### Aliases: isoMDS Shepard ### Keywords: multivariate ### ** Examples data(swiss) swiss.x <- as.matrix(swiss[, -1]) swiss.dist <- dist(swiss.x) swiss.mds <- isoMDS(swiss.dist) 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") ### 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") ### 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 data(iris3) 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]) ## 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 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.)) ### 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) ### Name: lm.ridge ### Title: Ridge Regression ### Aliases: lm.ridge plot.ridgelm print.ridgelm select select.ridgelm ### Keywords: models ### ** Examples data(longley) # not the same as the S-PLUS dataset names(longley)[1] <- "y" lm.ridge(y ~ ., longley) plot(lm.ridge(y ~ ., longley, lambda = seq(0,0.1,0.001))) select(lm.ridge(y ~ ., longley, lambda = seq(0,0.1,0.0001))) ### 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 minn38a <- array(0, c(3,4,7,2), lapply(minn38[, -5], levels)) minn38a[data.matrix(minn38[,-5])] <- minn38$fol fm <- loglm(~1 + 2 + 3 + 4, minn38a) # numerals as names. deviance(fm) ##[1] 3711.9 fm1 <- update(fm, .~.^2) fm2 <- update(fm, .~.^3, print = TRUE) ## 5 iterations: deviation 0.0750732 anova(fm, fm1, fm2) ## Not run: LR tests for hierarchical log-linear models ##D ##D Model 1: ##D ~ 1 + 2 + 3 + 4 ##D Model 2: ##D . ~ 1 + 2 + 3 + 4 + 1:2 + 1:3 + 1:4 + 2:3 + 2:4 + 3:4 ##D Model 3: ##D . ~ 1 + 2 + 3 + 4 + 1:2 + 1:3 + 1:4 + 2:3 + 2:4 + 3:4 + ##D 1:2:3 + 1:2:4 + 1:3:4 + 2:3:4 ##D ##D Deviance df Delta(Dev) Delta(df) P(> Delta(Dev) ##D Model 1 3711.915 155 ##D Model 2 220.043 108 3491.873 47 0.00000 ##D Model 3 47.745 36 172.298 72 0.00000 ##D Saturated 0.000 0 47.745 36 0.09114 ##D ## End(Not run) # Case 1. An array generated with xtabs. loglm(~ Type + Origin, xtabs(~ Type + Origin, Cars93)) ## Not run: Call: ##D loglm(formula = ~Type + Origin, data = xtabs(~Type + Origin, ##D Cars93)) ##D ##D Statistics: ##D X^2 df P(> X^2) ##D Likelihood Ratio 18.362 5 0.0025255 ##D Pearson 14.080 5 0.0151101 ##D ## End(Not run) # Case 2. Frequencies given as a vector in a data frame names(quine) ## [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.7 1368.7 c(fm$df, gm$df) # resid df do not! c(fm$df, gm$df.residual) # resid df do not! ## [1] 127 128 # The loglm residual degrees of freedom is wrong because of # a non-detectable redundancy in the model matrix. ### 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, len=20)) ### Name: lqs ### Title: Resistant Regression ### Aliases: lqs lqs.formula lqs.default lmsreg ltsreg ### Keywords: models robust ### ** Examples data(stackloss) set.seed(123) lqs(stack.loss ~ ., data = stackloss) lqs(stack.loss ~ ., data = stackloss, method = "S", nsamp = "exact") ### Name: mca ### Title: Multiple Correspondence Analysis ### Aliases: mca print.mca ### Keywords: category multivariate ### ** Examples farms.mca <- mca(farms, abbrev=TRUE) farms.mca plot(farms.mca) ### Name: menarche ### Title: Age of Menarche data ### Aliases: menarche ### Keywords: datasets ### ** Examples mprob <- glm(cbind(Menarche, Total - Menarche) ~ Age, binomial(link = probit), data = menarche) ### 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) # and predict at 130C unlist(predict(motor.wei, data.frame(temp=130), se.fit = TRUE)) motor.cox <- coxph(Surv(time, cens) ~ temp, motors) summary(motor.cox) # 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)) ) ### Name: muscle ### Title: Effect of Calcium Chloride on Muscle Contraction in Rat Hearts ### Aliases: muscle ### Keywords: datasets ### ** Examples 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) B 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)) attach(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) detach() require(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, ...) { lines(spline(x, y)) panel.xyplot(x, Muscle$logLength[subscripts], ...) }) ### 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 var(mvrnorm(n=1000, rep(0, 2), Sigma)) var(mvrnorm(n=1000, rep(0, 2), Sigma, empirical = TRUE)) ### 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) ### Name: nlschools ### Title: Eighth-Grade Pupils in the Netherlands ### Aliases: nlschools ### Keywords: datasets ### ** Examples library(nlme) nl1 <- nlschools attach(nl1) classMeans <- tapply(IQ, class, mean) nl1$IQave <- classMeans[as.character(class)] nl1$IQ <- nl1$IQ - nl1$IQave detach() cen <- c("IQ", "IQave", "SES") nl1[cen] <- scale(nl1[cen], center = TRUE, scale = FALSE) nl.lme <- lme(lang ~ IQ*COMB + IQave + SES, random = ~ IQ | class, data = nl1) summary(nl.lme) ### 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) npk.aov summary(npk.aov) alias(npk.aov) coef(npk.aov) options(contrasts = c("contr.treatment", "contr.poly")) npk.aov1 <- aov(yield ~ block + N + K, data = npk) summary.lm(npk.aov1) se.contrast(npk.aov1, list(N=="0", N=="1"), data = npk) model.tables(npk.aov1, type = "means", se = TRUE) ### 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) summary(oats.aov) summary(oats.aov, split = list(Nf=list(L=1, Dev=2:3))) 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) ### Name: parcoord ### Title: Parallel Coordinates Plot ### Aliases: parcoord ### Keywords: hplot ### ** Examples data(state) parcoord(state.x77[, c(7, 4, 6, 2, 5, 3)]) data(iris3) ir <- rbind(iris3[,,1], iris3[,,2], iris3[,,3]) parcoord(log(ir)[, c(3, 4, 2, 1)], col = 1 + (0:149)%/%50) ### 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) ### Name: plot.mca ### Title: Plot Method for Objects of Class 'mca' ### Aliases: plot.mca ### Keywords: hplot multivariate ### ** Examples plot(mca(farms, abbrev = TRUE)) ### 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 summary(house.plr) ## slightly worse fit from summary(update(house.plr, method = "probit")) ## although it is not really appropriate, can fit summary(update(house.plr, method = "cloglog")) predict(house.plr, housing, type = "p") addterm(house.plr, ~.^2, test = "Chisq") house.plr2 <- stepAIC(house.plr, ~.^2) house.plr2$anova anova(house.plr, house.plr2) house.plr <- update(house.plr, Hess=TRUE) pr <- profile(house.plr) confint(pr) plot(pr) pairs(pr) ### 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) predict(fit, bacteria, level = 0, type="response") predict(fit, bacteria, level = 1, type="response") ### Name: predict.lda ### Title: Classify Multivariate Observations by Linear Discrimination ### Aliases: predict.lda ### Keywords: multivariate ### ** Examples data(iris3) 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 ### Name: predict.lqs ### Title: Predict from an lqs Fit ### Aliases: predict.lqs ### Keywords: models ### ** Examples data(stackloss) set.seed(123) fm <- lqs(stack.loss ~ ., data = stackloss, method = "S", nsamp = "exact") predict(fm, stackloss) ### Name: predict.qda ### Title: Classify from Quadratic Discriminant Analysis ### Aliases: predict.qda ### Keywords: multivariate ### ** Examples data(iris3) 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 ### 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 data(iris3) 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 ### Name: rational ### Title: Rational Approximation ### Aliases: rational .rat ### Keywords: math ### ** Examples X <- matrix(runif(25), 5, 5) solve(X, X/5) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 2.0000e-01 3.7199e-17 1.2214e-16 5.7887e-17 -8.7841e-17 ## [2,] -1.1473e-16 2.0000e-01 7.0955e-17 2.0300e-17 -1.0566e-16 ## [3,] 2.7975e-16 1.3653e-17 2.0000e-01 -1.3397e-16 1.5577e-16 ## [4,] -2.9196e-16 2.0412e-17 1.5618e-16 2.0000e-01 -2.1921e-16 ## [5,] -3.6476e-17 -3.6430e-17 3.6432e-17 4.7690e-17 2.0000e-01 ## 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 ### 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 renumerate(.Last.value) ## ~ (1 + 2 + 3)^3 + a/b ### 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 data(stackloss) summary(rlm(stack.loss ~ ., stackloss)) rlm(stack.loss ~ ., stackloss, psi = psi.hampel, init = "lts") rlm(stack.loss ~ ., stackloss, psi = psi.bisquare) ### 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 data(Puromycin) 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))) rms.curv(Purfit1) ##Parameter effects: c^theta x sqrt(F) = 0.2121 ## Intrinsic: c^iota x sqrt(F) = 0.092 ### 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) ### Name: sammon ### Title: Sammon's Non-Linear Mapping ### Aliases: sammon ### Keywords: multivariate ### ** Examples data(swiss) swiss.x <- as.matrix(swiss[, -1]) swiss.sam <- sammon(dist(swiss.x)) plot(swiss.sam$points, type = "n") text(swiss.sam$points, labels = as.character(1:nrow(swiss.x))) ### 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 cpus1 <- cpus attach(cpus) for(v in names(cpus)[2:7]) cpus1[[v]] <- cut(cpus[[v]], unique(quantile(cpus[[v]])), include.lowest = TRUE) detach() 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 example(birthwt) birthwt.glm <- glm(low ~ ., family = binomial, data = bwt) birthwt.step <- stepAIC(birthwt.glm, trace = FALSE) birthwt.step$anova birthwt.step2 <- stepAIC(birthwt.glm, ~ .^2 + I(scale(age)^2) + I(scale(lwt)^2), trace = FALSE) birthwt.step2$anova quine.nb <- glm.nb(Days ~ .^4, data = quine) quine.nb2 <- stepAIC(quine.nb) quine.nb2$anova ### Name: summary.negbin ### Title: Summary Method Function for Objects of Class 'negbin' ### Aliases: summary.negbin print.summary.negbin ### Keywords: models ### ** Examples summary(glm.nb(Days ~ Eth*Age*Lrn*Sex, quine, link = log)) ### 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)) ## Not run: ##D Call: ##D rlm(formula = calls ~ year, data = phones, maxit = 50) ##D ##D Residuals: ##D Min 1Q Median 3Q Max ##D -18.31 -5.95 -1.68 26.46 173.77 ##D ##D Coefficients: ##D Value Std. Error t value ##D (Intercept) -102.622 26.553 -3.86 ##D year 2.041 0.429 4.76 ##D ##D Residual standard error: 9.03 on 22 degrees of freedom ##D ##D Correlation of Coefficients: ##D [1] -0.994 ##D ## End(Not run) ### 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)) theta.ml(quine$Days, fitted(quine.nb)) theta.mm(quine$Days, fitted(quine.nb), dfr = df.residual(quine.nb)) ## 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) summary(fit) attach(yeast) mu <- fitted(fit) theta.md(numbers, mu, dfr = 399, weights = fr) theta.ml(numbers, mu, weights = fr) theta.mm(numbers, mu, dfr = 399, weights = fr) detach() ### Name: ucv ### Title: Unbiased Cross-Validation for Bandwidth Selection ### Aliases: ucv ### Keywords: dplot ### ** Examples ucv(geyser$duration) ### Name: waders ### Title: Counts of Waders at 15 Sites in South Africa ### Aliases: waders ### Keywords: datasets ### ** Examples plot(corresp(waders, nf=2)) ### Name: whiteside ### Title: House Insulation: Whiteside's Data ### Aliases: whiteside ### Keywords: datasets ### ** Examples require(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) summary(gasA) gasBA <- lm(Gas ~ Insul/Temp - 1, whiteside) summary(gasBA) gasQ <- lm(Gas ~ Insul/(Temp + I(Temp^2)) - 1, whiteside) coef(summary(gasQ)) gasPR <- lm(Gas ~ Insul + Temp, whiteside) anova(gasPR, gasBA) options(contrasts = c("contr.treatment", "contr.poly")) gasBA1 <- lm(Gas ~ Insul*Temp, whiteside) coef(summary(gasBA1)) ### Name: width.SJ ### Title: Bandwidth Selection by Pilot Estimation of Derivatives ### Aliases: width.SJ ### Keywords: dplot ### ** Examples attach(geyser) width.SJ(duration, method = "dpi") width.SJ(duration) detach() width.SJ(galaxies, method = "dpi") width.SJ(galaxies) ### Name: wtloss ### Title: Weight Loss Data from an Obese Patient ### Aliases: wtloss ### Keywords: datasets ### ** Examples wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th), data = wtloss, start = list(b0=90, b1=95, th=120), trace = TRUE)