Confidence that Optimum Lies Within a Given Region I have an ordinary least squares regression equation. Through calculus or simulation, I can find the combination of explanatory values that maximize the estimated mean response; I can also establish a confidence interval for the estimated mean response. How would I go about determining (with a certain level of confidence) the region that contains the true maximum response?
 A: How many variables do you have?  If only two (or if you are willing to hold the rest constant), you could do it by brute-force simulation, exploiting the multivariate-normality of your estimates (or bootstrapping your estimates).  Bascially what you'd do would be to define a grid of your data (probably using expand.grid).  If your model is m, you'd then simulate $M$ replicates of $\beta$ by b = mvrnorm(1000,m$coef,vcov(m)).  You could then get fitted values across your grid X by XB = X%*%t(b).  You could see which maximize $y$ via apply(XB,2,which.max).  Then coerce your long-format grid into a wide format matrix, in which each entry is the count of the number of times that combination of x1 and x2 maximized the function based on the $M$ simulates.  You could them calculate a 2-dimensional kernel density across that matrix, and draw contour lines around the region that contains 95% of the density.  
I tried coding up an example with fake data, but I had trouble making a fake dataset where an OLS model gave a sufficiently uncertain maximum. In >2-dimensional space, defining the maximum would be trickier using brute force, and harder to interpret anyway.  In my experience, functionals of (multivariate) normally-distrbuted parameters aren't generally multivariate-normally distributed (though I haven't seen an analytical proof of this), so I don't know whether it'd be realistic to get a vector of means and a covariance matrix for the $X$ values which maximize $Y$.
