Graphing and Analyzing Probit Regression I am currently analyzing a data set having to do with oak tree mortality. I am trying to understand the correlation between dead crowns (dead = 0, live = 1) and the distance to (1) nearest dead crown and (2) distance to nearest bay tree. (See Below).
I have run a probit regression model:
Living<-sample(0:1,10,replace=T)
    Dead_Dist<-c(158.68,99.62,64.42,64.42,86.92,117.77,41.81,41.81,54.73,64.35)
    Bay_Dist<-c(92.47,179.92,317.73,365.23,58.70,193.23,330.36,123.14,88.65,72.34)
    mydata<-cbind(Living,Dead_Dist,Bay_Dist)
    mydata1<-data.frame(mydata)
    myprobit <- glm(Living ~ Bay_Dist + Dead_Dist, family = binomial(link = "probit"),data = mydata1)
    myprobit

But I am at a loss for how to plot and interpret the results, as I am fairly new to R and regressions.
If anyone has suggestions with how to proceed with this analysis, I would appreciate it!
 A: First, here are some fake results to consider using set.seed(8) before sampling Living:
Call:
glm(formula = Living ~ Bay_Dist + Dead_Dist, family = binomial(link = "probit"), 
data = mydata)

Deviance Residuals: 
Min       1Q   Median       3Q      Max  
-1.7826  -1.0012   0.7055   0.8691   1.4256  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.5650207  1.5138702   1.034    0.301
Bay_Dist    -0.0002185  0.0040025  -0.055    0.956
Dead_Dist   -0.0159294  0.0133800  -1.191    0.234

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 13.460  on 9  degrees of freedom
Residual deviance: 11.793  on 7  degrees of freedom
AIC: 17.793

In this sample it would seem the effects are pretty weak (small estimates, similar in size to their standard errors), but we can't say with any certainty that the effects really are weak without more data. Nothing is crossing any conventional significance thresholds (for whatever those are worth), but this is to be expected with such a small sample. Distance to the nearest dead crown is doing more for the model than distance to the nearest bay tree though (larger z).
One thing you could do is check product terms (squared/cubed predictors and interactions), but this will cost precious degrees of freedom in a small sample. As for plotting, I'd recommend ggplot2. I'm somewhat new to that package myself, so I'll offer a very simple plot:
require(ggplot2);ggplot(mydata,aes(x=Dead_Dist,y=Living))+
geom_point(position=position_jitter(width=0,height=.04))+    #Jittering helps reduce overlap
stat_smooth(method='glm',family=binomial(link='probit'))+
scale_y_continuous('Probability of dead crown')+
scale_x_continuous('Distance to nearest dead crown')




The blue line is just the model of Living~Dead_Dist; this version doesn't control Bay_Dist. The black dots are the observations. I've only jittered them vertically (they're all 0 or 1, and it's not hard to tell which) because their horizontal position is meaningful. The dark area is the 95% confidence band.
A: First how to interpret the results. In the Probit model, you model the probability of success $\pi=\Phi(x'\beta)$, where $\Phi$ is the cumulative normal distribution. In other words $\Phi^{-1}(\pi)=x'\beta$. So the interpretation of the regression coefficients are a little bit weird! For example let's look at the summary of your model
> summary(myprobit)

Call:
glm(formula = Living ~ Bay_Dist + Dead_Dist, family = binomial(link = "probit"), 
    data = mydata1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5314  -0.8864  -0.4853   0.9984   1.6558  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.557554   1.509486   1.032    0.302
Bay_Dist    -0.002704   0.003895  -0.694    0.488
Dead_Dist   -0.017399   0.014683  -1.185    0.236

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13.460  on 9  degrees of freedom
Residual deviance: 11.707  on 7  degrees of freedom
AIC: 17.707

Number of Fisher Scoring iterations: 5

Here the estimated coefficient for Dead_Dist is -0.017399. What does that means? Well based on what I wrote above (i.e. $\Phi^{-1}(\pi)=x'\beta$) it means that if you increase the nearest dead crown by one unit (while keeping Bay_Dist variable unchanged) then $\Phi^{-1}$ of the probability that we observe a dead crown deceases by 0.017399. This is as I said a little bit weird but we need to live with it! 
Now lets find the estimated probabilities and try to plot them in 3D.
> p=predict(myprobit, type = "response" )
> library(scatterplot3d)
> scatterplot3d(Bay_Dist,Dead_Dist,p,pch=19) 


