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.