Is it possible to model the conditional expectation of a binary outcome using an additively-separable link function? Logit and probit link functions aren't additively separable.  So, fitting a model using these link functions implies that the effect of one predictor in determining the outcome is not independent of the values of the other predictors.  This might make sense in a lot of contexts, but might be suspect in others.    
Take a regression of going to college (0,1) on income and gender, and suppose for a moment that the effect of income is the same for both men and women, but that men are more likely than women to go to college.  In a logit (or probit) model, the marginal effect of income will be different for a man versus a woman, because of a negative coefficient on the female dummy -- the model is biased.  If you included an interaction, it would seem like the effect of income was weaker for women than men, because the model would be overcompensating for the mis-specified link function.  This would be wrong, and misleading.
Of course, linear probability models are wrong for different reasons.
So, what sort of options exist for dealing with mis-specified link functions?  And is there an additively-seperable model for binary data?
 A: The logistic regression estimates conditional probabilities
$$\hat P(y_i\mid x_i, z_i) = \Lambda(\hat \alpha + \hat \beta x_i +\hat \gamma z_i)$$
with $y_i$ "going to college", $x_i$ "income" and $z_i$ "gender", with $z_i=1$ representing "male". Let $i$ index the observations related to males, and $j$ the observations related to females. We will write for compactness $\Lambda(\hat \alpha + \hat \beta x_i +\hat \gamma z_i) = \hat \Lambda_i$,and analogously for females, $\hat \Lambda_j$. Based on the question, we examine the case where we have $x_i =x_j ,\; \hat \Lambda_i > \hat \Lambda_j$, which implies $\hat \gamma >0$.
Since "income" is (treated as) a continuous regressor, the "marginal effect of income" for males is
$$m_i(x_i)\equiv \frac {\partial \hat P(y_i\mid x_i, z_i=1)}{\partial x_i} = \hat \Lambda_i(1-\hat \Lambda_i)\hat \beta$$
while for females is
$$m_j(x_j)\equiv \frac {\partial \hat P(y_j\mid x_j, z_j=0)}{\partial x_j} = \hat \Lambda_j(1-\hat \Lambda_j)\hat \beta$$
The expression $\Lambda (1-\Lambda)$ is maximized at $ \Lambda = 1/2$. It follows that for $x_i=x_j$, in general we will have $m_i(x_i) \neq m_j(x_j)$, although we won't necessarily have $m_i(x_i) > m_j(x_j)$ -this will depend on whether, for the given level of income, $\Lambda_i$ is closer to $1/2$ than $\Lambda_j$. 
Now, $\Lambda (1-\Lambda)$ is the conditional variance, and it is the sole source of the difference of the two "marginal effects". And naturally, we expect that $\hat \Lambda_i(1-\hat \Lambda_i) \neq \hat \Lambda_j(1-\hat \Lambda_j)$, because in the one the variable $z$ is set equal to $1$, while in the other it is set equal to $0$. In a sense, we have here two different conditional distributions. Why should we expect that they should have the same variance?  
Why the variance of the random indicator function "going to college or not" should be the same for males and for females, conditional on the same level of income?  The conditional probability itself, the expected value of the indicator function, is not equal, why should the variance suddenly become equal?  
We did perform something familiar, the operation of differentiation on some mathematical expression with respect to some variable. In mathematical terms, it does express mathematically how the expression changes as the variable changes, something we do call the "marginal effect". But it requires per-case interpretation as to what it reflects and represents in the context of the specific framework on which we have mapped and applied our mathematical tools.
A: One way to think of a GLM is to specify a linear predictor, which must be a linear function of covariates:
$$
\eta_i = X_i \beta,
$$
an inverse-link $g$, whose input is linear predictor and output is the conditional mean:
$$
\mu_i = E(Y_i | X_i) = g \left( \eta_i \right),
$$
and a variance function, which depends on only on the conditional mean, and is determined by the parametric family of the conditional distribution $Y_i | \eta_i$.  For a Bernoulli variable, the log-likelihood of $ Y | X$ will be
$$
\begin{align*}
L(Y|X) &=  \sum_i \log \left( \mu_i^{y_i} \left( 1-\mu_i \right)^{1-y_i} \right)  \\
 & =  \sum_i y_i \mbox{logit}(\mu_i) + \log ( 1 - \mu_i) \\
& = \sum_i y_i \mbox{logit}(g(\eta_i)) + \log ( 1 - g(\eta_i)) \quad \mathbf{(1)}
\end{align*}
$$
So you can see that the cancellation that occurs is one reason why it's attractive to take $g$ to be the inverse logit function $\exp(\eta_i)/(1+\exp(\eta_i))$. But otherwise, sure, you can certainly take $g$ to be any function you'd like, and in particular
$$g(\eta_i) = \eta_i \quad \mathbf{(2)}$$
which of course would give you linearity of $\mu_i$ in $\eta_i$.
You just can't solve it with off-the-shelf software (that I can find with a cursory google search).  It might not be an easy optimization problem, either, since it's not clear to me that the parameter space need be convex in $\beta$.
Another issue is that although I conjecture that for any design $\mathbf X$ there exists $\beta$ such that the likelihood exists and is finite (eg, as long as there's an intercept, you can always take rest of it to be zero), predictions out of sample may not always exist.
Non-parametric estimation of the link
It is also possible to estimate the model
$$\mu_i = h( X_i \beta)$$
for an unknown function $h$. This is an example of a single index projection pursuit model. (Contrast it to the generalized additive model, which assumes $\mu_i = h_1(x_{i1}) + h_2(x_{i2}) + \dots + h_p(x_{ip})$ for unknown functions $h_1, \dots , h_p$.)
The intro of "Nonparametric Estimation of the Link Function
Including Variable Selection" (Tutz and Petry 2010) does a reasonable job of providing background references to this.
