[R] Linear relative rate / excess relative risk models
Wollschlaeger, Daniel
wollschlaeger at uni-mainz.de
Tue Jul 29 13:10:29 CEST 2014
A while ago, I inquired about fitting excess relative risk models in R. This is a follow-up about what I ended up doing in case the question pops up again.
While I was not successful in using standard tools, switching to Bayesian modeling using rstan (mc-stan.org/rstan.html) worked better. The results closely match those from Epicure.
Using the data here: http://dwoll.de/err/dat.txt
The stan model fit below replicates the results from Epicure here: http://dwoll.de/err/epicure.log
Of course I am still interested in learning about other options or approaches.
Daniel
##-------
## rstan code for fitting an excess relative risk model with linear dose-response
## events = personYears * exp(beta0 + beta1*age + beta2*sex) * (1 +beta3*dose)
dat_code <- '
data {
int<lower=0> N;
vector[N] pyears;
vector[N] age;
vector[N] sex;
vector[N] dose;
int<lower=0> event[N];
}
transformed data {
vector[N] log_pyears;
log_pyears <- log(pyears);
}
parameters {
vector[2] beta;
real beta0; # baseline intercept param
real<lower=0> betaD; # dose param -> non-negative
}
model {
# beta0 unspecified -> uniform prior (improper)
beta ~ normal(0, 100); # flat normal prior for params
betaD ~ cauchy(0, 2.5); # ok even if not truncated, cf. stan reference
event ~ poisson_log(log_pyears + beta0 + beta[1]*age + beta[2]*sex +
log(1 + betaD*dose));
}
'
library(rstan)
stan_dat <- with(dat,
list(pyears=pyears,
age=age,
sex=sex,
N=length(pyears),
event=event,
dose=dose))
stanFit <- stan(model_code=dat_code, data=stan_dat,
thin=5, iter=10000, chains=2)
traceplot(stanFit)
stanFit
-----Original Message-----
From: Wollschlaeger, Daniel
Sent: Thursday, January 9, 2014 10:44 AM
To: David Winsemius
Cc: r-help at r-project.org
Subject: RE: AW: [R] Linear relative rate / excess relative risk models
Thanks for your suggestions! Here are links to simulated data and the Epicure syntax + reference fit:
http://dwoll.de/err/dat.txt
http://dwoll.de/err/epicure.log
The model tested in Epicure is
lambda = exp(alpha0 + alpha1*agePyr)*(1 + beta0*dosePyr*exp(beta1*agePyr))
with counts in variable event and offset pyears.
Many thanks, D
> -----Original Message-----
> From: David Winsemius [mailto:dwinsemius at comcast.net]
> Sent: Thursday, January 09, 2014 4:33 AM
> To: Wollschlaeger, Daniel
> Cc: r-help at r-project.org
> Subject: Re: AW: [R] Linear relative rate / excess relative risk models
>
>
> On Jan 8, 2014, at 3:22 PM, Wollschlaeger, Daniel wrote:
>
> > If I understand you correctly, that is exactly the approach taken by
> > Atkinson & Therneau: They get the baseline rates from published rate
> > tables from the general population, multiply them by the appropriate
> > person-time from their data to get expected counts, and use this as
> > offset.
> >
> > Unfortunately, I won't have comparable baseline rate tables. And
> > while I could fit a separate model only to the unexposed group for
> > expected counts, I'd prefer to fit both factors (lambda0 and 1+ERR)
> > simultaneously - as it is typically done in the existing literature.
>
> If you would describe your data situation more completely (ideally
> with a reproducible example) you might get a better answer. It's also
> considered polite on this mailing list to include the email chain, so
> appending original question:
>
> --
> David
>
> >
> > Best, Daniel
> >
> > ________________________________________
> > Von: David Winsemius [dwinsemius at comcast.net]
> > Gesendet: Mittwoch, 8. Januar 2014 19:06
> > An: Wollschlaeger, Daniel
> > Cc: r-help at r-project.org
> > Betreff: Re: [R] Linear relative rate / excess relative risk models
> >
> > I would fit a Poisson model to the dose-response data with offsets
> > for the baseline expecteds.
>
> David Winsemius, MD
> Alameda, CA, USA
>
> ============================
> My question is how I can fit linear relative rate models (= excess
> relative risk models, ERR) using R. In radiation epidemiology, ERR
> models are used to analyze dose-response relationships for event rate
> data and have the following form [1]:
>
> lambda = lambda0(z, alpha) * (1 + ERR(x, beta))
>
> * lambda is the event rate
> * lambda0 is the baseline rate function for non-exposed persons and
> depends on covariates z with parameters alpha
> * ERR is the excess relative risk function for exposed persons and
> depends on covariates x (among them dose) with parameters beta
> * lambda/lambda0 = 1 + ERR is the relative rate function
>
> Often, the covariates z are a subset of the covariates x (like sex and
> age). lambda is assumed to be log-linear in lambda0, and ERR typically
> has a linear (or lin-quadratic) dose term as well as a log-linear
> modifying term with other covariates:
>
> lambda0 = exp(alpha0 + alpha1*z1 + alpha2*z2 + ...)
> ERR = beta0*dose * exp(beta1*x1 + beta2*x2 + ...)
>
> The data is often grouped in form of life tables with the observed
> event counts and person-years (pyr) for each cell that results from
> categorizing and cross-classifying the covariates. The counts are
> assumed to have a Poisson-distribution with mean mu = lambda*pyr, and
> the usual Poisson-likelihood is used. The interest is less in lambda0,
> but in inference on the dose coefficient beta0 and on the modifier
> coefficients beta.
>
> In the literature, the specialized Epicure program is almost
> exclusively used. Last year, a similar question on R-sig-Epi [2] did
> not lead to a successful solution (I contacted the author). Atkinson &
> Therneau in [3] discuss excess risk models but get lambda0 separately
> from external data instead of fitting lambda0 as a log-linear term.
> Some R packages sound promising to me (eg., gnm, timereg) but I
> currently don't see how to correctly specify the model.
>
> Any help on how to approach ERR models in R is highly appreciated!
> With many thanks and best regards
More information about the R-help
mailing list