arpeters / survivalexercises Goto Github PK

@thomasnwilson & @ARPeters
I made a silly goof yesterday. Once we discovered that survSplit uses <= logic for its cutoff, I moved the cutoff in the wrong direction (from 250 to 251) when trying to mimic our count process model. I should have come down one survival unit (from 250 to 249). Once you do this, the models are nearly identical.

For next week, let's move on to Chpt 8 in K&K. Please do the practice and test exercises for next week. Thanks!

@thomasnwilson, we joked about it as you left, but in all seriousness, less than a minute after the door shut behind you, @ARPeters and I figured out the connection between the textbook reference to the Table 9.6-9.8 marginal probability calculations and the CIC formula (so, obviously, you were holding us back:). The mistake I was making during class was overlooking the lack of a subscript on the S(t_[f-1]) term in the incidence calculations, Ic(t_f). S(t_[f-1]) reflects the survival at the last time period from any and all causes (including the competing risks). Thus in Table 9.6, S(t_[f-1]) at t_f = 5, is actually 1/100 or .01, and that means Ic(t_f = 5) is .01 * 1/1 = .01 as is the cumulative incidence (since CIC(t_f < 5) = 0). This .01 estimate is the the same as the marginal probability estimate demonstrated in Table 9.8. Thought I'd better write this all down before I forget it. We can discuss more next week.

@wibeasley , @bard1536 , I have created a private repository where we can share code, and I have declared you two to be collaborators, but what is Thomas Wilson's profile name?

@thomasnwilson @ARPeters
I found my mistake from yesterday. I wasn't handling censoring correctly for those who ended participation before the final event time. Below is the key segment of the corrected code. Notice I take the minimum of the observation period and the unique event time into account when setting 'dir'. Thus, if an observation is censored between two unique event times, then 'dir' is set equal to the censoring time (not the next event time) and the last unique event time (e.g., say subject A drops out without an event at time = 15, and the event times that surround this censoring time are 12 and 16; then dir becomes 15 - 12 = 3, and not 16 - 12 = 4). Also, once the censoring time is reached in the sequence of unique event times (i.e., first time eventTimes[i] > SurvivalTime for a censored case), then you need to exit the for loop below (notice I changed the logic for entering the second to last 'if()' command so that censored cases can get inside and then exit the parent for loop).

Make these changes and the counting process (not point process- I was using the wrong terminology yesterday) glm results will look near identical to the coxph results. See if you can now take this a step further and replicate some of the stratified results from Chpts 4 & 5. Good luck!

createPTable <- function(d){
dNew <- d
for(i in 1:length(eventTimes)){
dNew[i,] <- d
dNew[i,"r"] <- i
dNew[i,"tr"] <- eventTimes[i]
dNew[i,"dir"] <- ifelse(i==1,min(d$SurvivalTime,eventTimes[i]),min(d$SurvivalTime,eventTimes[i]) - eventTimes[i-1])
dNew[i,"yir"] <- 0
if(d$SurvivalTime <= eventTimes[i]) {
if(d$Status %in% 1) dNew[i,"yir"] <- 1
break
}
}
return(dNew)
}

Just synced new code that solved that problem we experienced Tuesday. The fix we tried in class was correct, but I didn't realize we were replacing a data.frame name (ptProcessChemo) with itself after the looping code (cautionary tale: not best practice to continually write over an old dataset if you ever submit code line by line like I do). Thus, I was applying the correct code to an already corrupted data.frame. I should have created the data.frame anew and then applied the code. Anyway, it works now, although we still can't replicate K&K's table b/c our sample size is different from theirs.

Finally, I rewrite the Time1-Time3 variable code to make it more efficient. Probably worth browsing when you get a second (I commented out the old method).