Implementing Exponential General Linear Model


#1

I have the following model:

time_i | \lambda_i \sim Exp(\lambda_i), i = 1, \dots, n

log(\lambda_i) = - \beta_0 - \beta_1 conc_i, i = 1, \dots, n,

where Exp(\lambda) is the exponential distribution with rate parameter \lambda.

Also, linking the rate parameter \lambda to the covariate via log(λ) = −\beta_0 − \beta_1 x is equivalent to linking the mean of the response to the covariate via log(\mu) = \beta_0 + \beta_1 x.

I would appreciate it if someone could please demonstrate how to implement this model.

I’ve had a look at the following Poisson GLM as an example:

Y_i | \mu_i = Pois(\mu_i), i = 1, \dots, n

log(\mu_i) = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}, i = 1, \dots, n

```{stan output.var="PoissonGLMQR"}
  data{
    int<lower=1> n;     // number of observations
    int<lower=1> p;     // number of predictors
    matrix[n, p] x;     // design matrix
    int<lower=0> y[n];  // responses
  }
  
  transformed data{
    matrix [n, p] Q_ast;
    matrix [p, p] R_ast;
    matrix [p, p] R_ast_inverse;
    Q_ast = qr_Q(x)[,1:p] * sqrt(n - 1.0);
    R_ast = qr_R(x)[1:p,] / sqrt(n - 1.0);
    R_ast_inverse = inverse(R_ast);
  }
  
  parameters{
    vector[p] theta;  // regression coefficients for predictors
  }
  
  transformed parameters{
    vector[p] beta;
    vector[n] mu;
    mu = exp(Q_ast*theta);
    beta = R_ast_inverse * theta;
  }
  
  model{
     y  ~ poisson(mu);
  }
```

I understand that this example has used QR reparameterization (see section 9.2 of the Stan reference manual). However, since I’m new to Stan, I’m finding this quite intimidating, and I’m unsure of how to implement the exponential GLM in the same way. Does one even need to use the QR reparameterization for the exponential case?

Anyway, my naive attempt at the exponential GLM is the following adaptation of the poisson GLM code:

```{stan output.var="ExponentialGLM"}
  data{
    int<lower=1> n;     // number of observations
    int<lower=1> p;     // number of predictors
    matrix[n, p] x;     // design matrix
    real<lower=0> y[n];  // responses
  }
  
  transformed data{
    matrix [n, p] Q_ast;
    matrix [p, p] R_ast;
    matrix [p, p] R_ast_inverse;
    Q_ast = qr_Q(x)[,1:p] * sqrt(n - 1.0);
    R_ast = qr_R(x)[1:p,] / sqrt(n - 1.0);
    R_ast_inverse = inverse(R_ast);
  }
  
  parameters{
    vector[p] theta;  // regression coefficients for predictors
  }
  
  transformed parameters{
    vector[p] beta;
    vector[n] mu;
    mu = exp(Q_ast*theta);
    beta = (R_ast_inverse * theta);
  }
  
  model{
     y  ~ exponential(mu);
  }
```

In the above code, I attempted to use the fact that linking the rate parameter \lambda to the covariate via log(λ) = −\beta_0 − \beta_1 x is equivalent to linking the mean of the response to the covariate via log(\mu) = \beta_0 + \beta_1 x. This is why I coded mu instead of lambda.

But, as I said, I’m totally unsure if this is even the right way to go about coding such a model in Stan. All I did was simply try to adapt the poisson GLM example to the aforementioned exponential GLM.

I would appreciate it if someone could please demonstrate the exponential GLM and, if it isn’t too much trouble, please take the time to clarify my misunderstands above.


#2

If you are using the exponential distribution for the likelihood, the outcome should be a vector rather than an integer array.


#3

Thanks for the response, Ben; much appreciated. I changed it to real<lower=0> y[n]; and the model now compiles fine.

What about my implementation of the model itself? I’m assuming my attempt here is incorrect.


#4

Improper priors are never correct


#5

Thanks, ben. I think I’ve got it now.

The QR reparameterization version of the model is as follows:

```{stan output.var="ExponentialGLM2"}
  data{
    int<lower=1> n;     // number of observations
    int<lower=1> p;     // number of predictors
    matrix[n, p] x;     // design matrix
    real<lower=0> y[n];  // responses
  }

  transformed data{
    matrix [n, p] Q_ast;
    matrix [p, p] R_ast;
    matrix [p, p] R_ast_inverse;
    Q_ast = qr_Q(x)[,1:p] * sqrt(n - 1.0);
    R_ast = qr_R(x)[1:p,] / sqrt(n - 1.0);
    R_ast_inverse = inverse(R_ast);
  }

  parameters{
    vector[p] theta;  // regression coefficients for predictors
  }

  transformed parameters{
    vector[p] beta;
    vector[n] lambda;
    lambda = exp(-Q_ast*theta);
    beta = (R_ast_inverse * theta);
  }

  model{
     y  ~ exponential(lambda);
     // theta[1] ~ normal();
     // theta[2] ~ normal();
  }
```

As per your suggestion, I’m going to use normal priors. Does everything look fine now?


#6

Syntax-wise, yes.