Note
I've edited the example to be more intuitive and closer to my real data
Intro
I've got data on customers purchases and with it am trying to predict which customers are more likely to make next purchase at some time in the future. Data consist of customers' features like sex, age etc., and their prior purchase behavior like total spendings and number of orders, one row for every customer. The last two columns are the indicator of wether he have made next purchase or not, and number of days till purchase or till today, in case of no purchase.
Problem
I am building a Cox regression and then want to predict probability of next purchase for individual observations in, say, 30 days from last purchase.
Reproducible example:
library(survival)
library(rms)
library(pec)
library(ggplot2)
data(cost)
# split into train and test sets
set.seed(1)
ind <- sample(1:nrow(cost), 100)
test.set <- cost[ind, ]
train.set <- cost[-ind, ]
For Cox regression I use cph from rms package, for prediction - predictSurvProb from pec package as suggested in this discussion.
# fit Cox model
fit <- cph(Surv(time, status) ~ ., data = train.set, surv = TRUE)
# predict pobability of event in 30 days
test.set$predicted.probs <- 1 - predictSurvProb(fit, newdata = test.set, times = 1000)[, 1]
Thus, for every customer we have his probability of making a purchase in 1000 units of time. I want to validate prediction against real data.
Now to the question: what is the best/valid way to do it?
Here's what I've tried:
I expect that valid model would predict higher probabilities for customers who made their purchase earlier so correlation between probabilities and number of days to event' would be negative and strong (e.g. for customer who actualy made next purchase in 2 days, probability of buying in 30 days would be very high).
with(test.set, cor(predicted.probs, time))
# [1] -0.5221604
Also, probability for those who made purchase (status = 1) would be higher than for those who didn't.
with(test.set, by(predicted.probs, status, mean))
# status: 0
# [1] 0.2371247
# --------------
# status: 1
# [1] 0.4083586
And a graph to eyeball my assumptions:
qplot(data = test.set, x = time, y = predicted.probs, color = time)
Am I correct in my reasoning?
