EBLUP: Empirical Best Linear Unbiased Predictor

In qape: Quantile of Absolute Prediction Errors

Empirical Best Linear Unbiased Predictor

Description

The function computes the value of the EBLUP of the linear combination of the variable of interest under the linear mixed model estimated using REML.

Usage

EBLUP(YS, fixed.part, random.part, reg, con, gamma, weights, estMSE)

Arguments

YS	values of the variable of interest observed in the sample.
fixed.part	fixed-effects terms declared as in lmer object.
random.part	random-effects terms declared as in lmer object.
reg	the population matrix of auxiliary variables named in fixed.part and random.part.
con	the population 0-1 vector with 1s for elements in the sample and 0s for elements which are not in the sample.
gamma	the population vector which transpose multiplied by the population vector of the variable of interest gives the predicted characteristic. For example, if gamma is the population vector of 1s, the sum of the values of the variable of interest in the whole dataset is predicted.
weights	the population vector of weights, defined as in lmer object, allowing to include heteroscedasticity of random components in the mixed linear model.
estMSE	logical. If TRUE, the naive MSE estimator and its components are computed.

Details

The function computes the value of the EBLUP of the linear combination of the variable of interest based on the formula (21) in Zadlo (2017) (see Remark 5.1 in the paper for further explanations). Predicted values for unsampled population elements in subsets for which random effects are not observed in the sample are computed based only on fixed effects. The naive MSE estimator of the EBLUP, which is the sum of two components given by equations (31) and (32) in Zadlo (2017) p. 8094, where unknown parameters are replaced by their REML estimates, is also computed. The naive MSE estimator ignores the variability of EBLUP resulting from the estimation of variance components.

Value

The function returns a list with the following objects:

fixed.part	the fixed part of the formula of model.
random.part	the random part of the formula of model.
thetaP	the value of the predictor.
beta	the estimated vector of fixed effects.
Xbeta	the product of two matrices: the population model matrix of auxiliary variables X and the estimated vector of fixed effects.
sigma2R	the estimated variance parameter of the distribution of random components.
R	the estimated covariance matrix of random components for sampled elements.
G	the estimated covariance matrix of random effects.
model	the formula of the model (as in lmer object).
mEst	lmer object with the estimated model.
YS	the sample vector of the variable of interest.
reg	the population matrix of auxiliary variables named in fixed.part and random.part.
con	the population 0-1 vector with 1s for elements in the sample and 0s for elements which are not in the sample.
regS	the sample matrix of auxiliary variables named in fixed.part and random.part.
regR	the matrix of auxiliary variables named in fixed.part and random.part for population elements which are not observed in the sample.
gamma	the population vector which transpose multiplied by the population vector of the variable of interest gives the predicted characteristic.
gammaS	the subvector of gamma for sampled elements.
gammaR	the subvector of gamma for population elements which are not observed in the sample.
weights	the population vector of weights, defined as in lmer object, allowing to include the heteroscedasticity of random components in the mixed linear model.
Z	the population model matrix of auxiliary variables associated with random effects.
ZBlockNames	labels of blocks of random effects in Z matrix.
X	the population model matrix of auxiliary variables associated with fixed effects.
ZS	the submatrix of Z matrix where the number of rows equals the number of sampled elements and the number of columns equals the number of estimated random effects.
XR	the submatrix of X matrix (with the same number of columns) for population elements which are not observed in the sample.
ZR	the submatrix of Z matrix where the number of rows equals the number of population elements which are not observed in the sample and the number of columns equals the number of estimated random effects.
eS	the sample vector of estimated random components.
vS	the estimated vector of random effects.
g1	the first component of the naive MSE estimator (computed if estMSE = TRUE).
g2	the second component of the naive MSE estimator (computed if estMSE = TRUE).
neMSE	the naive MSE estimator (computed if estMSE = TRUE).

Author(s)

Alicja Wolny-Dominiak, Tomasz Zadlo

References

1. Henderson, C.R. (1950) Estimation of Genetic Parameters (Abstract). Annals of Mathematical Statistics 21, 309-310.
2. Royall, R.M. (1976) The Linear Least Squares Prediction Approach to Two-Stage Sampling. Journal of the American Statistical Association 71, 657-473.
3. Zadlo, T. (2017) On prediction of population and subpopulation characteristics for future periods, Communications in Statistics - Simulation and Computation 461(10), 8086-8104.

Examples

library(lme4)
library(Matrix)


### Prediction of the subpopulation mean based on the cross-sectional data

data(invData) 
# data from one period are considered: 
invData2018 <- invData[invData$year == 2018,] 
attach(invData2018)

N <- nrow(invData2018) # population size
n <- 100 # sample size
# subpopulation of interest: NUTS4type==2
Nd <- sum(NUTS4type == 2) # subpopulation size

set.seed(123456)
sampled_elements <- sample(N,n)
con <- rep(0,N)
con[sampled_elements] <- 1 # elements in the sample
YS <- investments[sampled_elements]
fixed.part <- 'newly_registered'
random.part <- '(1| NUTS2)' 
reg = invData2018[, -which(names(invData2018) == 'investments')]

gamma <- rep(0,N)
gamma[NUTS4type == 2] <- 1/Nd

weights <- rep(1,N) # homoscedastic random components
estMSE <- TRUE 

# Predicted value of the mean in the following subpopulation: NUTS4type==2
EBLUP(YS, fixed.part, random.part, reg, con, gamma, weights, estMSE)$thetaP

# All results
EBLUP(YS, fixed.part, random.part, reg, con, gamma, weights, estMSE)

detach(invData2018)

##########################################################

### Prediction of the subpopulation total based on the longitudinal data

data(invData)
attach(invData)

N <- nrow(invData[(year == 2013),]) # population size in the first period
n <- 38 # sample size in the first period
# subpopulation and time period of interest: NUTS2 == '02' & year == 2018
# subpopulation size in the period of interest:
Ndt <- sum(NUTS2 == '02' & year == 2018)

set.seed(123456)
sampled_elements_in_2013 <- sample(N,n)
con2013 <- rep(0,N)
con2013[sampled_elements_in_2013] <- 1 # elements in the sample in 2013

# balanced panel sample - the same elements in all 6 periods:
con <- rep(con2013,6)

YS <- investments[con == 1]
fixed.part <- 'newly_registered'
random.part <- '(newly_registered  | NUTS4)' 
reg <- invData[, -which(names(invData) == 'investments')]

gamma <- rep(0,nrow(invData))
gamma[NUTS2 == '02' & year == 2018] <- 1

weights <- rep(1,nrow(invData)) # homoscedastic random components
estMSE <- TRUE 

# Predicted value of the total 
# in the following subpopulation: NUTS4type == 2 
# in the following time period: year == 2018
EBLUP(YS, fixed.part, random.part, reg, con, gamma, weights, estMSE)$thetaP

# All results
EBLUP(YS, fixed.part, random.part, reg, con, gamma, weights, estMSE)

detach(invData)

qape documentation built on Aug. 21, 2023, 5:07 p.m.