Autoregressive effects of residuals

[paul-buerkner]

Autoregressive (AR) effects as they are currently implemented refer to effects of the response on itself. Some other packages (such as nlme) use the AR term to refer to autoregression of the residuals.

Thus, it would be nice to have AR effects of the residuals in addition to AR effects of the response.

[paul-buerkner]

AR effects of residuals are now implemented. They are replacing the previous AR effects to be in accordance with the nlme implementation of autocorrelation structures.

AR effects of the response can still be modeled but are now called ARR effects to avoid using the term 'AR' for two different things.

I will leave this issue open until I am confident that the implementation is stable.

[BERENZ]

Could you provide an example how to specify such model? I would like to test the implementation. :)

edit:
I found cor_arr however, there is no information about that in the description of cor_brms

[paul-buerkner]

This is an example model using a data set that comes with brms:

fit_brm <- brm(count ~ log_Age_c + log_Base4_c * Trt_c + (1|patient),
data = epilepsy, autocor = cor_arma(~visit|patient, 1, 1),
prior = c(set_prior("normal(0,0.5)", class = "ar"),
set_prior("normal(0,0.5)", class = "ma")))

This model defines both, AR1 and MA1 effects (of residuals). If you want AR1 effects only, just replace cor_arma(~visit|patient, 1, 1) with cor_arma(~visit|patient, 1) or simpler cor_ar(~visit|patient).

If you want to test AR effects of the response (previously called AR and now called ARR), you use
cor_arr(~visit|patient) for ARR1. Note that the name 'ARR' was my (probably not ideal) idea, and you will not find it in the literature.

The cor_brms is just the overall class of correlation structures and you should not worry about it. They only class that is currently relevant is cor_arma.

[BERENZ]

Hi,

results that I've got are different from those that supposed to be…(I've reinstalled package to github version 0.5.0.9000). It looks like the AR(1) still refers to the AR effect of the response not for the residuals. Below please find the code for the sae model that has two random effect - domain and time | domain (AR(1)).

library(sae2)
library(sae)
library(brms)
# Loading 'brms' package (version 0.5.0.9000). Useful instructions 
# can be found by typing help('brms'). A more detailed introduction 
# to the package is available through vignette('brms').

data(spacetime)      # Load data set from sae package
spacetime$sei <- spacetime$Var^.5
D <- nrow(spacetimeprox)            # number of domains
T <- length(unique(spacetime$Time)) # number of time instants

# Fit model with AR(1) time effects for each domain (domain is already present)
resultT.RY  <- eblupRY(Y ~ X1 + X2, D, T, vardir = diag(spacetime$Var),
                       data=spacetime, ids=spacetime$Area)

# Fit the same model with brms
fit_brm <- brm(Y | se(sei) ~  X1 + X2 + (1|Area),
               data = spacetime, 
               autocor = cor_ar(~ Time | Area),
               prior = c(set_prior("beta(2,2)", class = "ar"),
                         set_prior("normal(0,5)", class = "b")))
fit_brm_fitted <- fitted(fit_brm)

Results obtained from sae2 package (sig2_u refers to AR(1) effect, sig2_v refers to domain effect, rho is the autocorrelation). The same results I've obtained in metafor package.

> resultT.RY$fit
$model
[1] "T: Rao-Yu, REML"

$convergence
[1] TRUE

$iterations
[1] 40

$estcoef
                 beta std.error    tvalue       pvalue
(Intercept)  2.324605 0.4488063  5.179527 2.224491e-07
X1          -2.606220 0.8264804 -3.153396 1.613827e-03
X2          -1.830352 0.1989480 -9.200151 3.574481e-20

$estvarcomp
           estimate    std.error
sig2_u 0.0008172101 0.0004979779
sig2_v 0.0000000033 0.0009056348
rho    0.3331617406 0.7526800594

$goodness
          loglike restrictedloglike 
         62.57172          55.87707 

and here are results from brms package

> fit_brm
 Family: gaussian (identity) 
Formula: Y | se(sei) ~ X1 + X2 + (1 | Area) 
   Data: spacetime (Number of observations: 33) 
Samples: 2 chains, each with n.iter = 2000; n.warmup = 500; n.thin = 1; 
         total post-warmup samples = 3000
   WAIC: -115.39

Random Effects: 
~Area (Number of levels: 11) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.08      0.02     0.04     0.14        541 1.01

Correlation Structure: arma(~Time|Area, 1, 0, 0)
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
ar[1]     0.19      0.12     0.03     0.47       1194    1

Fixed Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     1.47      0.45     0.57     2.33        997    1
X1           -2.07      0.82    -3.59    -0.40       1050    1
X2           -0.53      0.21    -0.94    -0.14       1252    1

Samples were drawn using NUTS(diag_e). For each parameter, Eff.Sample is a 
crude measure of effective sample size, and Rhat is the potential scale 
reduction factor on split chains (at convergence, Rhat = 1).

[paul-buerkner]

Thanks for testing! brms is definitly fitting AR effects of residuals now when using cor_ar, so this cannot explain the differences.

I have another idea though. The eblupRY function (what I great name by the way), seems to fit 2 variances (sig_u and sig_v) and 1 correlation (rho), whereas the brms model only has 1 variance / sd (for Area) and 1 correlation (ar[1]). So I think that these two models are not equivalent.
From what you wrote, I think that the brms model should rather look like this

fit_brm <- brm(Y | se(sei) ~  X1 + X2 + (1|Area) + (1|Time:Area),
               data = spacetime, 
               autocor = cor_ar(~ Time | Area),
               prior = c(set_prior("normal(0,0.5)", class = "ar"),
                             set_prior("normal(0,5)", class = "b")))

with an additional random effect (1|Time:Area). I also changed the prior of the AR effect, as beta priors restrict the parameter to [0,1], but autocorrelations may be between [-1,1] or even lie outside this interval if the process is not stationary.

The results look very similar to that of eblupRY:

 Family: gaussian (identity) 
Formula: Y | se(sei) ~ X1 + X2 + (1 | Area) + (1 | Time:Area) 
   Data: spacetime (Number of observations: 33) 
Samples: 2 chains, each with n.iter = 2000; n.warmup = 500; n.thin = 1; 
         total post-warmup samples = 3000
   WAIC: -143.93

Random Effects: 
~Area (Number of levels: 11) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.02      0.02        0     0.08        286 1.01

~Time:Area (Number of levels: 33) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.03      0.01     0.01     0.04        551    1

Correlation Structure: arma(~Time|Area, 1, 0, 0)
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
ar[1]     0.21       0.4    -0.54     0.97        782    1

Fixed Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     2.18      0.51     1.15     3.15       1059 1.00
X1           -2.50      0.94    -4.36    -0.64       1163 1.00
X2           -1.63      0.35    -2.17    -0.75        365 1.01

Samples were drawn using NUTS(diag_e). For each parameter, Eff.Sample is a 
crude measure of effective sample size, and Rhat is the potential scale 
reduction factor on split chains (at convergence, Rhat = 1).

When using the autocor argument in brm, only the autocorrelations itself is added to the model (in this case ar[1]), and one has to specify the random effects seperately. So when I refer to the AR[1] effect, I always mean the autocorrelation not the residual variance (which seems to be sig2_u in the sae2 model). For the present brms model, the residual variance is modeled through the random effect of Time:Area.

[BERENZ]

Great! Thanks for clarification, I should noticed that my results from brm function don't have second variance component.

However, the I am not sure whether specifying (1|Time:Area) and cor_ar(~ Time | Area) separately would give the same results. In that model the second component is assumed to follow AR(1) process for residuals given by: u_{t} = \rho u_{t-1} + e where e \sim N(0,\sigma^2). Thus, the variance component and rho are estimated together.

The model Y | se(sei) ~ X1 + X2 + (1 | Area) + (1 | Time:Area) would result in:

1. random effect for Area
2. random effect in Time/Area
3. autocorrelation of residuals independent from (1) and (2).

Nice example of different variance structure provide proc mixed documentation. Toepliz matrix is nice for specifying AR effects, that's pitty that stan do not have that...

In addition, thanks for commend on prior. I've used beta(2,2) because I know that in this case \rho would be between (0,1). For other models I would use normal prior.

[paul-buerkner]

The autocorrelation is not independent from the random effects components, even though they are defined in seperate parts of the model specification. This is because fixed effects, random effects, and the autocorrelation components all go into the same linear regression for Y and thus will all be estimated together.

Since the random effect of Time:Area has the same number of levels as observations in the data, it corresponds to the residuals e, and leads to the equation u_{t} = \rho u_{t-1} + e with e \sim N(0,\sigma^2), when combined with an AR(1) autocorrelation (even if they are not explicitly defined together).

Thanks for the link! I will have a look into this. :-)

[BERENZ]

I've looked into the code that brm provide for the above model, in particular transformed parameters

transformed parameters { 
  vector[N] eta;  # linear predictor 
  row_vector[Karma] E[N];  # ARMA design matrix 
  vector[N] e;  # residuals 
  vector[N_1] r_1;  # REs 
  vector[N_2] r_2;  # REs 
  # compute linear predictor 
  eta <- X * b + b_Intercept; 
  E <- E_pre; 
  r_1 <- sd_1 * (pre_1);  # scale REs 
  r_2 <- sd_2 * (pre_2);  # scale REs 
  # if available add REs to linear predictor 
  for (n in 1:N) { 
    eta[n] <- eta[n] + Z_1[n] * r_1[J_1[n]] + Z_2[n] * r_2[J_2[n]]; 
    # calculation of ARMA effects 
    e[n] <- (Y[n]) - eta[n]; 
    for (i in 1:Karma) { 
      if (n + 1 - i > 0 && n < N && tgroup[n + 1] == tgroup[n + 1 - i]) { 
        E[n + 1, i] <- e[n + 1 - i]; 
      } 
    } 
    eta[n] <- (eta[n] + head(E[n], Kar) * ar); 
  } 
} 

and as I see here Z_1[n] * r_1[J_1[n]] refers to the first random effect, Z_2[n] * r_2[J_2[n]] refers to the second random effect and the AR is calculated for the residual given by e[n] <- (Y[n]) - eta[n]. That means the model calculate rho after simulating the random effects and sigma2_2 (for the second component) is not connected with the rho itself. While, AR(1) effect should have variance structure as in the proc mixed link.

I have found some example at Stan google group https://groups.google.com/forum/#!searchin/stan-users/AR/stan-users/Cq85_k0QcRo/3mEmkO3p5pgJ

data {
int<lower=0> N; // length of data set
vector[N] y;    // data
}

        parameters {
            real mu; 
            real<lower=0> sigma;
            real<lower=0,upper = 1> rho;
        }

        transformed parameters {
        // declare mu, Sigma, chol_Sigma
            vector[N] Mu;
            matrix[N,N] Sigma;
            matrix[N,N] chol_Sigma;

        // define mu
            for(i in 1:N)
                Mu[i] <- mu;

        // define Sigma
            for(i in 1:N)
            for(j in 1:N)
                Sigma[i,j] <- pow(sigma,2)*pow(rho,abs(i-j));

        // define chol_Sigma
            chol_Sigma <- cholesky_decompose(Sigma);
        }

        model {
         // prior distribution
                rho ~ uniform(0.0,1.0); 
         // likelihood
            y ~ multi_normal_cholesky(Mu, chol_Sigma);
        }

EDIT: Maybe there is misunderstanding between us in the sense what I refer as AR process in residuals. In mixed models the difference between observed values and fixed effects are residuals for which we assume various covariance structure. Then, after estimating the model with fixed and random effect we have conditional residuals which are free of the covariance structure and follow N(0,\sigma) where \sigma in general is unknown, in small area is sampling variance. Then, the way how the AR process is implemented now is not the same as is assumed in mixed model theory.

EDIT2: In November I will have finally enough time to work on that in Stan. Probably the best way is to code as on the Stan google group.

[paul-buerkner]

You are right. My implementation becomes problematic as soon as the user defines known SEs.

If this is not the case, the computation order should be correct (at least the models can accurately recover parameters from simulations).

I will have a deeper look at this in the next days and see how I can implement AR structures for known SEs. Thanks for taking the time to test this! :-)

[BERENZ]

I've adapted the code from Stan google group for SAE model. Hope it will be useful for you. :)
diag_matrix(rep_vector(sd_1,N)) could be used to specify user-defined matrices for covariance structure.

library(sae)
library(brms)

data_in <- list(N=nrow(milk),
                 Y = milk$yi,
                K = 4,
                X = model.matrix(~as.factor(MajorArea),milk),
                samperr = milk$SD^2)
model <- 
'
data { 
  int<lower=1> N;  
  vector[N] Y;  
  vector[N] samperr; 
  int<lower=1> K;  
  matrix[N, K] X;  
} 

parameters { 
  vector[K] b;  # fixed effects 
  real<lower=0> sd_1;  # RE standard deviation 
} 
transformed parameters { 
  vector[N] eta;
  matrix[N,N] Sigma;
  matrix[N,N] chol_Sigma;
  eta <- X * b; 
  Sigma <- diag_matrix(rep_vector(sd_1,N)) + diag_matrix(samperr);
  chol_Sigma <- cholesky_decompose(Sigma);

} 
model { 
  # prior specifications 
  sd_1 ~ cauchy(0,5); 
  # likelihood contribution 
  Y ~ multi_normal_cholesky(eta, chol_Sigma); 
} 
'

fit <- stan(model_code = model, 
            data= data_in)

[paul-buerkner]

I just pushed some commits to github allowing to fit AR1 processes in models with predefined SEs using AR1 covariance matrices. For the spacetime example, the following code should work.

fit_brm <- brm(Y | se(sei) ~ X1 + X2 + (1|Area),
               data = spacetime,
               autocor = cor_ar(~ Time | Area),
               prior = c(set_prior("normal(0,1)", class = "ar"),
                         set_prior("normal(0,5)", class = "b")))

For the time being, this implementation of AR1 effects will only be used in models with predefined SEs and I will run more tests in the next few days to make sure that the models work as intended.

[BERENZ]

Thanks for info! I will test it during this week. :)

[paul-buerkner]

MA(1) and ARMA(1,1) effects using covariance model formulation are now also implemented.

In addition, I added the logical argument cov to the cor_arma / cor_ar / cor_ma functions, indicating if covariance model formulation should be used. The code for the spacetime example with an AR(1) effect now looks like this:

fit_brm <- brm(Y | se(sei) ~ X1 + X2 + (1|Area),
               data = spacetime,
               autocor = cor_ar(~ Time | Area, cov = TRUE),
               prior = c(set_prior("normal(0,1)", class = "ar"),
                         set_prior("normal(0,5)", class = "b")))