Nothing
#' An S4 class that stores the outputs of the fitted IGLS model.
#'
#' An MLwiN model run via the IGLS estimation method is represented by an "mlwinfitIGLS" object
#'
#' @section An instance of the Class:
#' An instance is created by calling function \code{\link{runMLwiN}}.
#'
#' @slot Nobs Computes the number of complete observations.
#' @slot DataLength Total number of cases.
#' @slot Hierarchy For each higher level of a multilevel model, returns the number of units at that level, together with the minimum, mean and maximum number of lower-level units nested within units of the current level.
#' @slot D A vector specifying the type of distribution to be modelled, which can include \code{'Normal'}, \code{'Binomial'} \code{'Poisson'}, \code{'Multinomial'}, \code{'Multivariate Normal'}, or \code{'Mixed'}.
#' @slot Formula A formula object (or a character string) specifying a multilevel model.
#' @slot levID A character string (vector) of the specified level ID(s).
#' @slot contrasts A list of contrast matrices, one for each factor in the model.
#' @slot xlevels A list of levels for the factors in the model.
#' @slot FP Displays the fixed part estimates.
#' @slot RP Displays the random part estimates.
#' @slot FP.cov Displays a covariance matrix of the fixed part estimates.
#' @slot RP.cov Displays a covariance matrix of the random part estimates.
#' @slot elapsed.time Calculates the CPU time used for fitting the model.
#' @slot call The matched call.
#' @slot LIKE The deviance statistic (-2*log(like)).
#' @slot Converged Boolean indicating whether the model has converged
#' @slot Iterations Number of iterations that the model has run for
#' @slot Meth If \code{Meth = 0} estimation method is set to RIGLS. If \code{Meth = 1} estimation method is set to IGLS.
#' @slot residual If \code{resi.store} is \code{TRUE}, then the residual estimates at all levels are returned.
#' @slot data The data.frame that was used to fit the model.
#' @slot nonlinear A character vector specifying linearisation method used. The first element specifies marginal quasi-likelihood linearization (\code{N = 0}) or penalised quasi-likelihood linearization (\code{N = 1}); The second element specifies first (\code{M = 1}) or second (\code{M = 2}) order approximation.
#' @slot version The MLwiN version used to fit the model
#'
#' @author Zhang, Z., Charlton, C.M.J., Parker, R.M.A., Leckie, G., and Browne,
#' W.J. (2016) Centre for Multilevel Modelling, University of Bristol.
#'
#' @seealso
#' \code{\link{runMLwiN}}
#'
#' @examples
#' \dontrun{
#' library(R2MLwiN)
#' # NOTE: if MLwiN not saved in location R2MLwiN defaults to, specify path via:
#' # options(MLwiN_path = 'path/to/MLwiN vX.XX/')
#' # If using R2MLwiN via WINE, the path may look like this:
#' # options(MLwiN_path = '/home/USERNAME/.wine/drive_c/Program Files (x86)/MLwiN vX.XX/')
#'
#' ## Example: tutorial
#' data(tutorial, package = "R2MLwiN")
#'
#' (mymodel <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student),
#' data = tutorial))
#'
#' ##summary method
#' summary(mymodel)
#'
#' ##logLik method
#' logLik(mymodel)
#' }
#'
#' @name mlwinfitIGLS-class
#' @rdname mlwinfitIGLS-class
#' @exportClass mlwinfitIGLS
setClass(Class = "mlwinfitIGLS", representation = representation(version = "character", Nobs = "numeric", DataLength = "numeric",
Hierarchy = "ANY", D = "ANY", Formula = "ANY", levID = "character", contrasts = "list", xlevels = "list",
FP = "numeric", RP = "numeric", RP.cov = "matrix", FP.cov = "matrix", LIKE = "ANY", elapsed.time = "numeric",
call = "ANY", residual = "list", Converged = "logical", Iterations = "numeric", Meth = "numeric",
nonlinear = "numeric", data = "data.frame"))
#' Extract or Replace parts of "mlwinfitIGLS" objects
#' @param x data frame
#' @param i,j elements to extract or replace. For \code{[} and \code{[[}, these are \code{character}.
#' @param drop not used.
#' @param value a suitable replacement value.
#' @rdname extract-methods-igls
setMethod("[", "mlwinfitIGLS", function(x, i, j, drop) {
if (i == "version") {
return(x@version)
}
if (i == "Nobs") {
return(x@Nobs)
}
if (i == "DataLength") {
return(x@DataLength)
}
if (i == "Hierarchy") {
return(x@Hierarchy)
}
if (i == "D") {
return(x@D)
}
if (i == "Formula") {
return(x@Formula)
}
if (i == "levID") {
return(x@levID)
}
if (i == "contrasts") {
return(x@contrasts)
}
if (i == "xlevels") {
return(x@xlevels)
}
if (i == "FP") {
return(x@FP)
}
if (i == "RP") {
return(x@RP)
}
if (i == "FP.cov") {
return(x@FP.cov)
}
if (i == "RP.cov") {
return(x@RP.cov)
}
if (i == "elapsed.time") {
return(x@elapsed.time)
}
if (i == "call") {
return(x@call)
}
if (i == "LIKE") {
return(x@LIKE)
}
if (i == "Converged") {
return(x@Converged)
}
if (i == "Iterations") {
return(x@Iterations)
}
if (i == "Meth") {
return(x@Meth)
}
if (i == "nonlinear") {
return(x@nonlinear)
}
if (i == "residual") {
return(x@residual)
}
if (i == "data") {
return(x@data)
}
})
#' @rdname extract-methods-igls
setReplaceMethod("[", signature(x = "mlwinfitIGLS"), function(x, i, j, value) {
if (i == "version") {
x@version <- value
}
if (i == "Nobs") {
x@Nobs <- value
}
if (i == "DataLength") {
x@DataLength <- value
}
if (i == "Hierarchy") {
x@Hierarchy <- value
}
if (i == "D") {
x@D <- value
}
if (i == "Formula") {
x@Formula <- value
}
if (i == "levID") {
x@levID <- value
}
if (i == "contrasts") {
x@contrasts <- value
}
if (i == "xlevels") {
x@xlevels <- value
}
if (i == "FP") {
x@FP <- value
}
if (i == "RP") {
x@RP <- value
}
if (i == "FP.cov") {
x@FP.cov <- value
}
if (i == "RP.cov") {
x@RP.cov <- value
}
if (i == "elapsed.time") {
x@elapsed.time <- value
}
if (i == "call") {
x@call <- value
}
if (i == "LIKE") {
x@LIKE <- value
}
if (i == "Converged") {
x@Converged <- value
}
if (i == "Iterations") {
x@Iterations <- value
}
if (i == "Meth") {
x@Meth <- value
}
if (i == "nonlinear") {
x@nonlinear <- value
}
if (i == "residual") {
x@residual <- value
}
if (i == "data") {
x@data <- value
}
validObject(x)
return(x)
})
#' @rdname extract-methods-igls
setMethod("[[", "mlwinfitIGLS", function(x, i, j, drop) {
if (i == "version") {
return(x@version)
}
if (i == "Nobs") {
return(x@Nobs)
}
if (i == "DataLength") {
return(x@DataLength)
}
if (i == "Hierarchy") {
return(x@Hierarchy)
}
if (i == "D") {
return(x@D)
}
if (i == "Formula") {
return(x@Formula)
}
if (i == "levID") {
return(x@levID)
}
if (i == "contrasts") {
return(x@contrasts)
}
if (i == "xlevels") {
return(x@xlevels)
}
if (i == "FP") {
return(x@FP)
}
if (i == "RP") {
return(x@RP)
}
if (i == "FP.cov") {
return(x@FP.cov)
}
if (i == "RP.cov") {
return(x@RP.cov)
}
if (i == "elapsed.time") {
return(x@elapsed.time)
}
if (i == "call") {
return(x@call)
}
if (i == "LIKE") {
return(x@LIKE)
}
if (i == "Converged") {
return(x@Converged)
}
if (i == "Iterations") {
return(x@Iterations)
}
if (i == "Meth") {
return(x@Meth)
}
if (i == "nonlinear") {
return(x@nonlinear)
}
if (i == "residual") {
return(x@residual)
}
if (i == "data") {
return(x@data)
}
})
#' @rdname extract-methods-igls
setReplaceMethod("[[", signature(x = "mlwinfitIGLS"), function(x, i, j, value) {
if (i == "version") {
x@version <- value
}
if (i == "Nobs") {
x@Nobs <- value
}
if (i == "DataLength") {
x@DataLength <- value
}
if (i == "Hierarchy") {
x@Hierarchy <- value
}
if (i == "D") {
x@D <- value
}
if (i == "Formula") {
x@Formula <- value
}
if (i == "levID") {
x@levID <- value
}
if (i == "contrasts") {
x@contrasts <- value
}
if (i == "xlevels") {
x@xlevels <- value
}
if (i == "FP") {
x@FP <- value
}
if (i == "RP") {
x@RP <- value
}
if (i == "FP.cov") {
x@FP.cov <- value
}
if (i == "RP.cov") {
x@RP.cov <- value
}
if (i == "elapsed.time") {
x@elapsed.time <- value
}
if (i == "call") {
x@call <- value
}
if (i == "LIKE") {
x@LIKE <- value
}
if (i == "Converged") {
x@Converged <- value
}
if (i == "Iterations") {
x@Iterations <- value
}
if (i == "Meth") {
x@Meth <- value
}
if (i == "nonlinear") {
x@nonlinear <- value
}
if (i == "residual") {
x@residual <- value
}
if (i == "data") {
x@data <- value
}
validObject(x)
return(x)
})
#' Summarize "mlwinfitIGLS" objects
#' @param object an \code{\link{mlwinfitIGLS-class}} object
#' @param ... other parameters
#' @seealso \code{\link[stats4]{summary-methods}}
#' @export
setMethod("summary", signature(object = "mlwinfitIGLS"), function(object, ...) {
object
})
printIGLS <- function(x, digits = max(3, getOption("digits") - 2), signif.stars = getOption("show.signif.stars"),
...) {
object <- summary(x)
align2right <- function(titlename, ele) {
# for printing the table on the screen
all.ele <- c(titlename, ele)
len.all.ele <- nchar(all.ele)
max.len.ele <- max(len.all.ele)
for (j in 1:length(all.ele)) {
if (len.all.ele[j] < max.len.ele) {
len.diff <- max.len.ele - len.all.ele[j]
all.ele[j] <- paste(paste(rep(" ", len.diff), collapse = ""), all.ele[j], sep = "")
}
}
all.ele
}
align2left <- function(titlename, ele) {
# for printing the table on the screen
all.ele <- c(titlename, ele)
len.all.ele <- nchar(all.ele)
max.len.ele <- max(len.all.ele)
for (j in 1:length(all.ele)) {
if (len.all.ele[j] < max.len.ele) {
len.diff <- max.len.ele - len.all.ele[j]
all.ele[j] <- paste(all.ele[j], paste(rep(" ", len.diff), collapse = ""), sep = "")
}
}
all.ele
}
signifstar <- function(pval) {
symnum(pval, corr = FALSE, na = "N/A",
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "** ", "* ", ". ", " "))
}
cat("\n")
cat(paste(rep("-", 50), collapse = "*"), "\n")
cat(object@version, " multilevel model", paste("(", object@D[1], ")", sep = ""), "\n")
if (!is.null(object@Hierarchy))
print(object@Hierarchy)
cat("Estimation algorithm: ")
if (object@Meth == 1) {
cat("IGLS")
} else {
cat("RIGLS")
}
if (object@D[1] != "Normal" && object@D[1] != "Multivariate Normal") {
if (object@nonlinear[1] == 0) {
cat(" MQL")
}
if (object@nonlinear[1] == 1) {
cat(" PQL")
}
cat(object@nonlinear[2])
}
cat(" Elapsed time :", paste(round(object@elapsed.time, 2), "s", sep = ""), "\n")
cat("Number of obs: ", object@Nobs, paste0("(from total ", object@DataLength, ")"))
if (object@Converged) {
cat(" The model converged after", object@Iterations, "iterations.\n")
} else {
cat(" The model did not converge after", object@Iterations, "iterations.\n")
}
cat(paste("Log likelihood: ", round(-0.5 * object@LIKE, 1)), "\n")
cat(paste("Deviance statistic: ", round(object@LIKE, 1)), "\n")
cat(paste(rep("-", 50), collapse = "-"), "\n")
cat("The model formula:\n")
print(formula(object))
levID0 <- object@levID
levID.display <- ""
if (is.na(levID0[length(levID0)])) {
levID0 <- levID0[-length(levID0)]
}
for (i in 1:length(levID0)) {
levID.display <- paste(levID.display, "Level ", length(levID0) + 1 - i, ": ", levID0[i], " ", sep = "")
}
cat(levID.display, "\n")
cat(paste(rep("-", 50), collapse = "-"), "\n")
FP.names <- names(object@FP)
RP.names <- names(object@RP)
est.all <- coef(object)
sd.all <- sqrt(diag(vcov(object)))
zscore.all <- est.all / sd.all
pvalue.all <- 2 * stats::pnorm(abs(zscore.all), lower.tail = FALSE)
confint.all <- confint(object)
stars.all <- signifstar(pvalue.all)
cat("The fixed part estimates: ", "\n")
FP.names2 <- gsub("FP+\\_", "", FP.names)
printcol0 <- align2left(" ", FP.names2)
printcol1 <- align2right("Coef.", format(round(est.all[FP.names], digits), nsmall = digits))
printcol2 <- align2right("Std. Err.", format(round(sd.all[FP.names], digits), nsmall = digits))
printcol3 <- align2right("z", format(round(zscore.all[FP.names], 2), nsmall = 2))
printcol4 <- align2right("Pr(>|z|)", formatC(pvalue.all[FP.names]))
printcol4b <- align2right(" ", stars.all[FP.names])
printcol5 <- align2right("[95% Conf.", format(round(confint.all[FP.names, 1], digits), nsmall = digits))
printcol6 <- align2right("Interval]", format(round(confint.all[FP.names, 2], digits), nsmall = digits))
for (i in 1:(1+length(FP.names2))) {
cat(printcol0[i], " ", printcol1[i], " ", printcol2[i], " ", printcol3[i], " ", printcol4[i], " ", printcol4b[i],
" ", printcol5[i], " ", printcol6[i], "\n")
}
if (signif.stars) {
cat("Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ", "\n")
}
nlev <- length(object@levID)
if (is.na(object@levID[length(object@levID)])) {
mlwinlev <- (nlev - 1):1
levID2 <- levID0
} else {
mlwinlev <- nlev:1
levID2 <- object@levID
}
for (i in 1:length(mlwinlev)) {
RPx.pos <- grep(paste("RP", mlwinlev[i], sep = ""), RP.names)
if (length(RPx.pos) != 0) {
cat(paste(rep("-", 50), collapse = "-"), "\n")
RPx.names <- gsub(paste("RP+", mlwinlev[i], "+\\_", sep = ""), "", RP.names[RPx.pos])
printcol0 <- align2left(" ", RPx.names)
printcol1 <- align2right("Coef.", format(round(est.all[RP.names[RPx.pos]], digits), nsmall = digits))
printcol2 <- align2right("Std. Err.", format(round(sd.all[RP.names[RPx.pos]], digits), nsmall = digits))
cat("The random part estimates at the", levID2[i], "level:", "\n")
for (i in 1:(1+length(RPx.names))) {
cat(printcol0[i], " ", printcol1[i], " ", printcol2[i], "\n")
}
}
}
cat(paste(rep("-", 50), collapse = "*"), "\n")
}
#' Show objects of class "mlwinfitIGLS"
#' @param object an \code{\link{mlwinfitIGLS-class}} object
#' @seealso \code{\link[stats4]{show-methods}}
#' @export
setMethod("show", signature(object = "mlwinfitIGLS"), function(object) printIGLS(object))
updateMLwiN <- function(object, Formula., levID., estoptions., ..., keep.order = TRUE, evaluate = TRUE) {
update.formula2 <- function (old, new, ...)
{
C_updateform <- get("C_updateform", asNamespace("stats"), inherits=FALSE)
tmp <- .Call(C_updateform, stats::as.formula(old), stats::as.formula(new))
out <- stats::formula(terms.formula(tmp, simplify = FALSE))
return(out)
}
environment(update.formula2) <- environment(update.formula)
my.update.formula <- function(old, new, keep.order = TRUE, ...) {
env <- environment(stats::as.formula(old))
tmp <- update.formula2(stats::as.formula(old), stats::as.formula(new))
out <- formula(stats::terms.formula(tmp, simplify = FALSE, keep.order = keep.order))
environment(out) <- env
return(out)
}
if (is.null(newcall <- stats::getCall(object)))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$...
if (length(newcall$Formula))
newcall$Formula <- eval(newcall$Formula)
if (!missing(Formula.)) {
newcall$Formula <- my.update.formula(stats::as.formula(newcall$Formula), Formula., keep.order = keep.order)
}
if (!missing(levID.)) {
newcall$levID <- {
if (length(newcall$levID))
my.update.formula(stats::as.formula(newcall$levID), levID., keep.order = keep.order) else levID.
}
}
if (!missing(estoptions.)) {
newcall$estoptions <- {
if (length(newcall$estoptions))
my.update.formula(stats::as.formula(newcall$estoptions), estoptions., keep.order = keep.order) else estoptions.
}
}
if (length(extras)) {
existing <- !is.na(match(names(extras), names(newcall)))
for (a in names(extras)[existing]) newcall[[a]] <- extras[[a]]
if (any(!existing)) {
newcall <- c(as.list(newcall), extras[!existing])
newcall <- as.call(newcall)
}
}
if (evaluate)
eval(newcall, sys.parent()) else newcall
}
#' Update "mlwinfitIGLS" objects
#' @param object a valid \code{mlwinfitIGLS} class object with an R function call component named \code{call}, the expression used to create itself.
#' @param Formula. changes to the formula. This is a two sided formula where "." is substituted for existing components in the \code{Formula} component of \code{object$call}.
#' @param levID. changes to the specifications of level ID(s).
#' @param estoptions. changes to the specifications of a list of options used for estimating the model.
#' @param ... additional arguments to the call, or arguments with changed values.
#' @param keep.order a logical value indicating whether the terms should keep their positions.
#' @param evaluate if \code{TRUE} (the default) the new call is evaluated;
#' otherwise the call is returned as an unevaluated expression.
#' @return either a new updated \code{mlwinfitIGLS} class object, or else an unevaluated expression for creating such an object.
#' @seealso \code{\link[stats4]{update-methods}}
#' @export
setMethod("update", signature(object = "mlwinfitIGLS"), updateMLwiN)
#' Extract the coefficient vector from "mlwinfitIGLS" objects
#' @param object An \code{\link{mlwinfitIGLS-class}} object
#' @param ... Other arguments
#' @seealso \code{\link[stats4]{coef-methods}}
#' @export
setMethod("coef", signature(object = "mlwinfitIGLS"), function(object, ...) {
c(object@FP, object@RP)
})
#' Extract the approximate variance-covariance matrix from "mlwinfitIGLS" objects
#' @param object An \code{\link{mlwinfitIGLS-class}} object
#' @param ... Other arguments
#' @seealso \code{\link[stats4]{vcov-methods}}
#' @export
setMethod("vcov", signature(object = "mlwinfitIGLS"), function(object, ...) {
m <- matrix(0, nrow(object@FP.cov) + nrow(object@RP.cov), ncol(object@FP.cov) + ncol(object@RP.cov))
colnames(m) <- c(colnames(object@FP.cov), colnames(object@RP.cov))
rownames(m) <- c(rownames(object@FP.cov), rownames(object@RP.cov))
m[colnames(object@FP.cov), rownames(object@FP.cov)] <- object@FP.cov
m[colnames(object@RP.cov), rownames(object@RP.cov)] <- object@RP.cov
m
})
#' Returns the residual degrees-of-freedom extracted from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object.
#' @param ... Other arguments
#' @seealso \code{\link[stats]{nobs}}, \code{\link[stats]{coef}}
#' @method df.residual mlwinfitIGLS
#' @export
df.residual.mlwinfitIGLS <- function(object, ...) {
nobs(object) - length(coef(object))
}
#' Returns the fitted values from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object.
#' @param ... Other arguments.
#' @seealso \code{\link[stats]{fitted.values}}
#' @method fitted mlwinfitIGLS
#' @export
fitted.mlwinfitIGLS <- function(object, ...) {
predict(object, type = "response")
}
#' Returns the residual data from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object
#' @param ... Other arguments.
#' @seealso \code{\link[stats]{residuals}}
#' @method residuals mlwinfitIGLS
#' @export
residuals.mlwinfitIGLS <- function(object, ...) {
form <- Formula.translate(object@Formula, object@D, object@data)
if (!is.list(form$resp)) {
D <- object@D
indata <- object@data
tval <- fitted(object)
if (D[1] == "Binomial") {
tval <- tval * indata[, D[3]]
}
if (D[1] == "Poisson" || D[1] == "Negbinom") {
if (!is.na(D[3])) {
tval <- tval + indata[, D[3]]
}
}
object@data[[form$resp]] - tval
} else {
warning("residuals only implemented for univariate models")
NULL
}
}
#' Returns the predicted data from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object.
#' @param newdata data frame for which to evaluate predictions
#' @param params a character vector specifying the parameters to use in evaluating predictions.
#' If \code{NULL}, \code{names(object[["FP"]])} is used by default.
#' @param type when this has the value \code{"link"} (default) the linear predictor is returned.
#' When \code{type="terms"} each component of the linear predictor is returned seperately. When \code{type="response"} predictions on the scale of the response are returned.
#' @param se.fit logical. When this is \code{TRUE} (not default) standard error estimates are returned for each prediction.
#' @param terms if \code{type="terms"}, which terms (default is all terms), a character vector.
#' @param ... Other arguments.
#' @seealso \code{\link[stats]{predict}}
#' @method predict mlwinfitIGLS
#' @export
predict.mlwinfitIGLS <- function(object, newdata = NULL, params = NULL, type = "link", se.fit = FALSE,
terms = NULL, ...) {
if (is.null(newdata)) {
indata <- object@data
} else {
indata <- Formula.translate(object@Formula, object@D, newdata)$indata
if (!isTRUE("Intercept" %in% indata)) {
indata[["Intercept"]] <- rep(1, nrow(indata))
}
}
if (is.null(params)) {
fp.names <- names(FP <- object@FP)
} else {
fp.names <- params
}
x.names <- sub("FP_", "", fp.names)
if (type == "link") {
tval <- as.vector(as.matrix(indata[x.names]) %*% as.matrix(object@FP[fp.names])[, 1])
if (se.fit) {
# seval <- as.vector(sqrt(diag(as.matrix(indata[x.names]) %*% as.matrix(object@FP.cov[fp.names, fp.names]) %*%
# t(as.matrix(indata[x.names])))))
seval <- as.vector(sqrt(rowSums(as.matrix(indata[x.names]) %*% as.matrix(object@FP.cov[fp.names, fp.names]) *
indata[x.names])))
return(list(fit = tval, se.fit = seval))
} else {
return(tval)
}
}
if (type == "terms") {
if (!is.null(terms)) {
x.names <- terms
fp.names <- paste0("FP_", terms)
}
tval <- as.matrix(t(t(indata[x.names]) * object@FP[fp.names]))
if (se.fit) {
seval <- as.matrix(sqrt(t(t(indata[x.names]^2) * diag(object@FP.cov[fp.names, fp.names]))))
return(list(fit = tval, se.fit = seval))
} else {
return(tval)
}
}
if (type == "response") {
tval <- as.vector(as.matrix(indata[x.names]) %*% as.matrix(object@FP[fp.names])[, 1])
names(tval) <- 1:(length(tval))
D <- object@D
if (D[1] == "Normal" || D[1] == "Multivariate Normal") {
return(tval)
}
if (D[1] == "Binomial") {
if (D[2] == "logit") {
antilogit <- function(x) {
exp(x)/(1 + exp(x))
}
return(antilogit(tval) * indata[, D[3]])
}
if (D[2] == "probit") {
return(stats::pnorm(tval) * indata[, D[3]])
}
if (D[2] == "cloglog") {
anticloglog <- function(x) {
1 - exp(-exp(x))
}
return(anticloglog(tval) * indata[, D[3]])
}
}
if (D[1] == "Poisson") {
if (is.na(D[3])) {
return(exp(tval))
} else {
return(exp(tval + indata[, D[3]]))
}
}
if (D[1] == "Negbinom") {
if (is.na(D[3])) {
return(exp(tval))
} else {
return(exp(tval + indata[, D[3]]))
}
}
if (D[1] == "Mixed") {
}
if (D[1] == "Multinomial") {
}
warning("link function transformation not yet implemented")
return(NULL)
}
}
#' Returns the log-likelihood from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object.
#' @param ... Other arguments.
#' @seealso \code{\link[stats4]{logLik-methods}}
#' @export
setMethod("logLik", signature(object = "mlwinfitIGLS"), function(object, ...) {
D <- object@D
if (D[1] == "Normal" || D[1] == "Multivariate Normal") {
val <- -0.5 * deviance(object)
attr(val, "df") <- length(coef(object))
attr(val, "nobs") <- nobs(object)
class(val) <- "logLik"
return(val)
} else {
warning("logLik only available for Normal models")
return(NULL)
}
})
#' Returns the deviance from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object
#' @param ... Other arguments
#' @seealso \code{\link[stats]{deviance}}
#' @method deviance mlwinfitIGLS
#' @export
deviance.mlwinfitIGLS <- function(object, ...) {
D <- object@D
if (D[1] == "Normal" || D[1] == "Multivariate Normal") {
return(object@LIKE)
} else {
warning("deviance only available for Normal models")
return(NULL)
}
}
#' Returns the number of used observations from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object.
#' @param ... Other arguments.
#' @seealso \code{\link[stats]{nobs}}
#' @method nobs mlwinfitIGLS
#' @export
nobs.mlwinfitIGLS <- function(object, ...) {
object@Nobs
}
#' Summarize "mlwinfitIGLS" objects
#' @param object an \code{\link{mlwinfitIGLS-class}} object
#' @param ... other parameters
#' @method summary mlwinfitIGLS
#' @export
summary.mlwinfitIGLS <- function(object, ...) {
summary(object)
}
#' Summarize "mlwinfitIGLS" objects
#' @param x an \code{\link{mlwinfitIGLS-class}} object
#' @param digits the number of significant digits to use when printing.
#' @param signif.stars logical. If TRUE, 'significance stars' are printed for each coefficient.
#' @param ... other parameters
#' @seealso \code{\link[base]{print}}
#' @method print mlwinfitIGLS
#' @export
print.mlwinfitIGLS <- function(x, digits = max(3, getOption("digits") - 2), signif.stars = getOption("show.signif.stars"), ...) {
printIGLS(x, digits = digits, signif.stars = signif.stars)
}
#' Summarize "mlwinfitIGLS" objects
#' @param object an \code{\link{mlwinfitIGLS-class}} object
#' @param ... other parameters
#' @seealso \code{\link[methods]{show}}
#' @method show mlwinfitIGLS
#' @export
show.mlwinfitIGLS <- function(object, ...) {
show(object)
}
#' Update "mlwinfitIGLS" objects
#' @param object a valid \code{mlwinfitIGLS} class object with an R function call component named \code{call}, the expression used to create itself.
#' @param ... additional arguments to the call, or arguments with changed values.
#' @return either a new updated \code{mlwinfitIGLS} class object, or else an unevaluated expression for creating such an object.
#' @seealso \code{\link[stats]{update}}
#' @method update mlwinfitIGLS
#' @export
update.mlwinfitIGLS <- function(object, ...) {
update(object)
}
#' Extract the coefficient vector from "mlwinfitIGLS" objects
#' @param object An \code{\link{mlwinfitIGLS-class}} object
#' @param ... Other arguments
#' @seealso \code{\link[stats]{coef}}
#' @method coef mlwinfitIGLS
#' @export
coef.mlwinfitIGLS <- function(object, ...) {
coef(object)
}
#' Extract the approximate variance-covariance matrix from "mlwinfitIGLS" objects
#' @param object An \code{\link{mlwinfitIGLS-class}} object
#' @param ... Other arguments
#' @seealso \code{\link[stats]{vcov}}
#' @method vcov mlwinfitIGLS
#' @export
vcov.mlwinfitIGLS <- function(object, ...) {
vcov(object)
}
#' "mlwinfitIGLS" model formula
#' @param x See \code{\link[stats]{formula}}
#' @param env See \code{\link[stats]{formula}}
#' @param ... Other arguments; see \code{\link[stats]{formula}}
#' @method formula mlwinfitIGLS
#' @export
formula.mlwinfitIGLS <- function(x, env = parent.frame(), ...) {
stats::as.formula(x@Formula)
}
#' Returns the log-likelihood from "mlwinfitIGLS" objects.
#' @param object An \code{\link{mlwinfitIGLS-class}} object.
#' @param ... Other arguments.
#' @seealso \code{\link[stats]{logLik}}
#' @method logLik mlwinfitIGLS
#' @export
logLik.mlwinfitIGLS <- function(object, ...) {
logLik(object)
}
#' Extract coefficients and GOF measures from a statistical object (texreg package).
#' @param model An \code{\link{mlwinfitIGLS-class}} model.
#' @param include.nobs should the number of observations be reported?
#' @param include.loglik should the log-likelihood be reported?
#' @param include.deviance should the deviance be reported?
#' @param ... Other arguments.
#' @seealso \code{\link[texreg]{extract}}
#' @export
setMethod("extract", signature = className("mlwinfitIGLS", "R2MLwiN"), function(model, include.nobs = TRUE, include.loglik=TRUE, include.deviance=TRUE, ...) {
co <- coef(model)
se <- sqrt(diag(vcov(model)))
gof <- numeric()
gof.names <- character()
gof.decimal <- logical()
if (include.nobs == TRUE) {
gof <- c(gof, nobs(model))
gof.names <- c(gof.names, "Num.\\ obs.")
gof.decimal <- c(gof.decimal, FALSE)
}
if (include.loglik == TRUE) {
gof <- c(gof, as.numeric(logLik(model)))
gof.names <- c(gof.names, "Log Likelihood")
gof.decimal <- c(gof.decimal, TRUE)
}
if (include.deviance == TRUE) {
gof <- c(gof, deviance(model))
gof.names <- c(gof.names, "Deviance")
gof.decimal <- c(gof.decimal, TRUE)
}
tr <- texreg::createTexreg(
coef.names = names(co),
coef = co,
se = se,
gof.names = gof.names,
gof = gof,
gof.decimal = gof.decimal
)
return(tr)
})
#' Extract coefficients and GOF measures from a statistical object (memisc package).
#' @param obj An \code{\link{mlwinfitIGLS-class}} model.
#' @param alpha level of the confidence intervals; their coverage should be 1-alpha/2
#' @param ... Other arguments.
#' @seealso \code{\link[memisc]{getSummary}}
#' @method getSummary mlwinfitIGLS
#' @export
getSummary.mlwinfitIGLS <- function (obj, alpha = 0.05, ...) {
co <- t(rbind(coef(obj), sqrt(diag(vcov(obj)))))
z <- co[, 1]/co[,2]
p <- 2 * stats::pnorm(abs(z), lower.tail = FALSE)
ci <- confint(obj, level = 1-alpha)
co <- cbind(co, z, p, ci[, 1], ci[, 2])
colnames(co) <- c("est", "se", "stat", "p", "lwr", "upr")
N <- nobs(obj)
ll <- logLik(obj)
deviance <- deviance(obj)
sumstat <- c(
logLik = ll,
deviance = deviance,
N = N
)
list(coef=co, sumstat=sumstat, contrasts=obj@contrasts, xlevels=obj@xlevels, call=obj@call)
}
#' Summarises information about the components of a model from a statistical object (broom package).
#' @param x An \code{\link{mlwinfitIGLS-class}} model.
#' @param conf.int should the confidence interval be included?
#' @param conf.level confidence interval level
#' @param ... Other arguments.
#' @seealso \code{\link[generics]{tidy}}
#' @method tidy mlwinfitIGLS
#' @export
tidy.mlwinfitIGLS <- function(x, conf.int = FALSE, conf.level = .95, ...) {
est <- coef(x)
term <- names(est)
sd <- sqrt(diag(vcov(x)))
zscore <- est / sd
pval <- 2 * stats::pnorm(abs(zscore), lower.tail = FALSE)
group <- rep("", length(term))
group[grep("FP", term)] <- "fixed"
nlev <- length(x@levID)
if (is.na(x@levID[nlev])) {
mlwinlev <- (nlev - 1):1
} else {
mlwinlev <- nlev:1
}
for (i in 1:length(mlwinlev)) {
group[grep(paste0("RP", i), term)] <- x@levID[mlwinlev[i]]
}
ret <- tibble::tibble(term=term, estimate=est, std.error=sd, statistic=zscore, p.value=pval, group=group)
if (conf.int) {
conf <- confint(x, level = conf.level)
rownames(conf) <- NULL
colnames(conf) <- c("conf.low", "conf.high")
ret <- cbind(ret, conf)
}
ret
}
#' Augment data frame with information derived from the model fit (broom package).
#' @param x An \code{\link{mlwinfitIGLS-class}} model.
#' @param data original data onto which columns should be added
#' @param newdata new data to predict on, optional
#' @param type.predict Type of prediction to compute
#' @param type.residuals Type of residuals to compute
#' @param ... Other arguments.
#' @seealso \code{\link[generics]{augment}}
#' @method augment mlwinfitIGLS
#' @export
augment.mlwinfitIGLS <- function(x, data = x@frame, newdata = NULL, type.predict, type.residuals, ...) {
warning("augment method not yet implemented for mlwinfitIGLS objects")
NULL
}
#' Extract GOF measures from a statistical object (broom package).
#' @param x An \code{\link{mlwinfitIGLS-class}} model.
#' @param ... Other arguments.
#' @seealso \code{\link[generics]{glance}}
#' @method glance mlwinfitIGLS
#' @export
glance.mlwinfitIGLS <- function(x, ...) {
tibble::tibble(
logLik = stats::logLik(x),
AIC = stats::AIC(x),
BIC = stats::BIC(x),
df.residual = stats::df.residual(x),
nobs = stats::nobs(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.