Using monte carlo integration with bisection method to find a true value in R

I try to (re)produce the results of Austin, PC. (2010).. The question is to obtain an approximated $\beta$, which leads to a specific risk difference(difference of proportions) $p_0 - p_1 = 0.02$, where $p_0$ & $p_1$ are proportions of subjects in the control and trt groups that have the response $Y = 1$, respectively. There are 4 covariates, 1 treatment indicator, and 1 outcome variable:$x_1, x_2, x_3,$ and $x_4$, which

set.seed(123)
N = 10000
x1=rnorm(N);x2=rnorm(N);
x3=rbinom(N,1,0.5);x4=rbinom(N,1,0.5)


where $N$ represents the number of subjects in each dataset and there are total $s = 1000$ simulated datasets. For each subject, we compute the $Pr(Y = 1)$ if a subject is in the trt $(t=1)$ or control $(t=0)$is

$p_{i,1} = \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i} + \beta_{(k)})}}$, or $p_{i,0} = \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i})}}$

where the values of $\alpha_0=-1, \alpha_1=\alpha_2=\alpha_3=\alpha_4=1$ are fixed, but $\beta_{(k)}$ is approximated by monte carlo integrations with $k$th iterative process by using bisection method : $\bar{p}_1 = \quad \idotsint_{x_1} \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i} + \beta_{(k)})}} f(x_1)...f(x_4)dx_1...dx_4$ and $\bar{p}_0 = \quad \idotsint_{x_1} \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i})}} f(x_1)...f(x_4)dx_1...dx_4$.

$\bar{p}_0^{mc(k)} = \frac{1}{10000}\sum_{i=1}^{10000}p_{i,0}$ and $\bar{p}_1^{mc(k)} = \frac{1}{10000}\sum_{i=1}^{10000}p_{i,1}$, where $mc$ represents monte carlo. The empirical marginal risk difference of $s^{th}$ dataset is $\gamma_{(s)}^{(k)} = \bar{p}_0^{mc(k)} - \bar{p}_1^{mc(k)}$. We repeated over 1000 datasets, the mean empirical marginal risk difference after $k^{th}$ iterative process (bisection method) is $\gamma^{(k)} = \frac{1}{1000}\sum_{s=1}^{1000}\gamma_{(n)}^{(k)} =0.02$

How can I write a R code to obtain approximated $\beta$ for $p_0 - p_1 = \gamma^{(k)}=0.02$? many thanks in advance.

• Can't these integrals be computed exactly (symbolically)? That would solve the problem of coping with Monte-Carlo integrations. Iterated integrals of these inverse logit functions can be expressed in terms of exponentials and polylogarithms. Barring that, why not use an adaptive quadrature (numerical integration) routine? For the bisection (or any root-finding) method, you only need to update the integrals, which involves integration over small regions: that ought to be quite accurate and fast.
– whuber
Commented Jan 20, 2012 at 22:44
• To add to whuber's comment, integrate is the function you use in R to perform numerical integration. Commented Jan 21, 2012 at 11:04
• I also think to use "integrate" but it does not make sense to me when covariates are binomial. I still do not how can binomial variable can use integrate. thanks
– Tu.2
Commented Jan 21, 2012 at 17:19

You should edit your question to make it clearer; your notations and your aim are very confused. I think I guessed what you want, but this only a guess...

First step, given $\beta$, given

• the values of the $\alpha_i$'s (here $\alpha_i = 1$),
• the distribution of the $x_i$’s

compute an approximation of $\gamma = p_1 - p_0$ by Monte Carlo (we compute a point estimation and a confidence interval).

MC_gamma <- function(beta, alpha = c(-1,1,1,1,1) , N = 100)
{
x1 <- rnorm(N);
x2 <- rnorm(N);
x3 <- rbinom(N, 1, 0.5);
x4 <- rbinom(N, 1, 0.5);

a <- alpha[1] + alpha[2]*x1 +  alpha[3]*x2 + alpha[4]*x3 + alpha[5]*x4;
p1 <- 1/(1+exp(-a-beta));
p0 <- 1/(1+exp(-a));

dp <- p1-p0;
mu  <- mean(dp);
std <- sd(dp);
r <- c(mu,std,mu-1.96*std/sqrt(N),mu+1.96*std/sqrt(N));
names(r) <- c("gamma", "sd", "lo", "hi");
return(r);
}


A graph (we plot the estimation and the bounds of the CI):

t <- seq(-1,1,length=101)
sapply(t, MC_gamma) -> u
matplot(t,t(u[c(1,3,4),]), type="l", lty=1)


You can tell that the $\beta$ leading to $\gamma = 0.02$ is around $0.25$. On this graph, obtained with only 100 iterations, the size of the confidence interval is too big, but with $N = 100\,000$ iterations we can try to apply bisection method, in a second step:

bisection <- function(F,a,b,epsilon = 1e-2, max=100)
{
N <- 1;
Fa <- F(a);
while(abs(a-b)>epsilon & (N <- N+1)<max)
{
c <- (a+b)/2;
Fc <- F(c);
if( F(c)*F(a) < 0)
{
b<-c;
}
else
{
a<-c;
Fa <- Fc;
}
}
return((a+b)/2);
}

> set.seed(1)
> bisection( function(x) MC_gamma(x,N=1e5)[1]-0.02, -1, 1, 1e-5)
[1] 0.1171837
> MC_gamma(0.117,N=1e6)
gamma          sd          lo          hi
0.020095504 0.008059177 0.020079708 0.020111300


It is a bit difficult to apply bisection method to a non-deterministic function. I think there is a need for an improvement, e.g. by dynamically increasing the number of iterations when approaching the solution...

Edit: an obvious solution to this non-deterministic issue

An obvious solution is to draw $x_1, x_2, x_3, x_4$ only once, and to use them for all computations of $\gamma$ values.

MC_gamma2 <- function(beta, x, alpha)
{
a <- as.vector(x %*% alpha) ;
p1 <- 1/(1+exp(-a-beta));
p0 <- 1/(1+exp(-a));

dp <- p1-p0;
mu  <- mean(dp);
std <- sd(dp);
r <- c(mu,std,mu-1.96*std/sqrt(N),mu+1.96*std/sqrt(N));
names(r) <- c("gamma", "sd", "lo", "hi");
return(r);
}

> N <- 1e6
> set.seed(1)
> x <- matrix( c( rep(1,N), rnorm(2*N), rbinom(2*N,k,0.5)), ncol=5)
> bisection( function(t) MC_gamma2(t,x,c(-1,1,1,1,1))[1]-0.02, -1, 1, 1e-5)
[1] 0.1163597
> MC_gamma(0.11636,N=1e6)
gamma         sd         lo         hi
0.01999671 0.00800599 0.01998102 0.02001240


This is much better, it saves time and you have less issues with a potential "sawteeth" behaviour of the function. Of course you may want to repeat this operation a number of times and take the mean of the results.

Can you please tell me if this is the solution recommended in the paper? (I don’t have access to it)

Edit 2: (minor modifications of the above, and an other application suggested in the comments)

For 5 variables $x_i$ all normally distributed, we just have to write:

> set.seed(1)
> x <- matrix( c( rep(1,N), rnorm(5*N)), ncol=6)
> bisection( function(t) MC_gamma2(t,x,c(-1,1,1,1,1,1))[1]-0.02, -1, 1, 1e-5)
[1] 0.1496468

• Elvis, many thanks. In Austin PC (2010) paper, he has 5 different simulated scenarios. One of them has 5 covariates and each of covariates is standard normally distributed. He simulated 1000 datasets and each dataset has 10000 subjects. His beta = 0.858468 for gamma = 0.02. I don't know how he got this but I think your answer is based on one dataset and it should be very close to his answer. again many thanks
– Tu.2
Commented Jan 21, 2012 at 2:01
• Hi Elvis, I used your method to check with Austin's result when x1,...,x5 are standard normal distributed with N = 10000, for 1000 simulated datasets. My answer is around -0.1496281, which is far away from 0.858468. However, your method makes more senses to me than Austin's method. Thanks
– Tu.2
Commented Jan 21, 2012 at 7:56
• @Tu.2 I think there is a problem with the notations. $\beta > 0$ leads to values $p_1 > p_0$ and negative values of $\gamma = p_0 - p_1$, so for $\beta > 0$ you can’t have $\gamma = 0.02$. Maybe you switched the indices at some point, but this shouldn’t be a real problem (just a matter of reversing signs). You precised that the binomial is with size $k=1$ which makes sense (I didn’t get that these variables were treatments). I’ll edit my answer accordingly. Commented Jan 21, 2012 at 10:17
• @Tu.2 I changed the definition of $\gamma$ to $\gamma = p_1-p_0$... I think it is more natural. I added the application suggested in your second comment, and I have the same result (up to the sign) than you. I am worried by the differences with the results in the paper, I fear we are missing something. Commented Jan 21, 2012 at 10:46
• @Tu.2 Can you send me a copy of it to elvis.bistalanouille [at] gmail.com ? I am really curious to see what’s going wrong. Commented Jan 21, 2012 at 10:52