我有以下模型:
我想学习如何实现这个模型。
研究与尝试
我以以下 Poisson GLM 为例:
```{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);
}
```
我了解此示例使用了 QR 重新参数化(请参阅 Stan 参考手册的第 9.2 节)。但是,由于我是 Stan 的新手,我觉得这很吓人,而且我不确定如何以同样的方式实现指数 GLM。甚至需要对指数情况使用 QR 重新参数化吗?
无论如何,我对指数 GLM 的天真尝试是对 Poisson GLM 代码的以下改编:
```{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);
}
```
但是,我完全不确定这是否是在 Stan 中编写这样一个模型的正确方法。我所做的只是尝试使 Poisson GLM 示例适应上述指数 GLM。
我正在寻求指数 GLM 的演示,并澄清我上面的误解。
样本数据
conc time
1 6.1 0.8
2 4.2 3.5
3 0.5 12.4
4 8.8 1.1
5 1.5 8.9
6 9.2 2.4
7 8.5 0.1
8 8.7 0.4
9 6.7 3.5
10 6.5 8.3
11 6.3 2.6
12 6.7 1.5
13 0.2 16.6
14 8.7 0.1
15 7.5 1.3