1
$\begingroup$

I am trying to analyse the influence of temperature on flowering dates of certain plants using survival analysis. I have data from several measuring stations that measure daily mean temperature and give me the date on which the plant of interest has flowered that year.

I figured out how to apply the cox proportional hazard model to my data when looking at a single year, but I am not sure how to extend my analysis to the whole dataset (all years). Here is some sample code of what I achieved till now:

Generate example data (I can also supply the real dataset on request):

library(survival)

events <- data.frame(
    year = rep(1980:1989,5),
    stid = rep(c(rep('a', 10), rep('b', 10), rep('c', 10), rep('d', 10), rep('e', 10))),
    doy = sample(110:170, 50)
)

temp <- data.frame(
  value = round(rep(-(abs(1:365 - 180))^(1/1.5) + 30,50) + runif(365 * 50 , -5, 5),2),
  stid = rep(rep(c(rep('a', 365), rep('b', 365), rep('c', 365),  rep('d', 365),  rep('e', 365))),10),
  day_of_year = rep(1:365,50), 
  year = sort(rep(1980:1989,365*5))
)

Normal coxph analysis for a single year (tmerge requires the newest version of the survival package):

temp = temp[temp$year == 1980,]
events = events[events$year == 1980,]

  # Define the time-ranges for the counting dataset
  dat = tmerge(events, events, id = stid, tstart = doy*0+1, tstop = doy) 

  # Ad a time dependent covariate (tdc)
  dat = tmerge(dat, temp, id = stid, temperature = tdc(temp$day_of_year, temp$value))

  # Ad an event
  dat = tmerge(dat, events, id = stid, flowering = event(events$doy))


surv = Surv(dat$tstart, dat$tstop, dat$flowering)
surv_cox = coxph(surv ~ dat$temp)
summary(surv_cox)

plot(cox.zph(surv_cox))

On Request: Plot of an excerpt of the dataset:

Plot of mean daily temperature with flowering date of that year marked as red line

The actual dataset is much bigger. It goes over 40 years and dozens of measuring stations (plot shows only data from one station).

I want to investigate how the daily temperature influences the flowering date = position of the vertical red lines. I figured out this should be possible with a cox proportional hazard model, but I am not sure how that works for periodically recurrent events

$\endgroup$

migrated from stackoverflow.com Jul 14 '15 at 2:59

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ If you have questions about the appropriateness of a particular statistical analysis, you should post to Cross Validated instead where such matters are on topic. This is not a specific programming question and as such is not a good fit for Stack Overflow. $\endgroup$ – MrFlick Jul 9 '15 at 19:23
  • $\begingroup$ Generally Therneau advises not keeping vectors separate as arguments to Surv but rather this: surv_cox = coxph( Surv(tstart, tstop, flowering) ~ temp, data=dat). This keeps the environments straight. I'm not able to comment on appropriateness because I didn't understand the data setup and the size of the example is too large to fit on one page. Maybe if you plotted the time course of values and did a bit more explaining my feeble brain could figure it out. $\endgroup$ – DWin Jul 11 '15 at 5:39
  • $\begingroup$ I added a plot and some explanation, does this clarify what i want to do? Also, thanks for the tip with the data parameter. $\endgroup$ – Stefan F Jul 13 '15 at 13:21
1
$\begingroup$

Looks like you have a predictor that varies cyclically. My impression from very limited reading of the botanical literature is that something along the lines of cumsum(degree_days) (where the sums are calculated from a point in the middle of winter and a "degree-day" is the number of degrees above a threshold temperature) is used to predict budding and perhaps flowering. So if that were the case then you would first process the temperatures so they were cumulative in whatever sense was appropriate to your specialized domain and then use that as the time-dependent variable. The problem with the multiple years should be handles as a strata(year) entered using the formula. See ?coxph for examples of each of those. Perhaps:

dat$cumtemp <- with(dat, ave( (temp > crit_temp)*(temp-crit_temp), year, cumsum) )

surv_cox = coxph( Surv( tstart, tstop, flowering) ~ cumtemp + strata(year))
summary(surv_cox)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.