Visualizing a two way interaction for a binary dependent variable I have a dataset that has two continuous independent variables, months and tenure. The target variable is binary, 1/0. I want to visualize how the events change across months*tenure. How can I do that? I tried proc g3d in SAS:
proc g3d data=ds;
plot months*tenure=target;
run; 

But the results are very large because there are about million data points. What is the best way to visualize such data interaction?
 A: Interestingly enough, we don't already seem to have a question on visualizing interactions with binary dependent data.
The most interesting aspect of your question is that you have binary data. And a lot of it, so I'd recommend plotting the fits from a suitable model, for instance a logistic regression. Let's generate some toy data and fit a model (I'll use R, because I know it better):
set.seed(1)
nn <- 1000
months <- sample(x=0:60,size=nn,replace=TRUE)
tenure <- sample(x=0:10,size=nn,replace=TRUE)
target <- 1/(1+exp(-0.02*months-0.1*tenure-0.01*months*tenure+rnorm(nn,sd=4)))>.95
summary(target)

model <- glm(target~months*tenure,family="binomial")
summary(model)

Alright, now we need to decide over which range of months and tenure we wish to plot the model fit.    
mm <- 0:60
tt <- 0:10

I see two ways of plotting this. Probably best to try both and see which one works best.
The first way would be to plot a 3d surface:
foo <- expand.grid(mm,tt)
colnames(foo) <- c("months","tenure")

fit <- cbind(foo,predict(model,newdata=foo,type="response"))
colnames(fit) <- c("months","tenure","target")

fit.matrix <- outer(mm,tt,
  function(months,tenure)predict(model,
    newdata=data.frame(months=months,tenure=tenure),type="response"))

opar <- par(mai=c(0.5,0.5,0,0))
  persp(x=mm,y=tt,z=fit.matrix,
    xlab="Months",ylab="Tenure",zlab="Fitted target",
    ticktype="detailed",theta=-30)
par(opar)


The second possibility would be to plot the response to one variable, say months, for different values of the other variable. I like setting the other variable to the quartiles in the originally observed data:    
probs <- c(.25,.5,.75)
plot(range(mm),c(0,1),type="n",xlab="Months",ylab="Fitted probability of target")
for ( ii in seq_along(probs) ) {
  lines(mm,
    predict(model,newdata=data.frame(months=mm,tenure=quantile(tenure,probs[ii])),type="response"),
    lty=ii)
}
legend("bottomright",lty=seq_along(probs),
  legend=paste("Tenure:",quantile(tenure,probs)))


Or of course, the other way around:
plot(range(tt),c(0,1),type="n",xlab="Tenure",ylab="Fitted probability of target")
for ( ii in seq_along(probs) ) {
  lines(tt,
    predict(model,newdata=data.frame(months=quantile(months,probs[ii]),tenure=tt),type="response"),
    lty=ii)
}
legend("bottomright",lty=seq_along(probs),
  legend=paste("Months:",quantile(months,probs)))


A: Ok. here is what worked for me. I grouped the targets/total count at each month*tenure combination and was able to plot it.
