# Empirical Best Linear Unbiased Predictor 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. ```r 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. ## Returns 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**). ## 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. ## Author(s) Alicja Wolny-Dominiak, Tomasz Zadlo ## Examples ```r 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) ```

EBLUP function