[R-sig-ME] glmmTMB: model forumlation and prediction questions

[Daniel B]

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]]