library(lme4)
library(Matrix)
library(mvtnorm)
library(matrixcalc)
library(future.apply)
data(invData)
# data from one period are considered:
invData2018 <- invData[invData$year == 2018,]
attach(invData2018)
N <- nrow(invData2018) # population size
con <- rep(1,N)
con[c(379,380)] <- 0 # last two population elements are not observed
YS <- log(investments[con == 1]) # log-transformed values
backTrans <- function(x) exp(x) # back-transformation of the variable of interest
fixed.part <- 'log(newly_registered)'
random.part <- '(1|NUTS2)'
reg <- invData2018[, -which(names(invData2018) == 'investments')]
weights <- rep(1,N) # homoscedastic random components
# Characteristics to be predicted:
# values of the variable for last two population elements
thetaFun <- function(x) {x[c(379,380)]}
set.seed(123)
predictor <- plugInLMM(YS, fixed.part, random.part, reg, con, weights, backTrans, thetaFun)
predictor$thetaP
### Estimation of prediction accuracy
est_accuracy <- bootParFuture(predictor, 10, c(0.75,0.9))
# Estimation of prediction RMSE
est_accuracy$estRMSE
# Estimation of prediction QAPE
est_accuracy$estQAPE
# [,1] [,2]
# 75% 1370.823 180.0514
# 90% 1477.444 249.7517
####### Interpretations in case of prediction of investments
####### for population element no. 379:
### It is estimated that at least 75% of absolute prediction errors are
# smaller or equal 1370.823 milion Polish zloty
# and at least 25% of absolute prediction errors are
# greater or equal 1370.823 milion Polish zloty.
### It is estimated that at least 90% of absolute prediction errors are
# smaller or equal 1477.444 milion Polish zloty
# and at least 10% of absolute prediction errors are
# greater or equal 1477.444 milion Polish zloty.
detach(invData2018)
Run the code above in your browser using DataLab