logit <- function(x){log(x/(1-x))}

## variance in the distribution of r
v <- 3

conditional.prob <- function(x,r){
  exp(x+r)/(1+exp(x+r))
}

plot(0,0,type = "n",
     xlim = c(-3,3),
     ylim = c(-10,10),
     ylab = "Probability, logit scale",
     xlab = "x"
     )

random.effects <- c(-v,v,0, qnorm(c(0.025, 0.975, 0.05,0.95), 0,v)) 
for (r in random.effects)
  ## same as plotting x + r directly
  curve(logit(conditional.prob(x, r)), add = T)

## numeric
marginal.prob <- function(x){
  integrate(function(r) conditional.prob(x, r) * dnorm(r, 0, v),
            -3,3)$value 
}
marginal.prob <- Vectorize(marginal.prob)
curve(logit(marginal.prob(x)), add = T, lty = "32")