Given sample size, group means, and misc. post-hoc range statistics, can you suggest a good way to estimate variance through simulation? I previously asked this question about the validity of my solutions for for $SE$ given $n$, $\bar{X_i}$ and summary statistics from post-hoc multiple comparisons such as Fisher's $LSD$ and Tukey's $HSD$, but I would like a general approach that can be applied to other post-hoc tests as well as other statistics. My hope is that this approach would be less prone to arithmetic error and easier to implement than my trying to find an analytical solution each time I encounter a new statistic.
I would like to get some estimate of $SE$ from previous studies that I will use in a meta-analysis. In addition, I would like to be able to do this for any statistic.
Is there a general simulation / optimization approach to this problem?


*

*$\bar{X_1}$, $\bar{X_2}$, some statistic, and n are given

*I am looking for $SE\pm 10\%$

*As in the earlier post, over estimates of SE are o.k.


My idea is that the solution would be something like the following examples in which LSD is given:
se.est <- function(x1bar, x2bar, lsd_obs, n, se) {
       y <- c(rnorm(n, x1bar, se), rnorm(n, x2bar, se)).
       x <- c(rep(1,n), rep(2,n))
       mse <- sum(lm(y~factor(x))$residuals^2)/(2*n-2)
       lsd_est <- qt(0.975,n)*sqrt(2*mse/2)
       ans <- (lsd_obs - lsd_est)^2
       }

SE <- optimize(par = c(x1bar, x2bar, lsd_obs,n), fn = se.est)

or
for (i in 1:10000){
     se[i] <- rgamma(1, 1, 0.01)
     y[i] <- c(rnorm(n, x1bar, se), rnorm(n, x2bar, se)).
     x[i] <- c(rep(1,n), rep(2,n))
     mse[i] <- sum(lm(y~factor(x))$residuals^2)/(2*n-2)
     lsd_est[i] <- qt(0.975,n)*sqrt(2*mse/2)
     }

 se <- se[which.min((lsd_est - lsd_obs)^2)]

Your suggestions for a quick and effective approach would be appreciated; also, is one of the above preferable or substantially more efficient? if you think that an analytical solution is preferable, please see my previous post.
 A: With a small modification, the first version will be far more efficient. As it's written, se.est is a random function. Even at the same arguments, its value will change each time because of rnorm. This will mess up optimize. You should use the same random numbers each time se.est is called. Here's one way:
e <- c(rnorm(2*n,0,1))
se.est <- function(se,x1bar, x2bar, lsd_obs, n, e) {
       y <- se*e + c(rep(x1bar,n),rep(x2bar,n)
       x <- c(rep(1,n), rep(2,n))
       mse <- sum(lm(y~factor(x))$residuals^2)/(2*n-2)
       lsd_est <- qt(0.975,n)*sqrt(2*mse/2)
       ans <- (lsd_obs - lsd_est)^2
       }

SE <- optimize(f=se.est,interval=c(0.1,100),x1bar,x2bar,lsd_obs,n,e)

I didn't try to run this, so I may have made typos or syntax errors, but hopefully the idea is clear. If n is small, you should probably optimize multiple times with different e's, and use the something like the average or maybe the 90th percentile of the resulting SE's.
This may not be a concern, but if the original data was not normal, this method could give inaccurate results. Without assuming normality, you could write your problem as:
 max_y se(y) subject to mean1(y)=x1bar, mean(2)=x2bar, statistic(y)=s

where you allow y to be a sample of size n with any distribution. Without any more restrictions, the solution will often be infinite (it depends on the statistic). 
