robpredict {rsae}

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 class 'pred_model_b'
print(x, digits = max(3L, getOption("digits") - 3L),
    ...)
## S3 method for class 'pred_model_b'
plot(x, type = "e", sort = NULL, ...)
## S3 method for class 'pred_model_b'
residuals(object, ...)
## S3 method for class 'pred_model_b'
as.matrix(x, ...)
## S3 method for class 'pred_model_b'
head(x, n = 6L, ...)
## S3 method for class 'pred_model_b'
tail(x, n = 6L, keepnums = TRUE, addrownums, ...)

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).

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).

Value

An instance of the S3 class pred_model_b

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.

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. doi: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. doi: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. doi: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. doi:10.1002/cjs.10029

See Also

saemodel(), makedata(), fitsaemodel()

Examples

# 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)