[R-sig-ME] glmmTMB: model forumlation and prediction questions
Daniel B
dampfschlaghammer at gmail.com
Fri Mar 23 15:49:05 CET 2018
Hi,
I am trying to estimate a mixed effects model based on the assumption that
a dependent variable y is determined by a random effect x. The random
effects x exists for "groupcount" different groups (with cross correlation
between groups)
and across "n" time periods (with autocorrelation to its own lagged
values). I have included some R code, maybe this also helps in
understanding the model.
The formulation I found for glmmTMB is:
model_with_ar <- glmmTMB(y ~ ar1(timestep+0|group), data=outcome_table)
For comparison, I have included an alternative formulation that
disregards the AR1.
I have three important questions:
· Does the formulation of the model “model_with_ar” correspond to
the problem I described in the text?
· How can I do prediction that takes into account the AR1 nature of
the process and therefore?
This means the prediction for the next, unobservd time period should not
assume that all random effects are 0, but instead take into account the
latest realization of the random effects process (given by the BLUP for the
latest time period), the AR1 parameter and the estimated process variance
and covariance. I can of course just manually extract this information, but
is there a more elegant way to do that (especially when I want to apply
this dummy model to real data with many fixed and random effects)?
· How can I correctly identify the cross correlation between the
groups at a certain time period?
The identification of the auto correlation is easy as it is reported in the
model output. Should I do the estimation for the cross correlaton directly
based on the BLUPs for the random effects (as done in the R code)?
Thank you very much!
Regards
require(glmmTMB);require(lme4);require(tsDyn);require(data.table);require(ggplot2)
#Initialize parameters
n <- 1000 #time periods
groupcount <- 10 #number of groups
members_in_group <- 20 #members per group
epsilon_std <- 2 #Random Noise Standard Deviation
#Cross Sectional Variance Covariance Matrix
covariance <- 0.7 #Cross Sectional covariance
variance_per_group <- runif(1,0,1) #Different variance of random effect per
group
varcovmatrix <-
diag(groupcount)+covariance-covariance*diag(groupcount)+diag(groupcount)*variance_per_group
#AR1 Process Matrix (no off-diagonal elements, i.e. no correlation with
lags of other random effects)
arparameter <- 0.9 #Ar Parameter, identical per group
ar_corrmatrix <- diag(groupcount)*arparameter
#Draw realization of random effects over time and groups
re_table <-
data.table(group=as.factor(rep(1:groupcount,each=n)),timestep=as.factor(1:n),
x=as.vector(VAR.sim(B=ar_corrmatrix,n=n,include="none",varcov=varcovmatrix)))
#Draw members in group (maximum of) and select randomly 80% of portfolio
(in order to get unbalanced portfolios)
for (i in 1:members_in_group){
if (i==1) outcome_table <- re_table[sample(1:.N,0.8*.N)]
if (i>1) outcome_table <-
rbind(outcome_table,re_table[sample(1:.N,0.8*.N)])
}
#Add random noise
outcome_table[,y:=x+rnorm(.N)*epsilon_std]
#Old model without ar 1 random effects specification
model_old <- glmmTMB(y ~ (-1+group|timestep),
data=outcome_table,control=glmmTMBControl(optCtrl =
list(iter.max=1e4,eval.max=1e4)))
#New model with ar 1 random effects specification
model_with_ar <- glmmTMB(y ~ ar1(timestep+0|group), data=outcome_table)
#Extract random effects for both models
re_model_with_ar <- t(ranef(model_with_ar)$cond$group)
re_model_old <- ranef(model_old)$cond$timestep
#Derive "empiric" correlation matrix based on simulation
simulated_re_values <-
cor(dcast(re_table,timestep~group,value.var="x")[,!"timestep"])
#Differnce between correlation matrix based on extracted random effects and
"empiric" correlation matrix
delta_with_ar <- cor(re_model_with_ar)-simulated_re_values
delta_old <- cor(re_model_old)-simulated_re_values
sum(delta_with_ar);sum(delta_old)
sum(abs(delta_with_ar));sum(abs(delta_old))
sum(delta_with_ar^2);sum(delta_old^2)
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list