There is this open source R package: https://cran.r-project.org/web/packages/morse/vignettes/tutorial.html

which has models for survival data and for reproduction data. This is not exactly the solution to your problem as these are concentration-response models, whereas you seem to imply that you have qualitative treatment factors in your experimental setting.

However this could give you ideas on how to model your data. The package relies on JAGS which is similar to the Stan software mentioned in the previous response.

Bjoern and Guillaume,

Thanks so much for your insights! 

I wanted to add a few clarifications about my study, that relate to your suggestions.

The way I arranged my survival data, the timeline begins at day zero of offspring life. I basically aligned all offspring, whether born in the first or last clutch, on a day zero to event/censor timeline.  

My treatments are entirely discrete, that is, they are not increasing in level or intensity. Guillaume, I think this means that Morse won't work here, right?

I have a total of 73 breeding pairs, who produced a total of ~560 offspring. I threw in a KM curve, by treatment, for perspective. It seems clear (not that visual impressions always count for much) that there are strong differences between at least two of the treatments.

Adding a random effect of pair to a model with treatment as the main effect, raises the log likelihood from -698 to -578, quite a change. But, when I do this, all p-values rise to ~0.5 (and coefficients for treatments now overlap widely with zero). I assume that because pairs vary so extensively in their offspring survival between treatments and within, using a random pair effect serves to mask other effects (hence my idea of a random effect of pairs nested in treatments, which I cannot do in Survreg or Proc Lifereg). Does this make any sense? Any advice for a non-statistician?

Thanks again!

To answer your question, no Morse will not work as is in your case.

Your comment about between and within treatment variation makes sense to me. If I understand your point correctly you could even plot the KM curves for each pairs, use dotted/dashed lines per treatment and visualise the relative between and within treatment variation.

I know of no easy way to model a random effect of pairs nested in treatments except by using a software such as BUGS/JAGS or Stan. However, using these you would not have to write your own likelihood exactly but rather to write the different components of your model as somewhat independent pieces and the software would take care of the inference.

Admittedly, this may remain somewhat involved but you can find tons of help, for instance on https://groups.google.com/forum/#!categories/stan-users/general