title: Optimize function in r
body: <p>Here is my code:</p>
  <pre><code>cee = abs(qnorm(.5*0.1)) # Bonferroni threshold for achieving study-wide significance = 0.1
  p.value = (simAll %&gt;% select(&quot;p.value&quot;))
  p.value1 &lt;- as.numeric(unlist(p.value))
  # we use &quot;cee&quot; so R does not get confused with the function 'c'
  betahat = log(OR) # Reported OR  
  z = sign(betahat)*abs(qnorm(0.5*p.value1)) # Reported p-value = 5.7e-4, which we convert to a z-value
  
  ###################################################
  #              THE PROPOSED APPROACH              #
  ###################################################
  
  se = betahat/z # standard error of betahat
  
  mutilde1 = optimize(f=conditional.like,c(-20,20),maximum=T,z=z,cee=cee)$maximum # the conditional mle
  </code></pre>
  <p>The p.value is the p-values for 1000 simulations, same as OR, for the &quot;se“ part, I can get 1000 different se values there. But for the mutilde1 line, there is an error exist: &quot;Error in optimize(f = conditional.like, c(-20, 20), maximum = T, z = z,  :
  invalid function value in 'optimize'&quot;
  How can I fix the issue?
  The <code>conditional.like()</code> function:</p>
  <pre><code>conditional.like=function(mu,cee,z){
    like=dnorm(z-mu)/(pnorm(mu-cee)+pnorm(-cee-mu))
    return((abs(z)&gt;cee)*like) }
  </code></pre>
  <p>The simALL is a table looks like this (total 1000 lines)：</p>
  <pre><code># A tibble: 1,000 x 6
        id term  estimate std.error statistic  p.value
     &lt;int&gt; &lt;chr&gt;    &lt;dbl&gt;     &lt;dbl&gt;     &lt;dbl&gt;    &lt;dbl&gt;
   1     1 .x     0.226       0.127    1.78   0.0747  
   2     2 .x     0.137       0.127    1.08   0.280   
   3     3 .x     0.304       0.127    2.38   0.0171  
   4     4 .x     0.497       0.128    3.87   0.000111
  </code></pre>
  <p>OR (total 1000 lines):</p>
  <pre><code>&gt; OR
     [1] 1.5537098 1.0939850 1.4491432 1.6377551 1.1646904 1.3387534 1.6377551 1.5009351 1.7918552
  </code></pre>
  <p>Also, here is my overall code:</p>
  <pre><code>library(tidyverse)
  library(broom)
  # create a tibble with an id column for each simulation and x wrapped in list()
  sim &lt;- tibble(id = 1:1000,
                x = list(rbinom(1000,1,0.5))) %&gt;% 
    # to generate z, pr, y, k use map and map2 from the purrr package to loop over the list column x
    # `~ ... ` is similar to `function(.x) {...}`
    # `.x` represents the variable you are using map on
    mutate(z  = map(x, ~ log(1.3) * .x), 
           pr = map(z, ~ 1 / (1 + exp(-.x))),
           y  = map(pr, ~ rbinom(1000, 1, .x)),
           k  = map2(x, y, ~ glm(.y ~ .x, family=&quot;binomial&quot;)),
           # use broom::tidy to get the model summary in form of a tibble
           sum = map(k, broom::tidy)) %&gt;% 
    # select id and sum and unnest the tibbles
    select(id, sum) %&gt;% 
    unnest(cols = c(sum))
  simOR &lt;- sim %&gt;%
    # drop the intercepts and every .x with a p &lt; 0.05
    filter(term !=&quot;(Intercept)&quot;,
           p.value &lt; 0.05)
  
  sim
  j1=exp(simOR %&gt;% select(&quot;estimate&quot;))
  OR1=as.numeric(unlist(j1))
  mean(OR1)
  
  simAll &lt;- sim %&gt;% 
    filter(term !=&quot;(Intercept)&quot;)
  
  
  j &lt;- exp(simAll %&gt;% select(&quot;estimate&quot;))
  OR2 &lt;- as.numeric(unlist(j))
  mean(OR2)
  
  simOR2 &lt;- sim %&gt;%
    filter(term !=&quot;(Intercept)&quot;,
           p.value &lt; 0.005)
  
  j2 &lt;- exp(simOR2 %&gt;% select(&quot;estimate&quot;))
  OR3 &lt;- as.numeric(unlist(j2))
  mean(OR3)
  #op &lt;- par(mfrow = c(3, 1))
  hga=hist(OR2, main = NULL, freq = T, breaks = 10) #OR2:Overall OR
  hgb=hist(OR1, freq = T,col=2,breaks=10, main=&quot;OR:p-value&lt;0.05&quot;) #OR1:p-value&lt;0.05
  hgc=hist(OR3, freq = T,col=2,breaks=10, main=&quot;OR:p-value&lt;0.005&quot;) #OR3:p-value&lt;0.005
  plot(hga,col=rgb(0,1,0,0.5),main = &quot;OR&quot;,xlim=c(0.8,2),ylim=c(0,250))
  plot(hgb, add = TRUE,col=rgb(0,0,0.8,0.5),xlim=c(0.8,2),ylim=c(0,250))
  plot(hgc, add = TRUE,col=rgb(1,0,0,0.5),xlim=c(0.8,2))
  abline(v = mean(OR2), lwd = 4, col = 3)
  abline(v = mean(OR3), lwd = 4, col=2)
  text(1.65,240,&quot;1.31&quot;,col=1)
  arrows(1.5,240,1.31,240,length=0.1,col=1,lwd=2)
  abline(v = mean(OR1), lwd = 4, col=4)
  text(2.1,220,&quot;1.43&quot;,col=4)
  arrows(1.98,220,1.43,220,length=0.1,col=4,lwd=2)
  text(2.1,220,&quot;1.55&quot;,col=2)
  arrows(1.98,220,1.55,220,length=0.1,col=2,lwd=2)
  
  #########################################
  ## THE FUNCTIONS BELOW ARE USED TO OBTAIN THE 
  ## BIAS-CORRECTED ESTIMATES 
  #########################################
  
  conditional.like=function(mu,cee,z){
    like=dnorm(z-mu)/(pnorm(mu-cee)+pnorm(-cee-mu))
    return((abs(z)&gt;cee)*like) }
  
  conditional.like.z=function(mu,cee,z){
    return(conditional.like(mu,cee,z)*mu)
  }
  
  #########################################
  ## THE FUNCTIONS BELOW ARE USED TO OBTAIN THE 
  ## BIAS-CORRECTED CONFIDENCE INTERVAL 
  #########################################
  
  ptruncnorm.lower=function(z,mu,cee,alpha){
    A=pnorm(-cee+mu)+pnorm(-cee-mu)
    term1=pnorm(z-mu)
    term2=pnorm(-cee-mu)
    term3=pnorm(-cee-mu)+pnorm(z-mu)-pnorm(cee-mu)
    result=(1/A)*(term1*(z&lt;= -cee)+term2*(abs(z)&lt;cee)+term3*(z&gt;=cee))
    return(result-(alpha/2))
  }
  ptruncnorm.upper=function(z,mu,cee,alpha){
    A=pnorm(-cee+mu)+pnorm(-cee-mu)
    term1=pnorm(z-mu)
    term2=pnorm(-cee-mu)
    term3=pnorm(-cee-mu)+pnorm(z-mu)-pnorm(cee-mu)
    result=(1/A)*(term1*(z&lt;= -cee)+term2*(abs(z)&lt;cee)+term3*(z&gt;=cee))
    return(result-(1-alpha/2))
  }
  find.lowerz=function(mu,z,cee,alpha){
    lowerz=uniroot(ptruncnorm.lower,lower=-20,upper=20,mu=mu,cee=cee,alpha=alpha)$root
    return(lowerz-z)
  }
  
  find.upperz=function(mu,z,cee,alpha){
    upperz=uniroot(ptruncnorm.upper,lower=-20,upper=20,mu=mu,cee=cee,alpha=alpha)$root
    return(upperz-z)
  }
  
  getCI=function(z,cee,alpha){
    uppermu=uniroot(find.lowerz,interval=c(-15,15),cee=cee,z=z,alpha=alpha)$root
    lowermu=uniroot(find.upperz,interval=c(-15,15),cee=cee,z=z,alpha=alpha)$root
    out=list(lowermu,uppermu)
    names(out)=c(&quot;lowermu&quot;,&quot;uppermu&quot;)
    return(out)
  }
  source(&quot;GW-functions.R&quot;)# YOU READ IN THE FUNCTIONS FOR OUR METHOD
  cee=abs(qnorm(.5*0.1)) # Bonferroni threshold for achieving study-wide significance = 0.1
  p.value=(simAll %&gt;% select(&quot;p.value&quot;))
  p.value1 &lt;- as.numeric(unlist(p.value))
  # we use &quot;cee&quot; so R does not get confused with the function 'c'
  betahat=log(OR) # Reported OR  
  z=sign(betahat)*abs(qnorm(0.5*p.value1)) # Reported p-value = 5.7e-4, which we convert to a z-value
  
  ###################################################
  #              THE PROPOSED APPROACH              #
  ###################################################
  
  se=betahat/z # standard error of betahat
  
  mutilde1=optimize(f=conditional.like,c(-20,20),maximum=T,z=z,cee=cee)$maximum
  </code></pre>