Nothing
Nest <-
function(y, nested_model, full_model, method=c("lm", "logit"),
data=d, digits_d=NULL, ...) {
# a dot in a parameter name to an underscore
dots <- list(...)
if (!is.null(dots)) if (length(dots) > 0) {
change <- c("nested.model", "full.model", "digits.d")
for (i in 1:length(dots)) {
if (names(dots)[i] %in% change) {
nm <- gsub(".", "_", names(dots)[i], fixed=TRUE)
assign(nm, dots[[i]])
get(nm)
}
}
}
method <- match.arg(method)
my.vars <- as.list(seq_along(data))
names(my.vars) <- names(data)
# get response variable
rv <- eval(substitute(y), envir=my.vars, enclos=parent.frame())
y.name <- names(my.vars)[rv]
if (is.null(digits_d)) digits_d <- .getdigits(data[,y.name], 2)
options(digits_d=digits_d)
if (method == "logit") {
is.bin <- TRUE
if (is.character(data[,y.name]))
data[,y.name] <- factor(data[,y.name])
if (is.factor(data[,y.name])) {
if (nlevels(data[,y.name]) != 2) is.bin <- FALSE
}
else {
for (i in 1:nrow(data))
if (!is.na(data[i,y.name]))
if (data[i,y.name]!=0 && data[i,y.name]!=1) is.bin <- FALSE
}
if (!is.bin) {
cat("\n"); stop(call.=FALSE, "\n","------\n",
"Response variable: ", y.name, "\n",
"If numeric, can only have values of 0 or 1.\n",
"If a factor, can only have two levels.\n\n")
}
}
# get nested model
nest.vars <- eval(substitute(nested_model), envir=my.vars, enclos=parent.frame())
n.pred <- length(nest.vars)
x.name <- character(length=n.pred)
for (ivar in 1:n.pred) {
x.name[ivar] <- names(my.vars)[nest.vars[ivar]]
}
xn.preds <- x.name[1]
if (n.pred > 1)
for(ivar in 2:n.pred) xn.preds <- paste(xn.preds, "+", x.name[ivar])
n.formula <- as.formula(paste(y.name, "~", xn.preds))
# get full model and analyze
full.vars <- eval(substitute(full_model), envir=my.vars, enclos=parent.frame())
full.vars <- union(nest.vars, full.vars)
n.pred <- length(full.vars)
x.name <- character(length=n.pred)
for (ivar in 1:n.pred) x.name[ivar] <- names(my.vars)[full.vars[ivar]]
xf.preds <- x.name[1]
if (n.pred > 1)
for(ivar in 2:n.pred) xf.preds <- paste(xf.preds, "+", x.name[ivar])
f.formula <- as.formula(paste(y.name, "~", xf.preds))
# do the model comparison
if (method=="lm") {
lm.full <- lm(f.formula, data=data)
lm.nest <- lm(n.formula, data=lm.full$model)
av <- anova(lm.nest, lm.full)
tx <- character(length = 0)
models <- attr(av, "heading")[2]
models <- sub("Model 1", "Reduced Model", models)
models <- sub("Model 2", "Full Model ", models)
tx[length(tx)+1] <- models
txtmdl <- tx
tx <- character(length = 0)
if (is.null(options()$knitr.in.progress)) {
tx[length(tx)+1] <- "Analysis of Variance"
tx[length(tx)+1] <- ""
}
max.c1 <- nchar("Residual")
mr.df <- av$Df[2]; mr.ss <- av$`Sum of Sq`[2]; mr.ms <- mr.ss/mr.df
mr.f <- av$F[2]; mr.p <- av$`Pr(>F)`[2]
mf.df <- av$Res.Df[2]; mf.ss <- av$RSS[2]; mf.ms <- mf.ss/mf.df
# width of columns
d <- digits_d
max.ln <- integer(length=0)
max.ln[1] <- max(nchar(.fmti(mr.df)),nchar(.fmti(mf.df)),nchar("df"))
max.ln[2] <- max(nchar(.fmt(mr.ss,d)),nchar(.fmt(mf.ss,d)),nchar("Sum Sq"))
max.ln[3] <- max(nchar(.fmt(mr.f,d)),nchar("Mean Sq"))
max.ln[4] <- nchar(.fmt(mr.p,d))
max.ln <- max.ln + 2
df.lbl <- .fmtc(" df", max.ln[1]+1)
SS.lbl <- .fmtc(" Sum Sq", max.ln[2]+1)
MS.lbl <- .fmtc("Mean Sq", max.ln[3]+1)
fv.lbl <- .fmtc("F-value", max.ln[4]+3)
tx[length(tx)+1] <- paste(eval(format("", width=max.c1+max.ln[1]-6)),
df.lbl, SS.lbl, MS.lbl, fv.lbl, " p-value", sep="")
if (n.pred > 0) {
rlb <- .fmtc("Tested", max.c1, j="left")
md <- .fmti(mr.df, max.ln[1])
ms <- .fmt(mr.ss, digits_d, max.ln[2])
mm <- .fmt(mr.ms, digits_d, max.ln[3])
mf <- .fmt(mr.f, digits_d, max.ln[4]+2)
mp <- .fmt(mr.p, 3, 9)
tx[length(tx)+1] <- paste(rlb, md, ms, mm, mf, mp)
tst <- c(mr.df, mr.ss, mr.ms, mr.f, mr.p)
names(tst) <- c("df", "ss", "ms", "fvalue", "pvalue")
# Residuals
rlb <- .fmtc("Residual", max.c1, j="left")
df <- .fmti(mf.df, max.ln[1])
SS <- .fmt(mf.ss, digits_d, max.ln[2])
MS <- .fmt(mf.ms, digits_d, max.ln[3])
tx[length(tx)+1] <- paste(rlb, df, SS, MS)
rsd <- c(mf.df, mf.ss, mf.ms)
names(rsd) <- c("df", "ss", "ms")
tot <- NA
txtbl <- tx
}
} # end method="lm"
else {
lm.full <- suppressWarnings(glm(f.formula, data=data, family="binomial"))
lm.nest <- suppressWarnings(glm(n.formula, data=lm.full$model, family="binomial"))
av <- anova(lm.nest, lm.full, test="Chisq")
tx <- character(length = 0)
models <- attr(av, "heading")[2]
models <- sub("Model 1", "Reduced Model", models)
models <- sub("Model 2", "Full Model ", models)
tx[length(tx)+1] <- models
txtmdl <- tx
tx <- character(length = 0)
if (is.null(options()$knitr.in.progress)) {
tx[length(tx)+1] <- "Analysis of Deviance"
tx[length(tx)+1] <- ""
}
max.c1 <- nchar("Residual")
mr.df <- av$Df[2]; mr.dv <- av$`Deviance`[2]; mr.p <- av$`Pr(>Chi)`[2]
mf.df <- av$`Resid. Df`[2]; mf.dv <- av$`Resid. Dev`[2];
mt.df <- av$`Resid. Df`[1]; mt.dv <- av$`Resid. Dev`[1];
# width of columns
d <- digits_d
max.ln <- integer(length=0)
max.ln[1] <- max(nchar(.fmti(mr.df)),nchar(.fmti(mf.df)),nchar("df"))
max.ln[2] <- nchar(.fmt(mr.dv,d)) + 2
max.ln <- max.ln + 2
df.lbl <- .fmtc(" df", max.ln[1]+1)
dv.lbl <- .fmtc(" Deviance", max.ln[2]+1)
tx[length(tx)+1] <- paste(eval(format("", width=max.c1-max.ln[1]+4)),
df.lbl, dv.lbl, " p-value", sep="")
if (n.pred > 0) {
rlb <- .fmtc("Tested", max.c1, j="left")
md <- .fmti(mr.df, max.ln[1])
dv <- .fmt(mr.dv, digits_d, max.ln[2])
mp <- .fmt(mr.p, 3, 9)
tx[length(tx)+1] <- paste(rlb, md, dv, mp)
tst <- c(mr.df, mr.dv, mr.p)
names(tst) <- c("df", "dev", "pvalue")
# Residuals
rlb <- .fmtc("Residual", max.c1, j="left")
df <- .fmti(mf.df, max.ln[1])
dev <- .fmt(mf.dv, digits_d, max.ln[2])
tx[length(tx)+1] <- paste(rlb, df, dev)
rsd <- c(mf.df, mf.dv)
names(rsd) <- c("df", "dev")
# Total
rlb <- .fmtc("Total", max.c1, j="left")
df <- .fmti(mt.df, max.ln[1])
dev <- .fmt(mt.dv, digits_d, max.ln[2])
tx[length(tx)+1] <- ""
tx[length(tx)+1] <- paste(rlb, df, dev)
tot <- c(mf.df, mf.dv)
names(tot) <- c("df", "dev")
txtbl <- tx
}
} # end method="logit"
class(txtmdl) <- "out"
class(txtbl) <- "out"
output <- list(
fun_call=match.call(), out_models=txtmdl, out_anova=txtbl,
anova_tested=tst, anova_residual=rsd, anova_total=tot
)
class(output) <- "out_all"
return(output)
}
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.