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@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:dwinsem...@comcast.net] > Sent: Thursday, January 09, 2014 4:33 AM > To: Wollschlaeger, Daniel > Cc: r-help@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 [dwinsem...@comcast.net] > > Gesendet: Mittwoch, 8. Januar 2014 19:06 > > An: Wollschlaeger, Daniel > > Cc: r-help@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 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
