#' Misclassification in models using MCSIMEX
#'
#' Implementation of the misclassification MCSIMEX algorithm as described by Küchenhoff, Mwalili and Lesaffre.
#'
#' @aliases mcsimex print.mcsimex summary.mcsimex print.summary.mcsimex plot.mcsimex predict.mcsimex refit.mcsimex
#'
#' @param model the naive model, the misclassified variable must be a factor
#' @param SIMEXvariable vector of names of the variables for which the MCSIMEX-method should be applied
#' @param mc.matrix if one variable is misclassified it can be a matrix. If more
#' than one variable is misclassified it must be a list of the misclassification
#' matrices, names must match with the SIMEXvariable names, column- and row-names
#' must match with the factor levels. If a special misclassification is desired,
#' the name of a function can be specified (see details)
#' @param lambda vector of exponents for the misclassification matrix (without 0)
#' @param B number of iterations for each lambda
#' @param fitting.method \code{linear}, \code{quadratic} and \code{loglinear}
#' are implemented (first 4 letters are enough)
#' @param jackknife.estimation specifying the extrapolation method for jackknife
#' variance estimation. Can be set to \code{FALSE} if it should not be performed
#' @param asymptotic logical, indicating if asymptotic variance estimation should
#' be done, the option \code{x = TRUE} must be enabled in the naive model
#' @param x object of class 'mcsimex'
#' @param digits number of digits to be printed
#' @param object object of class 'mcsimex'
#' @param xlab optional name for the X-Axis
#' @param ylab vector containing the names for the Y-Axis
#' @param ask ogical. If \code{TRUE}, the user is asked for input, before a new figure is drawn
#' @param show vector of logicals indicating for which variables a plot should be produced
#' @param newdata optionally, a data frame in which to look for variables with
#' which to predict. If omitted, the fitted linear predictors are used
#' @param \dots arguments passed to other functions
#'
#' @return An object of class 'mcsimex' which contains:
#'
#' \item{coefficients}{corrected coefficients of the MCSIMEX model,}
#' \item{SIMEX.estimates}{the MCSIMEX-estimates of the coefficients for each lambda,}
#' \item{lambda}{the values of lambda,}
#' \item{model}{the naive model,}
#' \item{mc.matrix}{the misclassification matrix,}
#' \item{B}{the number of iterations,}
#' \item{extrapolation}{the model object of the extrapolation step,}
#' \item{fitting.method}{the fitting method used in the extrapolation step,}
#' \item{SIMEXvariable}{name of the SIMEXvariables,}
#' \item{call}{the function call,}
#' \item{variance.jackknife}{the jackknife variance estimates,}
#' \item{extrapolation.variance}{the model object of the variance extrapolation,}
#' \item{variance.jackknife.lambda}{the data set for the extrapolation,}
#' \item{variance.asymptotic}{the asymptotic variance estimates,}
#' \item{theta}{all estimated coefficients for each lambda and B,}
#' ...
#'
#' @details
#' If \code{mc.matrix} is a function the first argument of that function must
#' be the whole dataset used in the naive model, the second argument must be
#' the exponent (lambda) for the misclassification. The function must return
#' a \code{data.frame} containing the misclassified \code{SIMEXvariable}. An
#' example can be found below.
#'
#' Asymptotic variance estimation is only implemented for \code{lm} and \code{glm}
#'
#' The loglinear fit has the form \code{g(lambda, GAMMA) = exp(gamma0 + gamma1 * lambda)}.
#' It is realized via the \code{log()} function. To avoid negative values the
#' minimum +1 of the dataset is added and after the prediction later substracted
#' \code{exp(predict(...)) - min(data) - 1}.
#'
#' The 'log2' fit is fitted via the \code{nls()} function for direct fitting of
#' the model \code{y ~ exp(gamma.0 + gamma.1 * lambda)}. As starting values the
#' results of a LS-fit to a linear model with a log transformed response are used.
#' If \code{nls} does not converge, the model with the starting values is returned.
#'
#' \code{refit()} refits the object with a different extrapolation function.
#'
#' @references
#'Küchenhoff, H., Mwalili, S. M. and Lesaffre, E. (2006) A general method for
#'dealing with misclassification in regression: The Misclassification SIMEX.
#'\emph{Biometrics}, \bold{62}, 85 -- 96
#'
#' Küchenhoff, H., Lederer, W. and E. Lesaffre. (2006) Asymptotic Variance Estimation
#' for the Misclassification SIMEX.
#' \emph{Computational Statistics and Data Analysis}, \bold{51}, 6197 -- 6211
#'
#' Lederer, W. and Küchenhoff, H. (2006) A short introduction to the SIMEX and MCSIMEX. \emph{R News}, \bold{6(4)}, 26--31
#'
#' @author Wolfgang Lederer, \email{wolfgang.lederer@gmail.com}
#' @seealso \code{\link[simex]{misclass}}, \code{\link[simex]{simex}}
#'
#' @examples
#' x <- rnorm(200, 0, 1.142)
#' z <- rnorm(200, 0, 2)
#' y <- factor(rbinom(200, 1, (1 / (1 + exp(-1 * (-2 + 1.5 * x -0.5 * z))))))
#' Pi <- matrix(data = c(0.9, 0.1, 0.3, 0.7), nrow = 2, byrow = FALSE)
#' dimnames(Pi) <- list(levels(y), levels(y))
#' ystar <- misclass(data.frame(y), list(y = Pi), k = 1)[, 1]
#' naive.model <- glm(ystar ~ x + z, family = binomial, x = TRUE, y = TRUE)
#' true.model <- glm(y ~ x + z, family = binomial)
#' simex.model <- mcsimex(naive.model, mc.matrix = Pi, SIMEXvariable = "ystar")
#'
#' op <- par(mfrow = c(2, 3))
#' invisible(lapply(simex.model$theta, boxplot, notch = TRUE, outline = FALSE,
#' names = c(0.5, 1, 1.5, 2)))
#' plot(simex.model)
#' simex.model2 <- refit(simex.model, "line")
#' plot(simex.model2)
#' par(op)
#'
#' # example using polr from the MASS package
#' \dontrun{
#' if(require(MASS)) {
#' yord <- cut((1 / (1 + exp(-1 * (-2 + 1.5 * x -0.5 * z)))), 3, ordered=TRUE)
#' Pi3 <- matrix(data = c(0.8, 0.1, 0.1, 0.2, 0.7, 0.1, 0.1, 0.2, 0.7), nrow = 3, byrow = FALSE)
#' dimnames(Pi3) <- list(levels(yord), levels(yord))
#' ystarord <- misclass(data.frame(yord), list(yord = Pi3), k = 1)[, 1]
#' naive.ord.model <- polr(ystarord ~ x + z, Hess = TRUE)
#' simex.ord.model <- mcsimex(naive.ord.model, mc.matrix = Pi3,
#' SIMEXvariable = "ystarord", asymptotic=FALSE)
#' }
#' }
#'
#' # example for a function which can be supplied to the function mcsimex()
#' # "ystar" is the variable which is to be misclassified
#' # using the example above
#' \dontrun{
#' my.misclass <- function (datas, k) {
#' ystar <- datas$"ystar"
#' p1 <- matrix(data = c(0.75, 0.25, 0.25, 0.75), nrow = 2, byrow = FALSE)
#' colnames(p1) <- levels(ystar)
#' rownames(p1) <- levels(ystar)
#' p0 <- matrix(data = c(0.8, 0.2, 0.2, 0.8), nrow = 2, byrow = FALSE)
#'
#' colnames(p0) <- levels(ystar)
#' rownames(p0) <- levels(ystar)
#' ystar[datas$x < 0] <-
#' misclass(data.frame(ystar = ystar[datas$x < 0]), list(ystar = p1), k = k)[, 1]
#' ystar[datas$x > 0] <-
#' misclass(data.frame(ystar = ystar[datas$x > 0]), list(ystar = p0), k = k)[, 1]
#' ystar <- factor(ystar)
#' return(data.frame(ystar))}
#'
#' simex.model.differential <- mcsimex(naive.model, mc.matrix = "my.misclass", SIMEXvariable = "ystar")
#' }
#'
#' @keywords models
#' @export
mcsimex <- function(model,
SIMEXvariable,
mc.matrix,
lambda = c(0.5, 1,1.5, 2),
B = 100,
fitting.method = "quadratic",
jackknife.estimation = "quadratic",
asymptotic = TRUE) {
fitting.method <- substr(fitting.method, 1, 4)
if (!any(fitting.method == c("quad", "line", "nonl", "logl", "log2"))) {
warning("Fitting Method not implemented. Using: quadratic", call. = FALSE)
fitting.method <- "quad"
}
if (jackknife.estimation != FALSE)
jackknife.estimation <- substr(jackknife.estimation, 1, 4)
if (!any(jackknife.estimation == c("quad", "line", "nonl", "logl", FALSE))) {
warning("Fitting Method (jackknife) not implemented. Using: quadratic",
call. = FALSE)
jackknife.estimation <- "quad"
}
if (any(lambda <= 0)) {
warning("lambda should not contain 0 or negative values. 0 or negative values will be ignored",
call. = FALSE)
lambda <- lambda[lambda >= 0]
}
if (class(model)[1] == "polr" && !any(names(model) == "Hessian"))
stop("The option Hessian must be enabled in the naive model", call. = FALSE)
if (class(model)[1] == "polr" && asymptotic)
stop("Asymptotic estimation is not supported for polr models", call. = FALSE)
if (class(model)[1] == "coxph" && asymptotic)
stop("Asymptotic estimation is not supported for coxph models", call. = FALSE)
if (class(model)[1] == "coxph" && is.null(model$model))
stop("The option model = TRUE must be enabled for coxph models", call. = FALSE)
if (class(model)[1] == "coxph" && grep("Surv\\(", names(model$model)[1]) == 1){
timeEventMatrix <- as.matrix(model$model[[1]])
timeName <- sub("Surv\\(","",strsplit(names(model$model)[1], ", ")[[1]][1])
eventName <- sub("\\)","",strsplit(names(model$model)[1], ", ")[[1]][2])
colnames(timeEventMatrix) <- c(timeName, eventName)
model$model <- cbind(model$model, timeEventMatrix)
}
if (!any(names(model) == "x") && asymptotic && class(model)[1] != "polr")
stop("The option x must be enabled in the naive model for asymptotic variance estimation",
call. = FALSE)
if (is.matrix(mc.matrix)) {
mc.matrix <- list(mc.matrix)
names(mc.matrix) <- SIMEXvariable
}
if (is.list(mc.matrix)) {
if (!all(check.mc.matrix(mc.matrix)))
stop("mc.matrix may contain negative values for exponents smaller than 1")
if (length(mc.matrix) != length(SIMEXvariable))
stop("mc.matrix and SIMEXvariable do not match")
}
if (any(!sapply(as.data.frame(model$model), is.factor)[SIMEXvariable]))
stop("SIMEXvariable must be a factor")
cl <- match.call()
ncoef <- length(model$coefficients)
if (class(model)[1] == "polr")
ncoef <- ncoef + length(model$zeta)
ndes <- length(model$y)
nlambda <- length(lambda)
if (class(model)[1] == "polr")
p.names <- c(names(coef(model)), names(model$zeta))
else
p.names <- names(coef(model))
factors <- lapply(SIMEXvariable, levels)
estimates <- matrix(data = NA, length(lambda) + 1, ncoef)
theta <- matrix(data = NA, B, ncoef)
colnames(theta) <- p.names
theta.all <- vector(mode = "list", nlambda)
if (jackknife.estimation != FALSE) {
var.exp <- list()
var.exp[[1]] <- extract.covmat(model)
}
if (asymptotic) {
psi <- matrix(rep(0, ndes * ncoef), ncol = ncoef, nrow = ndes)
psi <- residuals(model, type = "response") * model$x
PSI <- psi
am <- list()
xi <- model$x
a <- list()
dh <- model$family$mu.eta(model$linear.predictors)
for (k in 1:ndes) a[[k]] <- dh[k] * xi[k, ] %*% t(xi[k, ])
a.mat <- matrix(unlist(a), nrow = length(a), byrow = TRUE)
ab <- matrix(colSums(a.mat), nrow = NROW(a[[1]]), byrow = FALSE)
am[[1]] <- -ab/ndes
}
if (class(model)[1] == "polr")
estimates[1, ] <- c(model$coefficients, model$zeta)
else
estimates[1, ] <- model$coefficients
for (i in 1:length(lambda)) {
if (jackknife.estimation != FALSE)
variance.est <- matrix(0, ncol = ncoef, nrow = ncoef)
if (asymptotic) {
psi <- matrix(0, ncol = ncoef, nrow = ndes)
a <- list()
for (k in 1:ndes) a[[k]] <- matrix(0, nrow = ncoef, ncol = ncoef)
}
for (j in 1:B) {
SIMEXdata <- data.frame(model$model)
# doing the misclassification
SIMEXv <- data.frame(SIMEXdata[, SIMEXvariable])
colnames(SIMEXv) <- SIMEXvariable
if (is.character(mc.matrix)) {
SIMEXdata[, SIMEXvariable] <- eval(call(mc.matrix, SIMEXdata, lambda[i]))
} else {
SIMEXdata[, SIMEXvariable] <- misclass(SIMEXv, mc.matrix, lambda[i])
}
# updating the model and calculating the estimates
model.SIMEX <- update(model, data = data.frame(SIMEXdata))
if (class(model)[1] == "polr")
theta[j, ] <- c(model.SIMEX$coefficients, model.SIMEX$zeta)
else
theta[j, ] <- model.SIMEX$coefficients
if (jackknife.estimation != FALSE) {
variance.est <- variance.est + extract.covmat(model.SIMEX)
}
if (asymptotic) {
xi <- model.SIMEX$x
psi <- psi + (residuals(model.SIMEX, type = "response") * xi)
dh <- model$family$mu.eta(model.SIMEX$linear.predictors)
for (k in 1:ndes) a[[k]] <- a[[k]] - dh[k] * xi[k, ] %*% t(xi[k, ])
}
}
# taking the mean of the estimate -> SIMEX estimate
estimates[i + 1, ] <- colMeans(theta)
theta.all[[i]] <- theta
if (jackknife.estimation != FALSE) {
variance.est <- variance.est/B
s2 <- cov(theta)
var.exp[[i + 1]] <- variance.est - s2
}
if (asymptotic) {
xiB <- psi/B
PSI <- cbind(PSI, xiB)
a.mat <- matrix(unlist(a), nrow = length(a), byrow = TRUE)
ab <- matrix(colSums(a.mat), nrow = NROW(a[[1]]), byrow = FALSE)
am[[i + 1]] <- ab/(B * ndes)
}
}
lambda <- c(0, lambda)
colnames(estimates) <- p.names
# fitting the extrapolation function
switch(fitting.method,
quad = extrapolation <- lm(estimates ~ lambda + I(lambda^2)),
line = extrapolation <- lm(estimates ~ lambda),
logl = extrapolation <- lm(I(log(t(t(estimates) + (abs(apply(estimates, 2, min)) + 1) * (apply(estimates, 2, min) <= 0)))) ~ lambda),
log2 = extrapolation <- fit.logl(lambda, p.names, estimates),
nonl = extrapolation <- fit.nls(lambda, p.names, estimates))
if (any(class(extrapolation) == "lm") && fitting.method == "log2")
fitting.method <- "logl"
# predicting the SIMEX estimate
SIMEX.estimate <- vector(mode = "numeric", length = ncoef)
switch(fitting.method, quad = SIMEX.estimate <- predict(extrapolation, newdata = data.frame(lambda = -1)),
line = SIMEX.estimate <- predict(extrapolation, newdata = data.frame(lambda = -1)),
nonl = for (i in 1:length(p.names)) SIMEX.estimate[i] <- predict(extrapolation[[p.names[i]]], newdata = data.frame(lambda = -1)),
logl = SIMEX.estimate <- exp(predict(extrapolation, newdata = data.frame(lambda = -1))) - (abs(apply(estimates, 2, min)) + 1) * (apply(estimates, 2, min) <= 0),
log2 = for (i in 1:length(p.names)) SIMEX.estimate[i] <- predict(extrapolation[[p.names[i]]], newdata = data.frame(lambda = -1)) - ((abs(apply(estimates, 2, min)) + 1) * (apply(estimates, 2, min) <= 0))[i])
if (jackknife.estimation != FALSE) {
variance.jackknife <- matrix(unlist(var.exp), ncol = ncoef^2, byrow = TRUE)
switch(jackknife.estimation,
quad = extrapolation.variance <- lm(variance.jackknife ~ lambda + I(lambda^2)),
line = extrapolation.variance <- lm(variance.jackknife ~ lambda),
logl = extrapolation.variance <- lm(I(log(t(t(variance.jackknife) + (abs(apply(variance.jackknife, 2, min)) + 1) * (apply(variance.jackknife, 2, min) <= 0)))) ~ lambda),
nonl = extrapolation.variance <- fit.nls(lambda, 1:NCOL(variance.jackknife), variance.jackknife))
variance.jackknife2 <- vector("numeric", ncoef^2)
switch(jackknife.estimation,
nonl = for (i in 1:NCOL(variance.jackknife)) variance.jackknife2[i] <- predict(extrapolation.variance[[i]], newdata = data.frame(lambda = -1)),
quad = variance.jackknife2 <- predict(extrapolation.variance, newdata = data.frame(lambda = -1)),
line = variance.jackknife2 <- predict(extrapolation.variance, newdata = data.frame(lambda = -1)),
logl = variance.jackknife2 <- exp(predict(extrapolation.variance, newdata = data.frame(lambda = -1))) - (abs(apply(variance.jackknife, 2, min)) + 1) * (apply(variance.jackknife, 2, min) <= 0))
variance.jackknife <- rbind(variance.jackknife2, variance.jackknife)
variance.jackknife.lambda <- cbind(c(-1, lambda), variance.jackknife)
variance.jackknife <- matrix(variance.jackknife[1, ], nrow = ncoef,
ncol = ncoef, byrow = TRUE)
dimnames(variance.jackknife) <- list(p.names, p.names)
}
if (asymptotic) {
c11 <- cov(PSI)
a11 <- diag.block(am)
a11.inv <- solve(a11)
sigma <- a11.inv %*% c11 %*% t(a11.inv)
s <- construct.s(ncoef, lambda, fitting.method, extrapolation)
d.inv <- solve(s %*% t(s))
sigma.gamma <- d.inv %*% s %*% sigma %*% t(s) %*% d.inv
g <- list()
switch(fitting.method,
quad = g <- c(1, -1, 1), line = g <- c(1, -1),
logl = for (i in 1:ncoef) g[[i]] <- c(exp(coef(extrapolation)[1,i] - coef(extrapolation)[2, i]), -exp(coef(extrapolation)[1,i] - coef(extrapolation)[2, i])),
log2 = for (i in 1:ncoef) g[[i]] <- c(exp(coef(extrapolation[[i]])[1] - coef(extrapolation[[i]])[2]), -exp(coef(extrapolation[[i]])[1] - coef(extrapolation[[i]])[2])), nonl = for (i in 1:ncoef) g[[i]] <- c(-1, -(coef(extrapolation[[i]])[3] - 1)^-1, coef(extrapolation[[i]])[2]/(coef(extrapolation[[i]])[3] - 1)^2))
g <- diag.block(g, ncoef)
variance.asymptotic <- (t(g) %*% sigma.gamma %*% g)/ndes
dimnames(variance.asymptotic) <- list(p.names, p.names)
}
# creating class 'mcsimex'
theta <- matrix(unlist(theta.all), nrow = B)
theta.all <- list()
for (i in 1:ncoef) theta.all[[p.names[i]]] <- data.frame(theta[, seq(i, ncoef * nlambda, by = ncoef)])
z <- cbind(lambda, estimates)
z <- rbind(c(-1, SIMEX.estimate), z) # returning the estimated values
colnames(z) <- c("lambda", p.names)
if(class(model)[1] == "polr")
coefs <- z[1, -1][1:length(coef(model))]
else
coefs <- z[1, -1]
erg <- list(coefficients = coefs, SIMEX.estimates = z, lambda = lambda,
model = model, mc.matrix = mc.matrix, B = B, extrapolation = extrapolation,
fitting.method = fitting.method, SIMEXvariable = SIMEXvariable,
call = cl, theta = theta.all)
class(erg) <- ("mcsimex")
if (class(model)[1] == "polr")
erg$zeta <- z[1,-1][(length(coef(model))+1):length(z[1,-1])]
type <- 'response'
if (class(model)[1] == "polr")
type <- 'probs'
if (class(model)[1] == "coxph"){
type <- 'lp' ## Which type of prediction makes sense for residuals?
}
fitted.values <- predict(erg, newdata = model$model[, -1, drop = FALSE],
type = type)
erg$fitted.values <- fitted.values
if (class(model)[1] == "polr" || class(model)[1] == "coxph")
erg$residuals <- NULL
else if (is.factor(model$model[, 1]))
erg$residuals <- as.numeric(levels(model$model[, 1]))[model$model[, 1]] - fitted.values
else erg$residuals <- model$model[, 1] - fitted.values
if (jackknife.estimation != FALSE) {
erg$extrapolation.variance <- extrapolation.variance
erg$variance.jackknife <- variance.jackknife
erg$variance.jackknife.lambda <- variance.jackknife.lambda
}
if (asymptotic) {
erg$PSI <- PSI
erg$c11 <- c11
erg$a11 <- a11
erg$sigma <- sigma
erg$sigma.gamma <- sigma.gamma
erg$g <- g
erg$s <- s
erg$variance.asymptotic <- variance.asymptotic
}
return(erg)
}
#' @describeIn mcsimex Plots of the simulation and extrapolation
#' @export
plot.mcsimex <- function(x,
xlab = expression((1 + lambda)),
ylab = colnames(b)[-1],
ask = FALSE,
show = rep(TRUE, NCOL(b) - 1),
...) {
old.par <- par(no.readonly = TRUE)
on.exit(par(old.par))
par(...)
if (ask)
par(ask = TRUE)
p.names <- names(coef(x))
b <- x$SIMEX.estimates
a <- seq(-1, max(b[, 1]), by = 0.01)
d <- matrix(data = NA, nrow = length(a), ncol = NCOL(b) - 1)
switch(x$fitting.method,
quad = d <- matrix(predict(x$extrapolation, newdata = data.frame(lambda = a)), nrow = length(a), ncol = NCOL(b) - 1),
line = d <- matrix(predict(x$extrapolation, newdata = data.frame(lambda = a)), nrow = length(a), ncol = NCOL(b) - 1),
nonl = for (i in 1:length(p.names)) d[, i] <- predict(x$extrapolation[[p.names[i]]],
newdata = data.frame(lambda = a)),
log2 = for (i in 1:length(p.names)) d[, i] <- predict(x$extrapolation[[p.names[i]]], newdata = data.frame(lambda = a)) -
((abs(apply(x$SIMEX.estimates[-1, -1], 2, min)) + 1) * (apply(x$SIMEX.estimates[-1, -1], 2, min) <= 0))[i],
logl = d <- t(t(exp(predict(x$extrapolation, newdata = data.frame(lambda = a)))) - ((abs(apply(x$SIMEX.estimates[-1,
-1], 2, min)) + 1) * (apply(x$SIMEX.estimates[-1, -1], 2, min) <= 0))))
for (i in 2:NCOL(b)) {
if (show[i - 1]) {
plot(b[, 1] + 1, b[, i], xlab = xlab, ylab = ylab[i - 1], type = "n")
points(b[-1, 1] + 1, b[-1, i], pch = 19)
points(b[1, 1] + 1, b[1, i])
lines(a[a > 0] + 1, d[a > 0, (i - 1)])
lines(a[a < 0] + 1, d[a < 0, (i - 1)], lty = 2)
}
}
}
#' @describeIn mcsimex Predict with mcsimex correction
#' @export
predict.mcsimex <- function(object, newdata, ...) {
new.object <- object$model
new.object$coefficients <- object$coefficients
if (class(new.object)[1] == "polr") {
new.object$zeta <- object$zeta
new.object$fitted.values <- object$fitted.values
}
if (missing(newdata)) {
predict(new.object, ...)
} else {
predict(new.object, newdata = data.frame(newdata), ...)
}
}
#' @describeIn mcsimex Nice printing
#' @export
print.mcsimex <- function(x, digits = max(3, getOption("digits") - 3),
...) {
cat("\nNaive model:\n", deparse(x$model$call), "\n", sep = "")
cat("\nSIMEX-Variables: ")
cat(x$SIMEXvariable, sep = ", ")
cat("\nNumber of Simulations: ", paste(x$B), "\n\n", sep = "")
if (length(coef(x))) {
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2,
quote = FALSE)
} else cat("No coefficients\n")
if (length(x$zeta)) {
cat("Intercepts:\n")
print.default(format(x$zeta, digits = digits), print.gap = 2,
quote = FALSE)
}
cat("\n")
return(invisible(x))
}
#' @describeIn mcsimex Print summary nicely
#' @export
print.summary.mcsimex <- function(x, digits = max(3, getOption("digits") -
3), ...) {
cat("Call:\n")
print(x$call)
cat("\nNaive model: \n")
print(x$naive.model)
cat("\nSimex variable :", x$SIMEXvariable, "\n")
cat("Misclassification matrix: \n")
if (is.character(x$mc.matrix))
print(x$mc.matrix) else lapply(x$mc.matrix, print)
cat("\nNumber of iterations: ", x$B, "\n")
if (!is.null(x$residuals)) {
cat("\nResiduals: \n")
print(summary(x$residuals), digits)
}
cat("\nCoefficients: \n")
if (any(names(x$coefficients) == "asymptotic")) {
cat("\nAsymptotic variance: \n")
printCoefmat(x$coefficients$asymptotic, digits = digits)
}
if (any(names(x$coefficients) == "jackknife")) {
cat("\nJackknife variance: \n")
printCoefmat(x$coefficients$jackknife, digits = digits)
}
return(invisible(x))
}
#' @describeIn mcsimex Summary for mcsimex
#' @export
summary.mcsimex <- function(object, ...) {
if (class(object$model)[1] == "polr")
est <- c(coef(object), object$zeta)
else
est <- coef(object)
p.names <- names(est)
est.table <- list()
if (is.null(resid(object)))
n <- object$model$n
else
n <- length(resid(object))
p <- length(p.names)
rdf <- n - p
if (any(names(object) == "variance.jackknife")) {
se <- sqrt(diag(object$variance.jackknife))
tval <- est/se
pval <- 2 * pt(abs(tval), rdf, lower.tail = FALSE)
est.table[["jackknife"]] <- cbind(est, se, tval, pval)
dimnames(est.table[["jackknife"]]) <- list(p.names, c("Estimate",
"Std. Error", "t value", "Pr(>|t|)"))
}
if (any(names(object) == "variance.asymptotic")) {
se <- sqrt(diag(object$variance.asymptotic))
tval <- est/se
pval <- 2 * pt(abs(tval), rdf, lower.tail = FALSE)
est.table[["asymptotic"]] <- cbind(est, se, tval, pval)
dimnames(est.table[["asymptotic"]]) <- list(p.names, c("Estimate",
"Std. Error", "t value", "Pr(>|t|)"))
}
ans <- list()
class(ans) <- "summary.mcsimex"
ans$coefficients <- est.table
ans$residuals <- resid(object)
ans$call <- object$call
ans$B <- object$B
ans$naive.model <- object$model$call
ans$SIMEXvariable <- object$SIMEXvariable
ans$mc.matrix <- object$mc.matrix
return(ans)
}
#' @describeIn mcsimex Refits the model with a different extrapolation function
#' @export
refit.mcsimex <- function(object, fitting.method = "quadratic", jackknife.estimation = "quadratic",
asymptotic = TRUE, ...)
.refit(object, fitting.method = fitting.method, jackknife.estimation = jackknife.estimation,
asymptotic = asymptotic, ...)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.