The expected number of missing values is the mean of the probabilities
$$f_X(\alpha, \beta) = \frac{1}{n}\sum_i \frac{1}{1 + \exp(-\alpha-X_i\beta^\prime)}$$
where I have written $\alpha$ for the constant term, $\beta$ as a row vector of $d$ coefficients, and followed the convention of arranging the variable values in columns of the design matrix $X$ with $d$ columns (say) and $n$ rows so that $X_i$ is the $i^\text{th}$ row, representing one case.
This is a smooth function of the $d+1$ parameters $\alpha, \beta$ whose range is the full interval $(0,1)$. Therefore, given a target proportion $p$ there exists a $d$-dimensional manifold $f_X^{-1}(p)$ of parameter values for which
$$f_X(\alpha, \beta) =p.$$
I interpret the question as a request for a way to find $\alpha, \beta$ that satisfy this equation.
One way to find a definite solution is to pick a one-dimensional path $t\to\gamma(t)$ in the $(\alpha, \beta)$ space and compute its intersection with $f^{-1}(p)$, which amounts to finding a univariate root. When the path is smooth and intersects $f^{-1}(p)$ transversally, a standard root-finding method such as Newton-Raphson will usually work very well to solve the one-variable equation
$$f(\gamma(t)) - p = 0.$$
For instance, among the infinitely many choices available, you might fix the coefficients $\beta$ and vary $\alpha$, giving the path $t\to (t, \beta)$, or you might fix all the coefficients initially and use the path $t\to (t\alpha, t\beta)$ to scale them uniformly.
As a worked example of the first method (varying only the intercept), consider the following R code. After creating some sample data x and choosing a fixed set of coefficients beta, it sets up a function f to compute $t\to f(\gamma(t))$ and then invokes a root-finding routine uniroot to find zeros of $f(\gamma(t))-p$. Any such zero will be a usable value of $\alpha$, thereby fully specifying the missingness model (which can later be used to generate simulated datasets).
n <- 10; d <- 3 # Set the dimensions of the data
set.seed(17) # Create a reproducible starting point
x <- matrix(rnorm(n*d), n, d) # Create independent variables
#
# Determine a model with a given expected proportion of responses
# missing at random.
#
beta <- 1:d # Fix the coefficients at any desired values
f <- function(t) { # Define a path through parameter space
sapply(t, function(y) mean(1 / (1 + exp(-y -x %*% beta))))
# (sapply makes this function vectorizable)
}
#
# Find parameters (alpha, beta) yielding any specified proportions `p`.
#
p <- seq(.1, .9, .2) # Specify the target proportions
results <- sapply(p, function(p) {
alpha <- uniroot(function(t) f(t) - p, c(-1e6, 1e6), tol = .Machine$double.eps^0.5)$root
c(alpha, f(alpha))})
dimnames(results) <- list(c("alpha", "f(alpha)"), p=p)
print(results)
The output tabulates some solutions for various values of $p$ along with a check of their validity:
p
0.1 0.3 0.5 0.7 0.9
alpha -5.102146 -3.45776 -2.251828 -0.7976692 1.27534
f(alpha) 0.100000 0.30000 0.500000 0.7000000 0.90000
Modify p as desired. Change f and/or beta to search for other sets of parameters.
