Nothing
#' @importFrom Rfast nth colVars
#############################################################################################
#' @name local_influence.mixpoissonreg
#' @title Local Influence Diagnostics for Mixed Poisson Regression Models
#' @description This function provides local influence diagnostic quantities. Currently the conformal normal and normal curvatures are available
#' under several perturbation schemes. The default is the conformal normal curvature since
#' it takes values on \eqn{[0,1]} and other nice properties (see Zhu and Lee, 2001 and Poon and Poon, 1999 for further details).
#' @param model a \code{mixpoissonreg} object.
#' @param perturbation a list or vector of perturbation schemes to be returned. The currently available schemes are
#' "case_weights", "hidden_variable", "mean_explanatory", "precision_explanatory", "simultaneous_explanatory". See Barreto-Souza and Simas (2016) for further details.
#' @param curvature the curvature to be returned, 'conformal' for the conformal normal curvature (see Zhu and Lee, 2001 and Poon and Poon, 1999) or
#' 'normal' (see Zhu and Lee, 2001 and Cook, 1986).
#' @param direction the 'max.eigen' returns the eigenvector associated to the largest eigenvalue of the perturbation matrix. The 'canonical' considers
#' the curvatures under the canonical directions, which is known as "total local curvature" (see Lesaffre and Verbeke, 1998). For conformal
#' normal curvatures both of them coincide. The default is 'canonical'.
#' @param parameters the parameter to which the local influence will be computed. The options are 'all', 'mean' and 'precision'.
#' This argument affects the 'case_weights' and 'hidden_variable' perturbation schemes. The default is 'all'.
#' @param mean.covariates a list or vector of characters containing the mean-explanatory variables to be used in the 'mean-explanatory' and 'simultaneous-explanatory'
#' perturbation schemes. If NULL, the 'mean-explanatory' and 'simultaneous-explanatory' perturbation schemes will be computed by perturbing all
#' mean-related covariates. The default is NULL.
#' @param precision.covariates a list or vector of characters containing the precision-explanatory variables to be used in the 'precision-explanatory'
#' and 'simultaneous-explanatory'
#' perturbation schemes. If NULL, the 'precision-explanatory' and 'simultaneous-explanatory' perturbation schemes will be computed by perturbing all
#' precision-related covariates. The default is NULL.
#' @param ... other graphical arguments to be passed.
#' @return a list containing the resulting perturbation schemes as elements. Each returned element has an attribute 'benchmark', which for
#' the conformal normal curvature, it is computed following Zhu and Lee (2001), and for normal curvature it is computed following Verbeke and ...
#' If the 'direction' is 'max.eigen' the 'benchmark' attribute is NA.
#'
#' The 'mean_explanatory', 'precision_explanatory' and 'simultaneous_explanatory' elements of the list contain an attribute 'covariates' indicating
#' which covariates were used in the perturbation schemes.
#' @details
#' \code{local_influence.mixpoissonreg} provides local influence diagnostics for mixed Poisson regression models for all perturbation schemes considered in
#' Barreto-Souza and Simas (2016), for normal and conformal normal curvatures. Further, it is also provides results for the canonical directions, which is called
#' the total local influence (see Lesaffre and Verbeke, 1998), as well as for the direction of largest curvature, which is the direction of the eigenvector of the
#' perturbation matrix associated to the largest eigenvalue.
#'
#' \code{local_influence_plot.mixpoissonreg} provides a plot of the local influence diagnostics. Each plot corresponds to a perturbation scheme. The first plot considers
#' the 'case-weights' perturbation; the second plot considers the 'hidden-variable' perturbation (which was introduced in Barreto-Souza and Simas, 2016); the third plot
#' considers the mean-explanatory perturbation; the fourth plot considers the precision-explanatory perturbation; the fifth plot considers the simultanous-explanatory perturbation.
#'
#' For both \code{local_influence.mixpoissonreg} and \code{local_influence_plot.mixpoissonreg}, one can select which covariates will be perturbed in the 'mean-explanatory',
#' 'precision-explanatory' and 'simultaneous-explanatory' perturbation schemes. These are chosen in the 'mean.covariates' and 'precision.covariates' arguments.
#'
#' If one considers the total local influence, then Zhu and Lee (2001) provides benchmark for influential observations for all perturbation schemes. These are returned as
#' attributes in the returned list from \code{local_influence.mixpoissonreg}. When using the \code{local_influence_plot.mixpoissonreg}, only points above the benchmark
#' will be displayed. One can also set the option 'draw_benchmark' to TRUE to plot the benchmark line.
#'
#'
#'
#' @references
#' DOI:10.1007/s11222-015-9601-6 \doi{10.1007/s11222-015-9601-6}(Barreto-Souza and Simas; 2016)
#'
#' Cook, R. D. (1986) *Assessment of Local Influence.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 48, pp.133-169. \href{https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01398.x}{https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01398.x}
#'
#' Lesaffre, E. and Verbeke, G. (1998) *Local Influence in Linear Mixed Models*. Biometrics, 54, pp. 570-582.
#'
#' Poon, W.-Y. and Poon, Y.S. (1999) *Conformal normal curvature and assessment of local influence.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 61, pp.51-61. \href{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00162}{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00162}
#'
#' Zhu, H.-T. and Lee, S.-Y. (2001) *Local influence for incomplete data models.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 63, pp.111-126. \href{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00279}{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00279}
#' @examples
#' daysabs_progML <- mixpoissonregML(daysabs ~ prog | prog, data = Attendance)
#' local_influence(daysabs_progML)
#'
#' @export
local_influence.mixpoissonreg <- function(model, perturbation = c("case_weights", "hidden_variable",
"mean_explanatory", "precision_explanatory",
"simultaneous_explanatory"), curvature = c("conformal", "normal"),
direction = c("canonical", "max.eigen"), parameters = c("all", "mean", "precision"),
mean.covariates = NULL, precision.covariates = NULL, ...){
loc_infl <- list()
parameters <- rlang::arg_match(parameters)
direction <- rlang::arg_match(direction)
curvature <- rlang::arg_match(curvature)
modeltype <- model$modeltype
if(is.null(model$x)){
stop("x component not found. fit the model again with argument x = TRUE")
}
if(is.null(model$w)){
stop("w component not found. fit the model again with argument w = TRUE")
}
n = model$nobs
x = model$x
w = model$w
beta = model$coefficients$mean
alpha = model$coefficients$precision
n_beta = length(model$coefficients$mean)
n_alpha = length(model$coefficients$precision)
if(model$intercept[1] & n_beta == 1){
warning("Mean explanatory and simultaneous explanatory should not be considered since
there is only the intercept for the mean.")
}
if(model$intercept[2] & n_alpha == 1){
warning("Precision explanatory and simultaneous explanatory should not be considered since
there is only the intercept for the precision parameter.")
}
mu = model$fitted.values
phi = model$fitted.precisions
if(is.null(model$y)){
y = model$residuals + model$fitted.values
} else{
y = model$y
}
lambda = lambda_r(y,mu,phi,model$modeltype)
a = mu*lambda - y
link_precision <- build_links_mpreg(model$link.precision)
eta_prec <- link_precision$linkfun(phi)
dphideta <- link_precision$mu.eta(eta_prec)
qsi0 <- switch(model$modeltype,
"NB" = {-1},
"PIG" = {-1/2}
)
B<- switch(model$modeltype,
"PIG" = {function(x) {-sqrt((-2*x))}},
"NB" = {function(x) {-log(-x)}}
)
dD<-switch(model$modeltype,
"PIG" = {function(x) {1/(2*x)}},
"NB" = {function(x) {log(x)+1-digamma(x)}}
)
kappa = kappa_r(y,mu,phi, model$modeltype)
b = dphideta*(B(qsi0)-kappa - dD(phi)-qsi0*lambda)
Wbeta <- obs_fisher_weight_matrix_mixpoisson(model, parameters = "mean")
Walpha <- obs_fisher_weight_matrix_mixpoisson(model, parameters = "precision")
A = diag(c(a))
B1 = diag(c(b))
for(pert in perturbation){
switch(pert,
"case_weights" = {
Bmatrix = (1 - (parameters == "precision")) * (A%*%x) %*% solve(t(x)%*%Wbeta%*%x) %*% t(A%*%x) +
(1 - (parameters == "mean")) * (B1%*%w)%*%solve(t(w)%*%Walpha%*%w)%*% t(B1%*%w)
},
"hidden_variable" = {
P1 = diag(-c(mu*lambda))
P2 = diag(c(phi*(lambda*qsi0-kappa)) )
Bmatrix = (1 - (parameters == "precision")) * (P1%*%x) %*% solve(t(x)%*%Wbeta%*%x) %*% t(P1%*%x) +
(1 - (parameters == "mean")) * (P2%*%w)%*%solve(t(w)%*%Walpha%*%w)%*% t(P2%*%w)
},
"mean_explanatory" = {
S_x = sqrt(Rfast::colVars(x))
if(!is.null(mean.covariates)){
expl_mean <- names(model$coefficients$mean) %in% mean.covariates
S_x = S_x * expl_mean
}
if(length(S_x) > 1){
S_x = diag(S_x)
}
ones_nbeta = rep(1,n_beta)
ones_n = rep(1,n)
T1 = A%*%ones_n%*%ones_nbeta%*%S_x+Wbeta%*%x%*%beta%*%ones_nbeta%*%S_x
Bmatrix = T1 %*% solve(t(x)%*%Wbeta%*%x) %*% t(T1)
},
"precision_explanatory" = {
S_w = sqrt(Rfast::colVars(w))
if(!is.null(precision.covariates)){
expl_precision <- names(model$coefficients$precision) %in% precision.covariates
S_w = S_w * expl_precision
}
if(length(S_w) > 1){
S_w = diag(S_w)
}
ones_nalpha = rep(1,n_alpha)
ones_n = rep(1,n)
T2 = B1%*%ones_n%*%ones_nalpha%*%S_w+Walpha%*%w%*%alpha%*%ones_nalpha%*%S_w
Bmatrix = T2 %*% solve(t(w)%*%Walpha%*%w) %*% t(T2)
},
"simultaneous_explanatory" = {
S_x = sqrt(Rfast::colVars(x))
if(!is.null(mean.covariates)){
expl_mean <- names(model$coefficients$mean) %in% mean.covariates
S_x = S_x * expl_mean
}
S_w = sqrt(Rfast::colVars(w))
if(!is.null(precision.covariates)){
expl_precision <- names(model$coefficients$precision) %in% precision.covariates
S_w = S_w * expl_precision
}
if(length(S_x) > 1){
S_x = diag(S_x)
}
if(length(S_w) > 1){
S_w = diag(S_w)
}
ones_nbeta = rep(1,n_beta)
ones_nalpha = rep(1,n_alpha)
ones_n = rep(1,n)
T1 = A%*%ones_n%*%ones_nbeta%*%S_x+Wbeta%*%x%*%beta%*%ones_nbeta%*%S_x
T2 = B1%*%ones_n%*%ones_nalpha%*%S_w+Walpha%*%w%*%alpha%*%ones_nalpha%*%S_w
Bmatrix = T1 %*% solve(t(x)%*%Wbeta%*%x) %*% t(T1) + T2 %*% solve(t(w)%*%Walpha%*%w) %*% t(T2)
}
)
Bmatrix <- switch(curvature,
"conformal" = {Bmatrix/sum(diag(Bmatrix))},
"normal" = Bmatrix)
if(pert == "precision_explanatory" & model$intercept[2] & n_alpha == 1){
loc_infl[[pert]] = rep(NA, n)
} else if(pert == "mean_explanatory" & model$intercept[1] & n_beta == 1){
loc_infl[[pert]] = rep(NA, n)
} else if(pert == "simultaneous_explanatory"){
if( (model$intercept[1] & n_beta == 1) | (model$intercept[2] & n_alpha == 1)){
loc_infl[[pert]] = rep(NA, n)
} else {
loc_infl[[pert]] = switch(direction,
"max.eigen" = {tryCatch(eigen(Bmatrix, symmetric = TRUE)$vec[,1], error = function(e){rep(NA, n)})},
"canonical" = {diag(Bmatrix)})
}
} else {
loc_infl[[pert]] = switch(direction,
"max.eigen" = {tryCatch(eigen(Bmatrix, symmetric = TRUE)$vec[,1], error = function(e){rep(NA, n)})},
"canonical" = {diag(Bmatrix)})
}
names(loc_infl[[pert]]) = 1:n
benchmark = switch(curvature,
"conformal" = {
ifelse(direction == "canonical", 1/n + 2*stats::sd(loc_infl[[pert]]), NA) #Zhu and Lee (2001)
},
"normal" = {
ifelse(direction == "canonical", 2 * mean(loc_infl[[pert]]), NA) # Verbeke and Molenberghs (2000, sect. 11.3)
}
)
attr(loc_infl[[pert]], "benchmark") = benchmark
if(pert == "mean_explanatory"){
attr(loc_infl[[pert]], "covariates") = mean.covariates
if(is.null(attr(loc_infl[[pert]], "covariates"))){
attr(loc_infl[[pert]], "covariates") = "all"
}
}
if(pert == "precision_explanatory"){
attr(loc_infl[[pert]], "covariates") = precision.covariates
if(is.null(attr(loc_infl[[pert]], "covariates"))){
attr(loc_infl[[pert]], "covariates") = "all"
}
}
if(pert == "simultaneous_explanatory"){
attr(loc_infl[[pert]], "covariates") = list("mean" = mean.covariates, "precision" = precision.covariates)
if(is.null(attr(loc_infl[[pert]], "covariates")$mean)){
attr(loc_infl[[pert]], "covariates")$mean = "all"
}
if(is.null(attr(loc_infl[[pert]], "covariates")$precision)){
attr(loc_infl[[pert]], "covariates")$precision = "all"
}
}
}
loc_infl
}
#############################################################################################
#' @name local_influence_plot.mixpoissonreg
#' @title Local Influence Plot Diagnostics for Mixed Poisson Regression Models
#' @description Local influence plots for mixed Poisson regression models. Currently the conformal normal and normal curvatures are available
#' under several perturbation schemes. The default is the conformal normal curvature since
#' it takes values on \eqn{[0,1]} and other nice properties (see Zhu and Lee, 2001 and Poon and Poon, 1999 for further details).
#' @param model a \code{mixpoissonreg} object.
#' @param which a list or vector indicating which plots should be displayed. If a subset of the plots is required, specify a subset of the numbers 1:5, see caption below (and the 'Details') for the different kinds.
#' @param caption captions to appear above the plots; character vector or list of valid graphics annotations. Can be set to "" or NA to suppress all captions.
#' @param sub.caption common title-above the figures if there are more than one. If NULL, as by default, a possible abbreviated version of deparse(x$call) is used.
#' @param curvature the curvature to be returned, 'conformal' for the conformal normal curvature (see Zhu and Lee, 2001 and Poon and Poon, 1999) or
#' 'normal' (see Zhu and Lee, 2001 and Cook, 1986).
#' @param direction the 'max.eigen' returns the eigenvector associated to the largest eigenvalue of the perturbation matrix. The 'canonical' considers
#' the curvatures under the canonical directions, which is known as "total local curvature" (see Lesaffre and Verbeke, 1998). For conformal
#' normal curvatures both of them coincide. The default is 'canonical'.
#' @param parameters the parameter to which the local influence will be computed. The options are 'all', 'mean' and 'precision'.
#' This argument affects the 'case_weights' and 'hidden_variable' perturbation schemes. The default is 'all'.
#' @param mean.covariates a list or vector of characters containing the mean-explanatory variables to be used in the 'mean-explanatory' and 'simultaneous-explanatory'
#' perturbation schemes. If NULL, the 'mean-explanatory' and 'simultaneous-explanatory' perturbation schemes will be computed by perturbing all
#' mean-related covariates. The default is NULL.
#' @param precision.covariates a list or vector of characters containing the precision-explanatory variables to be used in the 'precision-explanatory'
#' and 'simultaneous-explanatory'
#' perturbation schemes. If NULL, the 'precision-explanatory' and 'simultaneous-explanatory' perturbation schemes will be computed by perturbing all
#' precision-related covariates. The default is NULL.
#' @param detect.influential logical. Indicates whether the benchmark should be used to detect influential observations and identify them on the plot. If there is no benchmark available,
#' the top 'n.influential' observations will be identified in the plot by their indexes.
#' @param n.influential interger. The maximum number of influential observations to be identified on the plot.
#' @param draw.benchmark logical. Indicates whether a horizontal line identifying the benchmark should be drawn.
#' @param lty.benchmark the line type of the benchmark if drawn.
#' @param type_plot what type of plot should be drawn. The default is 'h'.
#' @param ask logical; if \code{TRUE}, the user is asked before each plot.
#' @param main character; title to be placed at each plot additionally (and above) all captions.
#' @param labels.id vector of labels, from which the labels for extreme points will be chosen. The default uses the observation numbers.
#' @param cex.id magnification of point labels.
#' @param cex.caption controls the size of caption.
#' @param cex.oma.main controls the size of the sub.caption only if that is above the figures when there is more than one.
#' @param include.modeltype logical. Indicates whether the model type ('NB' or 'PIG') should be displayed on the captions.
#' @param ... other graphical arguments to be passed.
#' @return Called for its side effects.
#' @details
#' \code{local_influence.mixpoissonreg} provides local influence diagnostics for mixed Poisson regression models for all perturbation schemes considered in
#' Barreto-Souza and Simas (2016), for normal and conformal normal curvatures. Further, it is also provides results for the canonical directions, which is called
#' the total local influence (see Lesaffre and Verbeke, 1998), as well as for the direction of largest curvature, which is the direction of the eigenvector of the
#' perturbation matrix associated to the largest eigenvalue.
#'
#' \code{local_influence_plot.mixpoissonreg} provides a plot of the local influence diagnostics. Each plot corresponds to a perturbation scheme. The first plot considers
#' the 'case-weights' perturbation; the second plot considers the 'hidden-variable' perturbation (which was introduced in Barreto-Souza and Simas, 2016); the third plot
#' considers the mean-explanatory perturbation; the fourth plot considers the precision-explanatory perturbation; the fifth plot considers the simultanous-explanatory perturbation.
#'
#' For both \code{local_influence.mixpoissonreg} and \code{local_influence_plot.mixpoissonreg}, one can select which covariates will be perturbed in the 'mean-explanatory',
#' 'precision-explanatory' and 'simultaneous-explanatory' perturbation schemes. These are chosen in the 'mean.covariates' and 'precision.covariates' arguments.
#'
#' If one considers the total local influence, then Zhu and Lee (2001) provides benchmark for influential observations for all perturbation schemes. These are returned as
#' attributes in the returned list from \code{local_influence.mixpoissonreg}. When using the \code{local_influence_plot.mixpoissonreg}, only points above the benchmark
#' will be displayed. One can also set the option 'draw_benchmark' to TRUE to plot the benchmark line.
#' @references
#' DOI:10.1007/s11222-015-9601-6 \doi{10.1007/s11222-015-9601-6}(Barreto-Souza and Simas; 2016)
#'
#' Cook, R. D. (1986) *Assessment of Local Influence.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 48, pp.133-169. \href{https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01398.x}{https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01398.x}
#'
#' Lesaffre, E. and Verbeke, G. (1998) *Local Influence in Linear Mixed Models*. Biometrics, 54, pp. 570-582.
#'
#' Poon, W.-Y. and Poon, Y.S. (1999) *Conformal normal curvature and assessment of local influence.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 61, pp.51-61. \href{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00162}{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00162}
#'
#' Zhu, H.-T. and Lee, S.-Y. (2001) *Local influence for incomplete data models.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 63, pp.111-126. \href{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00279}{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00279}
#' @examples
#' daysabs_progML <- mixpoissonregML(daysabs ~ prog | prog, data = Attendance)
#' local_influence_plot(daysabs_progML)
#'
#' @export
local_influence_plot.mixpoissonreg <- function(model, which = c(1,2,3,4),
caption = list("Case Weights Perturbation",
"Hidden Variable Perturbation",
"Mean Explanatory Perturbation",
"Precision Explanatory Perturbation",
"Simultaneous Explanatory Perturbation"),
sub.caption = NULL,
detect.influential = TRUE, n.influential = 5,
draw.benchmark = FALSE, lty.benchmark = 2,
type_plot = "h",
curvature = c("conformal", "normal"),
direction = c("canonical", "max.eigen"), parameters = c("all", "mean", "precision"),
mean.covariates = NULL, precision.covariates = NULL, main = "",
ask = prod(graphics::par("mfcol")) <
length(which) && grDevices::dev.interactive(),
labels.id = names(stats::residuals(model)),
cex.id = 0.75,
cex.oma.main = 1.25,
cex.caption = 1,
include.modeltype = TRUE,
...){
direction <- rlang::arg_match(direction)
curvature <- rlang::arg_match(curvature)
n_beta = length(model$coefficients$mean)
n_alpha = length(model$coefficients$precision)
if(model$intercept[1] & n_beta == 1 & any(c(3,5)%in% which)){
warning("Removing mean explanatory and simultaneous explanatory perturbations since
there is only the intercept for the mean.")
which <- setdiff(which, c(3,5))
}
if(model$intercept[2] & n_alpha == 1 & any(c(4,5)%in% which)){
warning("Removing precision explanatory and simultaneous explanatory perturbations since
there is only the intercept for the precision parameter.")
which <- setdiff(which, c(4,5))
}
getCaption <- function(k) if (length(caption) < k)
NA_character_
else {
if(include.modeltype){
grDevices::as.graphicsAnnot(paste0(caption[[k]], " - ", model$modeltype, " Regression"))
} else {
grDevices::as.graphicsAnnot(caption[[k]])
}
}
if (is.null(sub.caption)) {
cal <- model$call
if (!is.na(m.f <- match("formula", names(cal)))) {
cal <- cal[c(1, m.f)]
names(cal)[2L] <- ""
}
cc <- deparse(cal, 80)
nc <- nchar(cc[1L], "c")
abbr <- length(cc) > 1 || nc > 75
sub.caption <- if (abbr)
paste(substr(cc[1L], 1L, min(75L, nc)), "...")
else cc[1L]
}
ylab_infl <- switch(curvature,
"conformal" = {
yl <- switch(direction,
"canonical" = "Total Local Influence (Conformal)",
"max.eigen" = "Largest Curvature Direction (Conformal)")
yl
},
"normal" = {
yl <- switch(direction,
"canonical" = "Total Local Influence (Normal)",
"max.eigen" = "Largest Curvature Direction (Normal)")
}
)
one.fig <- prod(graphics::par("mfcol")) == 1
loc_infl <- suppressWarnings(local_influence(model, curvature = curvature,
direction = direction, parameters = parameters,
mean.covariates = mean.covariates,
precision.covariates = precision.covariates))
if (ask) {
oask <- grDevices::devAskNewPage(TRUE)
on.exit(grDevices::devAskNewPage(oask))
}
pert <- c("case_weights", "hidden_variable",
"mean_explanatory", "precision_explanatory",
"simultaneous_explanatory")
for(i in 1:length(pert)){
if(i %in% which){
xlab <- paste0("Index\n ", sub.caption)
if(detect.influential) {
ylim <- range(loc_infl[[pert[i]]], na.rm = TRUE)
ylim <- grDevices::extendrange(r = ylim, f = 0.08)
}
graphics::plot(loc_infl[[pert[i]]], type = type_plot, main = main, ylab = ylab_infl, ylim = ylim, ...)
graphics::mtext(getCaption(i), side = 3, cex = cex.caption)
if (one.fig)
graphics::title(sub = sub.caption, ...)
if(detect.influential){
bm <- attr(loc_infl[[pert[i]]], "benchmark")
infl_points <- as.vector(Rfast::nth(abs(loc_infl[[pert[i]]]), k = n.influential,
num.of.nths = n.influential,
index.return = TRUE, descending = TRUE))
if(!is.na(bm)){
infl_points <- (loc_infl[[pert[i]]][infl_points] > bm)
idx_x <- as.integer(names(which(infl_points == TRUE)))
idx_y <- loc_infl[[pert[i]]][idx_x]
graphics::text(idx_x, idx_y, labels = labels.id[idx_x], cex = cex.id, xpd = TRUE, pos = 3, offset = 0.2)
} else{
idx_x_pos <- infl_points[which(loc_infl[[pert[i]]][infl_points] >= 0)]
idx_x_neg <- setdiff(infl_points, idx_x_pos)
idx_y_pos <- loc_infl[[pert[i]]][idx_x_pos]
idx_y_neg <- loc_infl[[pert[i]]][idx_x_neg]
if(length(idx_x_pos)>0){
graphics::text(idx_x_pos, idx_y_pos, labels = labels.id[idx_x_pos], cex = cex.id, xpd = TRUE, pos = 3, offset = 0.2)
}
if(length(idx_x_neg)>0){
graphics::text(idx_x_neg, idx_y_neg, labels = labels.id[idx_x_neg], cex = cex.id, xpd = TRUE, pos = 1, offset = 0.2)
}
}
}
if(draw.benchmark){
bm <- attr(loc_infl[[pert[i]]], "benchmark")
if(!is.na(bm)){
graphics::abline(a = bm, b = 0, lty = lty.benchmark)
}
}
grDevices::dev.flush()
}
}
if (!one.fig && graphics::par("oma")[3L] >= 1)
graphics::mtext(sub.caption, outer = TRUE, cex = 1.25)
invisible()
}
#############################################################################################
#' @name local_influence
#' @aliases local_influence local_influence_plot
#' @title Local Influence Diagnostics
#' @param model an object for which the local influence is desired
#' @param ... further arguments passed to or from other methods.
#' @return The *local_influence* method returns a list containing the resulting perturbation schemes as elements.
#' The *local_influence_plot* is called for its side effects.
#' @details
#' \code{local_influence} is a generic function to return local influence diagnostics under different perturbation schemes and different directions.
#' \code{local_influence_plot} is a generic function to provide friendly plots of such diagnostics.
#'
#' Local influence diagnostics were first introduced by Cook (1986), where several perturbation schemes were introduced and normal curvatures were obtained. Poon and Poon (1999)
#' introduced the conformal normal curvature, which has nice properties and takes values on the unit interval \eqn{[0,1]}. Zhu and Lee (2001) following Cook (1986) and Poon and Poon (1999)
#' introduced normal and conformal normal curvatures for EM-based models.
#' @references
#' Cook, R. D. (1986) *Assessment of Local Influence.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 48, pp.133-169. \href{https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01398.x}{https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01398.x}
#'
#' Poon, W.-Y. and Poon, Y.S. (1999) *Conformal normal curvature and assessment of local influence.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 61, pp.51-61. \href{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00162}{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00162}
#'
#' Zhu, H.-T. and Lee, S.-Y. (2001) *Local influence for incomplete data models.* Journal of the Royal Statistical Society. Series B (Methodological), Vol. 63, pp.111-126. \href{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00279}{https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00279}
#' @seealso \code{\link{local_influence.mixpoissonreg}}, \code{\link{local_influence_plot.mixpoissonreg}},
#' \code{\link{local_influence_autoplot.mixpoissonreg}}
#'
#' @rdname local_influence
#' @export
local_influence <- function(model, ...){
UseMethod("local_influence", model)
}
#' @rdname local_influence
#' @export
local_influence_plot <- function(model, ...){
UseMethod("local_influence_plot", model)
}
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.