robpredict: Robust Prediction of Random Effects, Fixed Effects, and Area-Specific Means

Description

Function robpredict() predicts the area-level means by (1) the empirical best linear unbiased predictor (EBLUP) or (2) a robust prediction method which is due to Copt and Victoria-Feser (2009). In addition, the function computes the mean square prediction error (MSPE) of the predicted area-level means by a parametric bootstrap method.

Usage

robpredict(fit, areameans = NULL, k = NULL, reps = NULL, seed = 1024,
    progress_bar = TRUE)
# S3 method for pred_model_b
print(x, digits = max(3L, getOption("digits") - 3L),
    ...)
# S3 method for pred_model_b
plot(x, type = "e", sort = NULL, ...)
# S3 method for pred_model_b
residuals(object, ...)
# S3 method for pred_model_b
as.matrix(x, ...)
# S3 method for pred_model_b
head(x, n = 6L, ...)
# S3 method for pred_model_b
tail(x, n = 6L, keepnums = TRUE, addrownums, ...)

Value

An instance of the S3 class pred_model_b

Arguments

fit: an object of class fit_model_b; a fitted SAE model.
areameans: [matrix] or [NULL] area-level means of dimension (g, p), where g and p denote, respectively, the number of areas and number of fixed-effects terms in the regression model (incl. intercept). By default, areadata = NULL which implies that the predictions are base on the data used in fitting the model (not new data).
k: [numeric] or [NULL] robustness tuning constant (of the Huber psi-function) for robust prediction. Note that k does not necessarily have to be the same as the k that has been used in fitsaemodel(). By default, k = NULL which implies that the tuning constant specified in fitsaemodel() is used.
reps: [integer] or [NULL] number of bootstrap replicates for the computation of the mean squared prediction error (MSPE). If reps = NULL the MSPE is not computed.
seed: [integer] a positive integer used as argument seed in set.seed() to specify the random seed.
progress_bar: [logical] whether a progress bar is displayed for the parametric bootstrap; see NOTE below.
x: an object of class "pred_model_b".
digits: [integer] number of digits to be printed by.
type: [character] type of plot method: "e" (error bars; default) or "l" (lines).
sort: [character] or [NULL] if sort = "means", the predicted means are plotted in ascending order (default: sort = NULL); similarly, with sort = "fixef" and sort = "ranef" the predicted means are sorted along the fixed effects or the random effects, respectively.
object: an object of class fit_model_b.
n: [integer] vector of length up to dim(x), i.e., number of areas.
keepnums: in each dimension, if no names in that dimension are present, create them using the indices included in that dimension. Ignored if dim(x) is NULL or its length 1.
addrownums: deprecated - keepnums should be used instead.
...: additional arguments (not used).

NOTE

Users of Rgui.exe on Windows are recommended to call robpredict() with argument progress_bar = FALSE because Rgui.exe does not handle calls to txtProgressBar() well (the execution time of the same job increases and it tends to stall the execution of R). Users of R-Studio and Rterm.exe are not affected.

Details

Function robpredict() computes predictions of the area-level means and---if required---an estimate of the area-specific mean square prediction error (MSPE).

Prediction of area-level means

Case 1: If areameans is a matrix with area-level means of the explanatory variables, then the computation of the fixed effects effects are based on areameans.
Case 2: If areameans = NULL, then the predictions are based on the sample data that have been used to fit the model.

Mean square prediction error

If reps = NULL, the number of bootstrap replicates is not specified; hence, MSPE is not computed.
If reps is a positive integer and areameans is not NULL (see Case 1 above), then a (robust) parametric bootstrap estimate of MSPE is computed as proposed by Sinha and Rao (2009); see also Lahiri (2003) and Hall.

Robustness

The EBLUP obtains if k = NULL, i.e., if the robustness tuning constant k is unspecified.
Robust predictions of the area-level means are computed if k is a nonnegative real number. Small values of k imply that outliers are heavily downweighted; formally, the EBLUP corresponds to choosing the tuning constant k equal to infinity. The value of the tuning constant k specified in robpredict() can be different from the tuning constant k used in fitting the model. The robust prediction method is due to Copt and Victoria-Feser (2009); see also Heritier et al. (2009, 113-114) and differs from the method in Sinha and Rao (2009).

References

Copt, S. and M.-P. Victoria-Feser (2009). Robust Predictions in Mixed Linear Models, Research Report, University of Geneva.

Lahiri, P. (2003). On the impact of bootstrap in survey sampling and small area estimation. Statistical Science 18, 199--210. tools:::Rd_expr_doi("https://doi.org/10.1214/ss/1063994975")

Hall, P. and T. Maiti (2006). On parametric bootstrap methods for small area prediction. Journal of the Royal Statistical Society. Series B 68, 221--238. tools:::Rd_expr_doi("https://doi.org/10.1111/j.1467-9868.2006.00541.x")

Heritier, S., Cantoni, E., Copt, S., and M.-P. Victoria-Feser (2009). Robust methods in biostatistics. New York: John Wiley and Sons.

Schoch, T. (2012). Robust Unit-Level Small Area Estimation: A Fast Algorithm for Large Datasets. Austrian Journal of Statistics 41, 243--265. tools:::Rd_expr_doi("https://doi.org/10.17713/ajs.v41i4.1548")

Sinha, S.K. and J.N.K. Rao (2009). Robust small area estimation. Canadian Journal of Statistics 37, 381--399. tools:::Rd_expr_doi("https://doi.org/10.1002/cjs.10029")

Examples

Run this code

# use the landsat data
head(landsat)

# set up the model
model <- saemodel(formula = HACorn ~ PixelsCorn + PixelsSoybeans,
    area = ~CountyName,
    data = subset(landsat, subset = (outlier == FALSE)))

# summary of the model
summary(model)

# Huber M-estimate with robustness tuning constant k = 2
m <- fitsaemodel("huberm", model, k = 2)
m

# summary of the fitted model/ estimates
summary(m)

# robust prediction of the random effects and the area-level means (robust
# EBLUP) using the counts-specific means (landsat_means)
head(landsat_means)

# for robust prediction, we use the robustness tuning constant 'k = 1.8'
m_predicted <- robpredict(m, landsat_means, k = 1.8)
head(m_predicted)

# extract prediction as matrix
as.matrix(m_predicted)

# extract residuals from the predictions
head(residuals(m_predicted))

# prediction incl. MSPE; parametric bootstrap with only 'reps = 10'
# replications (for demonstration purposes; in practice, 'reps' should be
# considerably larger)
m_predicted_mspe <- robpredict(m, landsat_means, k = 1.8, reps = 10,
                               progress_bar = FALSE)
head(m_predicted_mspe)

Run the code above in your browser using DataLab