Difference between logit and probit models What is the difference between Logit and Probit model?
I'm more interested here in knowing when to use logistic regression, and when to use Probit.
If there is any literature which defines it using R, that would be helpful as well.
 A: One of the most well-known difference between logit and probit is the (theoretical) regression residuals distribution: normal for probit, logistic for logit (please see: Koop G. An Introduction to Econometrics Chichester, Wiley: 2008: 280).
A: A standard linear model (e.g., a simple regression model) can be thought of as having two 'parts'. These are called the structural component and the random component. For example:
$$
Y=\beta_0+\beta_1X+\varepsilon  \\
\text{where }  \varepsilon\sim\mathcal{N}(0,\sigma^2)
$$
The first two terms (that is, $\beta_0+\beta_1X$) constitute the structural component, and the $\varepsilon$ (which indicates a normally distributed error term) is the random component. When the response variable is not normally distributed (for example, if your response variable is binary) this approach may no longer be valid. The generalized linear model (GLiM) was developed to address such cases, and logit and probit models are special cases of GLiMs that are appropriate for binary variables (or multi-category response variables with some adaptations to the process). A GLiM has three parts, a structural component, a link function, and a response distribution. For example:
$$
g(\mu)=\beta_0+\beta_1X
$$
Here $\beta_0+\beta_1X$ is again the structural component, $g()$ is the link function, and $\mu$ is a mean of a conditional response distribution at a given point in the covariate space. The way we think about the structural component here doesn't really differ from how we think about it with standard linear models; in fact, that's one of the great advantages of GLiMs. Because for many distributions the variance is a function of the mean, having fit a conditional mean (and given that you stipulated a response distribution), you have automatically accounted for the analog of the random component in a linear model (N.B.: this can be more complicated in practice).  
The link function is the key to GLiMs: since the distribution of the response variable is non-normal, it's what lets us connect the structural component to the response--it 'links' them (hence the name). It's also the key to your question, since the logit and probit are links (as @vinux explained), and understanding link functions will allow us to intelligently choose when to use which one. Although there can be many link functions that can be acceptable, often there is one that is special. Without wanting to get too far into the weeds (this can get very technical) the predicted mean, $\mu$, will not necessarily be mathematically the same as the response distribution's canonical location parameter; the link function that does equate them is the canonical link function. The advantage of this "is that a minimal sufficient statistic for $\beta$ exists" (German Rodriguez). The canonical link for binary response data (more specifically, the binomial distribution) is the logit.  However, there are lots of functions that can map the structural component onto the interval $(0,1)$, and thus be acceptable; the probit is also popular, but there are yet other options that are sometimes used (such as the complementary log log, $\ln(-\ln(1-\mu))$, often called 'cloglog'). Thus, there are lots of possible link functions and the choice of link function can be very important.  The choice should be made based on some combination of:  


*

*Knowledge of the response distribution,  

*Theoretical considerations, and  

*Empirical fit to the data.  


Having covered a little of conceptual background needed to understand these ideas more clearly (forgive me), I will explain how these considerations can be used to guide your choice of link. (Let me note that I think @David's comment accurately captures why different links are chosen in practice.)  To start with, if your response variable is the outcome of a Bernoulli trial (that is, $0$ or $1$), your response distribution will be binomial, and what you are actually modeling is the probability of an observation being a $1$ (that is, $\pi(Y=1)$). As a result, any function that maps the real number line, $(-\infty,+\infty)$, to the interval $(0,1)$ will work.  
From the point of view of your substantive theory, if you are thinking of your covariates as directly connected to the probability of success, then you would typically choose logistic regression because it is the canonical link. However, consider the following example: You are asked to model high_Blood_Pressure as a function of some covariates. Blood pressure itself is normally distributed in the population (I don't actually know that, but it seems reasonable prima facie), nonetheless, clinicians dichotomized it during the study (that is, they only recorded 'high-BP' or 'normal'). In this case, probit would be preferable a-priori for theoretical reasons. This is what @Elvis meant by "your binary outcome depends on a hidden Gaussian variable". Another consideration is that both logit and probit are symmetrical, if you believe that the probability of success rises slowly from zero, but then tapers off more quickly as it approaches one, the cloglog is called for, etc.  
Lastly, note that the empirical fit of the model to the data is unlikely to be of assistance in selecting a link, unless the shapes of the link functions in question differ substantially (of which, the logit and probit do not). For instance, consider the following simulation:  
set.seed(1)
probLower = vector(length=1000)

for(i in 1:1000){      
    x = rnorm(1000)
    y = rbinom(n=1000, size=1, prob=pnorm(x))

    logitModel  = glm(y~x, family=binomial(link="logit"))
    probitModel = glm(y~x, family=binomial(link="probit"))

    probLower[i] = deviance(probitModel)<deviance(logitModel)
}

sum(probLower)/1000
[1] 0.695

Even when we know the data were generated by a probit model, and we have 1000 data points, the probit model only yields a better fit 70% of the time, and even then, often by only a trivial amount. Consider the last iteration:  
deviance(probitModel)
[1] 1025.759
deviance(logitModel)
[1] 1026.366
deviance(logitModel)-deviance(probitModel)
[1] 0.6076806

The reason for this is simply that the logit and probit link functions yield very similar outputs when given the same inputs.  
 
The logit and probit functions are practically identical, except that the logit is slightly further from the bounds when they 'turn the corner', as @vinux stated. (Note that to get the logit and the probit to align optimally, the logit's $\beta_1$ must be $\approx 1.7$ times the corresponding slope value for the probit. In addition, I could have shifted the cloglog over slightly so that they would lay on top of each other more, but I left it to the side to keep the figure more readable.) Notice that the cloglog is asymmetrical whereas the others are not; it starts pulling away from 0 earlier, but more slowly, and approaches close to 1 and then turns sharply.  
A couple more things can be said about link functions. First, considering the identity function ($g(\eta)=\eta$) as a link function allows us to understand the standard linear model as a special case of the generalized linear model (that is, the response distribution is normal, and the link is the identity function). It's also important to recognize that whatever transformation the link instantiates is properly applied to the parameter governing the response distribution (that is, $\mu$), not the actual response data. Finally, because in practice we never have the underlying parameter to transform, in discussions of these models, often what is considered to be the actual link is left implicit and the model is represented by the inverse of the link function applied to the structural component instead. That is:
$$
\mu=g^{-1}(\beta_0+\beta_1X)
$$
For instance, logistic regression is usually represented:
$$
\pi(Y)=\frac{\exp(\beta_0+\beta_1X)}{1+\exp(\beta_0+\beta_1X)}
$$
instead of:
$$
\ln\left(\frac{\pi(Y)}{1-\pi(Y)}\right)=\beta_0+\beta_1X
$$  
For a quick and clear, but solid, overview of the generalized linear model, see chapter 10 of Fitzmaurice, Laird, & Ware (2004), (on which I leaned for parts of this answer, although since this is my own adaptation of that--and other--material, any mistakes would be my own). For how to fit these models in R, check out the documentation for the function ?glm in the base package.  
(One final note added later:) I occasionally hear people say that you shouldn't use the probit, because it can't be interpreted. This is not true, although the interpretation of the betas is less intuitive. With logistic regression, a one unit change in $X_1$ is associated with a $\beta_1$ change in the log odds of 'success' (alternatively, an $\exp(\beta_1)$-fold change in the odds), all else being equal. With a probit, this would be a change of $\beta_1\text{ }z$'s. (Think of two observations in a dataset with $z$-scores of 1 and 2, for example.) To convert these into predicted probabilities, you can pass them through the normal CDF, or look them up on a $z$-table.  
(+1 to both @vinux and @Elvis. Here I have tried to provide a broader framework within which to think about these things and then using that to address the choice between logit and probit.)
A: In addition to vinux’ answer, which already tells the most important:


*

*the coefficients $\beta$ in the logit regression have natural interpretations in terms of odds ratio;

*the probistic regression is the natural model when you think that your binary outcome depends of a hidden gaussian variable $Z = X' \beta + \epsilon\ $ [eq. 1] with $\epsilon \sim \mathcal N(0,1)$ in a deterministic manner: $Y = 1$ exactly when $Z > 0$. 

*More generally, and more naturally, probistic regression is the more natural model if you think that the outcome is $1$ exactly when some $Z_0 = X' \beta_0 + \epsilon_0$ exceeds a threshold $c$, with $\epsilon \sim \mathcal N(0,\sigma^2)$.  It is easy to see that this can be reduced to the aforementioned case: just rescale $Z_0$ as $Z = {1\over \sigma}(Z_0-c)$; it’s easy to check that equation [eq. 1] still holds (rescale the coefficients and translate the intercept). These models have been defended, for example, in medical contexts, where $Z_0$ would be an unobserved continuous variable, and $Y$ eg a disease which appears when $Z_0$ exceeds some "pathological threshold".
Both logit and probit models are only models. "All models are wrong, some are useful", as Box once said! Both models will allow you to detect the existence of an effect of $X$ on the outcome $Y$; except in some very special cases, none of them will be "really true", and their interpretation should be done with cautiousness. 
A: Below, I explain an estimator that nests probit and logit as special cases and where one can test which is more appropriate.
Both probit and logit can be nested in a latent variable model, 
$$ y_i^* = x_i \beta + \varepsilon_i,\quad \varepsilon_i \sim G(\cdot), $$
where the observed component is 
$$ y_i = \mathbb{1}(y_i^* > 0). $$
If you choose $G$ to be the normal cdf, you get probit, if you choose the logistic cdf, you get logit. Either way, the likelihood function takes the form 
$$ \ell(\beta) = y_i \log G(x_i\beta) + (1-y_i) \log[1-G(x_i\beta)].$$
However, if you are concerned about which assumption you have made, you can use the Klein & Spady (1993; Econometrica) estimator. This estimator allows you to be fully flexible in your specification of the cdf, $G$, and you could then even subsequently test the validity of normality or logisticness (?). 
In Klein & Spady, the criterion function is instead 
$$ \ell(\beta) = y_i \log \hat{G}(x_i\beta) + (1-y_i) \log[1-\hat{G}(x_i\beta)],$$
where $\hat{G}(\cdot)$ is a nonparametric estimate of the cdf, for example estimated using a Nadaraya-Watson kernel regression estimator, 
$$ \hat{G}(z) = \sum_{i=1}^N y_i \frac{ K\left( \frac{z - x_i\beta}{h} \right)}{\sum_{j=1}^N K\left( \frac{z - x_j\beta}{h} \right)}, $$
where $K$ is called the "Kernel" (typically, the Gaussian cdf or a triangular kernel is chosen), and $h$ is a "bandwidth". There are plugin values to pick for the latter but it can be a lot more complicated and it can make the outer optimization over $\beta$ more complicated if $h$ changes in every step ($h$ balances the so-called bias-variance tradeoff). 
Improvements: Ichimura has suggested that the kernel regression, $\hat{G}$, should leave out the $i$th observation; otherwise, the choice of $h$ may be complicated by a problem with over-fitting in sample (too high variance). 
Discussion: One drawback with the Klein-Spady estimator is that it may get stuck in local minima. This is because the $G$ cdf adapts to the given $\beta$-parameters. I know of several students who have tried implementing it and have had problems achieving convergence and avoiding numerical issues. Hence, it is not an easy estimator to work with. Moreover, inference on the estimated parameters is complicated by the semi-parametric specification for $G$. 
A: They are very similar.
In both models, the probability that $Y=1$ given $X$ can be seen as the probability that a random hidden variable $S$ (with a certain fixed distribution) is below a certain threshold that depends linearly on $X$ :
$$P(Y=1|X)=P(S<\beta X)$$
Or equivalently :
$$P(Y=1|X)=P(\beta X-S>0)$$
Then it's all a matter of what you choose for the distribution of $S$ :


*

*in logistic regression, $S$ has a logistic distribution.

*in probit regression, $S$ has a normal distribution.


Variance is unimportant since it is automatically compensated by multiplying $\beta$ by a constant. Mean is unimportant as well if you use an intercept.
This can be seen as a threshold effect. Some invisible outcome $E=\beta X-S$ is a linear function of $X$ with some noise $-S$ added like in linear regression, and we get a 0/1 outcome by saying:


*

*when $E>0$, outcome is $Y=1$

*when $E<0$, outcome is $Y=0$


The differences between logistic and probit lies in the difference between the logistic and the normal distributions. There ain't that much. Once adjusted, they look like it :

Logistic has heavier tail. This may impact a little how events of small (<1%) or high (>99%) probability are fitted. Practically, the difference is not even noticeable in most situations : logit and probit predict essentially the same thing. See http://scholarworks.rit.edu/cgi/viewcontent.cgi?article=2237&context=article
"Philosophically", logistic regression can be justified by being equivalent to the principle of maximum entropy : http://www.win-vector.com/blog/2011/09/the-equivalence-of-logistic-regression-and-maximum-entropy-models/
In terms of calculation : logistic is simpler since the cumulative distribution of the logistic distribution has a closed formula unlike the normal distribution. But normal distributions have good properties when you go to multi-dimensional, this is why probit is often preferred in advanced cases.
A: Regarding your statement
I'm more interested here in knowing when to use logistic regression, and when to use probit
There are already many answers here that bring up things to consider when choosing between the two but there is one important consideration that hasn't been stated yet: When your interest is in looking at within-cluster associations in binary data using mixed effects logistic or probit models, there is a theoretical grounding for preferring the probit model. This is, of course, assuming that there is no a priori reason for preferring the logistic model (e.g. if you're doing a simulation and know it to be the true model). 
First, To see why this is true first note that both of these models can be viewed as thresholded continuous regression models. As an example consider the simple linear mixed effects model for the observation $i$ within cluster $j$: 
$$ y^{\star}_{ij} = \mu + \eta_{j} + \varepsilon_{ij} $$ 
where $\eta_j \sim N(0,\sigma^2)$ is the cluster $j$ random effect and $\varepsilon_{ij}$ is the error term. Then both the logistic and probit regression models are equivalently formulated as being generated from this model and thresholding at 0: 
$$
y_{ij} = 
\begin{cases} 1 & \text{if} \ \ \ y^{\star}_{ij}≥0\\
\\
0 &\text{if} \ \ \ y^{\star}_{ij}<0
\end{cases}
$$
If the $\varepsilon_{ij}$ term is normally distributed, you have a probit regression and if it is logistically distributed you have a logistic regression model. Since the scale is not identified, these residuals errors are specified as standard normal and standard logistic, respectively. 
Pearson (1900) showed that that if multivariate normal data were generated and thresholded to be categorical, the correlations between the underlying variables were still statistically identified - these correlations are termed polychoric correlations and, specific to the binary case, they are termed tetrachoric correlations. This means that, in a probit model, the intraclass correlation coefficient of the underlying normally distributed variables: 
$$ {\rm ICC} = \frac{ \hat{\sigma}^{2} }{\hat{\sigma}^{2} + 1 } $$
is identified which means that in the probit case you can fully characterize the joint distribution of the underlying latent variables. 
In the logistic model the random effect variance in the logistic model is still identified but it does not fully characterize the dependence structure (and therefore the joint distribution), since it is a mixture between a normal and a logistic random variable that does not have the property that it is fully specified by its mean and covariance matrix. Noting this odd parametric assumption for the underlying latent variables makes interpretation of the random effects in the logistic model less clear to interpret in general. 
A: An important point that has not been addressed in the previous (excellent) answers is the actual estimation step. Multinomial logit models have a PDF that is easy to integrate, leading to a closed-form expression of the choice probability. The density function of the normal distribution is not so easily integrated, so probit models typically require simulation. So while both models are abstractions of real world situations, logit is usually faster to use on larger problems (multiple alternatives or large datasets).
To see this more clearly, the probability of a particular outcome being selected is a function of the $x$ predictor variables and the $\varepsilon$ error terms (following Train) 
$$
P = \int I[\varepsilon > -\beta'x] f(\varepsilon)d\varepsilon
$$
Where $I$ is an indicator function, 1 if selected and zero otherwise.
Evaluating this integral depends heavily on the assumption of $f(x)$. In a logit model, this is a logistic function, and a normal distribution in the probit model. For a logit model, this becomes
$$
P=\int_{\varepsilon=-\beta'x}^{\infty} f(\varepsilon)d\varepsilon\\
= 1- F(-\beta'x) = 1-\dfrac{1}{\exp(\beta'x)}
$$ 
No such convenient form exists for probit models.
A: @Benoit Sanchez and @gung's graphs emphasize how little there is to distinguish the link functions, except with very large numbers of observations and/or in the extreme tails.  The conversion ${\rm logit}(p) = 1.77\ {\rm probit}(p)$ never has an error of more than $0.1$ over the range $0.1 \le p \le 0.9$.  Over $0.01 \le p \le 0.99$, the conversion is still a good approximation.  Using such a conversion factor, one can get a good approximation to the odds ratio from a model that has a probit link.
The graph uses the following R code:
trellis.par.set(clip=list(panel='off'))
p <- (1:99)/100
logit <- make.link('logit')$linkfun
probit <- make.link('probit')$linkfun
gph <- lattice::xyplot(logit(p)~probit(p), type=c('p','r'),
             scales=list(x=list(alternating=1, tck=c(0.6,0))))
b <- coef(lm(logit(p) ~ 0 + probit(p)))  ## Slope of line
gph1 <- update(gph, main = list(paste("Slope of line is", round(b,2)), cex=1))
probs <- c(.01,.04,.1, .25,.5,.75,.9,.96,.99)
gph1+latticeExtra::layer(panel.axis(side='top', at=probit(probs), 
                         labels=paste(probs),outside=T, rot=0))



A: Directly from 'Discrete Choice Methods with Simulation', from Kenneth Train, in the context of discrete choice models:

The logit model is limited in three important ways. It cannot represent
random taste variation. It exhibits restrictive substitution patterns due to
the IIA property. And it cannot be used with panel data when unobserved
factors are correlated over time for each decision maker. GEV models
relax the second of these restrictions, but not the other two. Probit models
deal with all three. They can handle random taste variation, they allow
any pattern of substitution, and they are applicable to panel data with
temporally correlated errors.
The only limitation of probit models is that they require normal distributions for all unobserved components of utility. In many, perhaps most
situations, normal distributions provide an adequate representation of
the random components. However, in some situations, normal distributions are inappropriate and can lead to perverse forecasts.

A: They mainly differ in the link function.
In Logit:
$\Pr(Y=1 \mid X) = [1 + e^{-X'\beta}]^{-1} $
In Probit:
$\Pr(Y=1 \mid X) = \Phi(X'\beta)$   (Cumulative standard normal pdf)
In other way, logistic has slightly flatter tails. i.e the probit curve approaches the axes more quickly than the logit curve.
Logit has easier interpretation than probit. Logistic regression can be interpreted as modelling log odds (i.e those who smoke >25 cigarettes a day are 6 times more likely to die before 65 years of age). Usually people start the modelling with logit. You could use the likelihood value of each model to decide for logit vs probit.
A: What I am going to say in no way invalidates what has been said thus far. I just want to point out that probit models do not suffer from IIA (Independence of Irrelevant alternatives) assumptions, and the logit model does. 
To use an example from Train's excellent book. If I have a logit that predicts whether I am going to ride the blue bus or drive in my car, adding red bus would draw from both car and blue bus proportionally. But using a probit model you can avoid this problem. In essence, instead of drawing from both proportionally, you may draw more from blue bus as they are closer substitutes. 
The sacrifice you make is that there is no closed form solutions, as pointed out above. Probit tends to be my goto when I am worried about IIA issues. That's not to say that there aren't ways to get around IIA in a logit framework (GEV distributions). But I've always looked at these sorts of models as a clunky way around the problem. With the computational speeds that you can get, I would say go with probit.
A: I offer a practical answer to the question, that only focuses on "when to use logistic regression, and when to use probit", without getting into statistical details, but rather focusing on decisions based on statistics. The answer depends on two main things: do you have a disciplinary preference, and do you only care for which model better fits your data?
Basic difference
Both logit and probit models provide statistical models that give the probability that a dependent response variable would be 0 or 1. They are very similar and often given practically idential results, but because they use different functions to calculate the probabilities, their results are sometimes slightly different.
Disciplinary preference 
Some academic disciplines generally prefer one or the other. If you are going to publish or present your results to an academic discipline with a specific traditional preference, then let that dictate your choice so that your findings would be more readily acceptable. For example (from Methods Consultants), 

Logit – also known as logistic regression – is more popular in health
  sciences like epidemiology partly because coefficients can be
  interpreted in terms of odds ratios. Probit models can be generalized
  to account for non-constant error variances in more advanced
  econometric settings (known as heteroskedastic probit models) and
  hence are used in some contexts by economists and political
  scientists.

The point is that the differences in results are so minor that the ability for your general audience to understand your results outweigh the minor differences between the two approaches.
If all you care about is better fit...
If your research is in a discipline that does not prefer one or the other, then my study of this question (which is better, logit or probit) has led me to conclude that it is generally better to use probit, since it almost always will give a statistical fit to data that is equal or superior to that of the logit model. The most notable exception when logit models give a better fit is in the case of "extreme independent variables" (which I explain below).
My conclusion is based almost entirely (after searching numerous other sources) on Hahn, E.D. & Soyer, R., 2005. Probit and logit models: Differences in the multivariate realm. Available at: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.329.4866&rep=rep1&type=pdf. Here is my summary of the practical decision conclusions from this article concerning whether logit versus probit multivariate models provide a better fit to the data (these conclusions also apply to univariate models, but they only simulated effects for two independent variables):


*

*In most scenarios, the logit and probit models fit the data equally well, with the following two exceptions.

*Logit is definitely better in the case of "extreme independent variables". These are independent variables where one particularly large or small value will overwhelmingly often determine whether the dependent variable is a 0 or a 1, overriding the effects of most other variables. Hahn and Soyer formally define it thus (p. 4):

An extreme independent variable level involves the conﬂuence of three
  events. First, an extreme independent variable level occurs at the
  upper or lower extreme of an independent variable. For example, say
  the independent variable x were to take on the values 1, 2, and 3.2.
  The extreme independent variable level would involve the values at x =
  3.2 (or x = 1). Second, a substantial proportion (e.g., 60%) of the total n must be at this level. Third, the probability of success at
  this level should itself be extreme (e.g., greater than 99%).



*

*Probit is better in the case of "random effects models" with moderate or large sample sizes (it is equal to logit for small sample sizes). For fixed effects models, probit and logit are equally good. I don't really understand what Hahn and Soyer mean by "random effects models" in their article. Although many definitions are offered (as in this Stack Exchange question), the definition of the term is in fact ambiguous and inconsistent. But since logit is never superior to probit in this regard, the point is rendered moot by simply choosing probit.


Based on Hahn and Soyer's analysis, my conclusion is to always use probit models except in the case of extreme independent variables, in which case logit should be chosen. Extreme independent variables are not all that common, and should be quite easy to recognize. With this rule of thumb, it doesn't matter whether the model is a random effects model or not. In cases where a model is a random effects model (where probit is preferred) but there are extreme independent variables (where logit is preferred), although Hahn and Soyer didn't comment on this, my impression from their article is that the effect of extreme independent variables are more dominant, and so logit would be preferred.
