Nothing
# File MASS/profiles.q copyright (C) 1996 D. M. Bates and W. N. Venables.
#
# port to R by B. D. Ripley copyright (C) 1998
#
# corrections copyright (C) 2000,3,6,7 B. D. Ripley
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 or 3 of the License
# (at your option).
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
profile.glm <- function(fitted, which = 1:p, alpha = 0.01,
maxsteps = 10, del = zmax/5, trace = FALSE, ...)
{
Pnames <- names(B0 <- coef(fitted))
nonA <- !is.na(B0)
pv0 <- t(as.matrix(B0))
p <- length(Pnames)
if(is.character(which)) which <- match(which, Pnames)
summ <- summary(fitted)
std.err <- summ$coefficients[, "Std. Error", drop = FALSE] # unaliased only
mf <- model.frame(fitted)
Y <- model.response(mf) ### change
n <- NROW(Y) ### change: NROW() rather than length()
O <- model.offset(mf)
if(!length(O)) O <- rep(0, n)
W <- model.weights(mf)
if(length(W) == 0L) W <- rep(1, n)
OriginalDeviance <- deviance(fitted)
DispersionParameter <- summ$dispersion
X <- model.matrix(fitted)
fam <- family(fitted)
switch(fam$family,
binomial = ,
poisson = ,
`Negative Binomial` = {
zmax <- sqrt(qchisq(1 - alpha, 1))
profName <- "z"
}
,
gaussian = ,
quasi = ,
inverse.gaussian = ,
quasibinomial = ,
quasipoisson = ,
{
zmax <- sqrt(qf(1 - alpha, 1, n - p))
profName <- "tau"
}
)
prof <- vector("list", length=length(which))
names(prof) <- Pnames[which]
for(i in which) {
if(!nonA[i]) next
zi <- 0
pvi <- pv0
a <- nonA; a[i] <- FALSE
Xi <- X[, a, drop = FALSE]
pi <- Pnames[i]
for(sgn in c(-1, 1)) {
if(trace)
message("\nParameter: ", pi, " ",
c("down", "up")[(sgn + 1)/2 + 1])
step <- 0
z <- 0
## LP is the linear predictor including offset.
## we need to take care to avoid aliased-out coefficients.
LP <- X[, nonA, drop = FALSE] %*% B0[nonA] + O
while((step <- step + 1) < maxsteps && abs(z) < zmax) {
bi <- B0[i] + sgn * step * del * std.err[Pnames[i], 1]
o <- O + X[, i] * bi
## call to glm.fit.null not needed from 1.4.1 on
fm <- glm.fit(x = Xi, y = Y, weights = W, etastart = LP,
offset = o, family = fam,
control = fitted$control)
LP <- Xi %*% fm$coefficients + o
ri <- pv0
ri[, names(coef(fm))] <- coef(fm)
ri[, pi] <- bi
pvi <- rbind(pvi, ri)
zz <- (fm$deviance - OriginalDeviance)/DispersionParameter
if(zz > - 1e-3) zz <- max(zz, 0)
else stop("profiling has found a better solution, so original fit had not converged")
z <- sgn * sqrt(zz)
zi <- c(zi, z)
}
}
si <- order(zi)
prof[[pi]] <- structure(data.frame(zi[si]), names = profName)
prof[[pi]]$par.vals <- pvi[si, ,drop=FALSE]
}
val <- structure(prof, original.fit = fitted, summary = summ)
class(val) <- c("profile.glm", "profile")
val
}
plot.profile <-
## R version: non-Trellis-based replacement for plot.profile
function(x, ...)
{
nulls <- sapply(x, is.null)
if (all(nulls)) return(NULL)
x <- x[!nulls]
nm <- names(x)
nr <- ceiling(sqrt(length(nm)))
oldpar <- par(mfrow = c(nr, nr))
on.exit(par(oldpar))
for(nm in names(x)) {
tau <- x[[nm]][[1L]]
parval <- x[[nm]][[2L]][, nm]
dev.hold()
plot(parval, tau, xlab = nm, ylab = "tau", type="n")
## allow for profiling failures
if(sum(tau == 0) == 1) points(parval[tau == 0], 0, pch = 3)
splineVals <- spline(parval, tau)
lines(splineVals$x, splineVals$y)
dev.flush()
}
}
pairs.profile <-
## Another plot method for profile objects showing pairwise traces.
## Recommended only for diagnostic purposes.
function(x, colours = 2:3, ...)
{
parvals <- lapply(x, "[[", "par.vals")
rng <- apply(do.call("rbind", parvals), 2L, range, na.rm = TRUE)
Pnames <- colnames(rng)
npar <- length(Pnames)
coefs <- coef(attr(x, "original.fit"))
form <- paste(as.character(formula(attr(x, "original.fit")))[c(2, 1, 3)],
collapse = "")
oldpar <- par(mar = c(0, 0, 0, 0), mfrow = c(1, 1),
oma = c(3, 3, 6, 3), las = 1)
on.exit(par(oldpar))
##
## The following dodge ensures that the plot region is square
##
fin <- par("fin")
dif <- (fin[2L] - fin[1L])/2
if(dif > 0) adj <- c(dif, 0, dif, 0)
else adj <- c(0, - dif, 0, - dif)
par(omi = par("omi") + adj)
##
##
cex <- 1 + 1/npar
frame()
mtext(form, side = 3, line = 3, cex = 1.5, outer = TRUE)
del <- 1/npar
for(i in 1L:npar) {
ci <- npar - i
pi <- Pnames[i]
for(j in 1L:npar) {
dev.hold()
pj <- Pnames[j]
par(fig = del * c(j - 1, j, ci, ci + 1))
if(i == j) {
par(new=TRUE)
plot(rng[, pj], rng[, pi], axes = FALSE,
xlab = "", ylab = "", type = "n")
op <- par(usr = c(-1, 1, -1, 1))
text(0, 0, pi, cex = cex, adj = 0.5)
par(op)
} else {
col <- colours
if(i < j) col <- col[2:1]
if(!is.null(parvals[[pj]])) {
par(new=TRUE)
plot(spline(x <- parvals[[pj]][, pj],
y <- parvals[[pj]][, pi]),
type = "l", xlim = rng[, pj],
ylim = rng[, pi], axes = FALSE,
xlab = "", ylab = "", col = col[2L])
pu <- par("usr")
smidge <- 2/100 * (pu[4L] - pu[3L])
segments(x, pmax(pu[3L], y - smidge), x,
pmin(pu[4L], y + smidge))
} else
plot(rng[, pj], rng[, pi], axes = FALSE,
xlab = "", ylab = "", type = "n")
if(!is.null(parvals[[pi]])) {
lines(x <- parvals[[pi]][, pj], y <- parvals[[pi]][, pi],
type = "l", col = col[1L])
pu <- par("usr")
smidge <- 2/100 * (pu[2L] - pu[1L])
segments(pmax(pu[1L], x - smidge), y, pmin(pu[2L], x + smidge), y)
}
points(coefs[pj], coefs[pi], pch = 3, cex = 3)
}
if(i == npar) axis(1)
if(j == 1) axis(2)
if(i == 1) axis(3)
if(j == npar) axis(4)
dev.flush()
}
}
par(fig = c(0, 1, 0, 1))
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.