Simulating values from an estimated value and confidence interval in R

Question

I am using the estimated county-level poverty measure from the Small Area Income and Poverty Estimates (SAIPE) as the dependent variable in a regression analysis. This value is itself the result of a model and comes complete with upper and lower 90% confidence interval bounds. To be clear, each of my 3000+ observations has an estimated value and its own confidence interval based on the sample size for that county (targeted at 2.5% of population)

I am wondering about the best way to incorporate the uncertainty in my dependent variable into my regression model.

One way I have imagined is doing a random draw from a distribution using the estimated value and confidence interval for each observation. I would then re-run my regression using this simulated value for the dependent variable in my regression analysis and compare model outcomes to the model using the estimated value. By simulating new values and comparing many times I would gain an understanding of how sensitive my findings are to the estimates on the dependent variable.

A key component of this is pulling a random number based on the estimated value and the upper and lower confidence intervals. The data is left and right censored at 0 and 100 and the estimated value is not centered within the confidence interval that is: abs(estimate-cl) != abs(estimate-cu)

I was intrigued by the discussion here: Sampling random numbers from a distribution with asymmetric confidence intervals generated by a bootstrapped estimate

but a modified version of the code just generates the estimated value Here is an example using the first 6 records from the 2008 SAIPE

sample.size<-c(1258.850,4405.300,745.900,539.725,1444.850,273.02)
POV08L90<-c(8.77,8.24,18.08,13.90,10.55,22.45)
POV08H90<-c(12.54,11.18,24.82,20.87,15.27,33.83)
POV08<-c(10.7,9.9,24.5,18.5,13.1,33.6)
test.data<-data.frame(sample.size,POV08L90,POV08H90,POV08)

gammaGenerate<-function(dat){
  for(i in 1:length(dat$sample.size)){
    n<-dat[i,"sample.size"]
    cl<-dat[i,"POV08L90"]
    cu<-dat[i,"POV08H90"]
    barx<-dat[i,"POV08"]
    talpha = qt(p=0.95,df=n-1)
    s = (cu - cl)*sqrt(n)/(2*talpha)
    kappa = 6*s*s*n*( cl - barx + talpha*s/sqrt(n) )
    gamma.shape = 4/(kappa*kappa)
    gamma.scale = s/sqrt(gamma.shape)
    gamma.shift = barx - gamma.shape*gamma.scale
    print(c(barx,(rgamma(n = 5, shape = gamma.shape) + gamma.shift)))    
  }
}
gammaGenerate(test.data)

Any help you can offer--either directing me to a better method of dealing with the uncertainty in my dependent variable, or an explanation for why my rgamma always lands at 0 would be very welcome.

Answer 1

The problem this is going to cause you in your regression is that different observations will have different uncertainties. This is known as "heteroskedasticity". Ignoring it in a regression context leads to inefficient estimates, although if the degree of heteroskedasticity is small, the inefficiency is small too. You can also get biased estimates of covariances etc. The rule of thumb I learned ("rule of thumb" = "no citation" in this response) was that if max(error variance)/min(error variance) < about 3, there wasn't likely to be much gained by attempting to correct for it.

Typically you would correct for this by using a weighted regression, e.g., weighted least squares, where, in the case of least squares, the weights are proportional to the inverse of the variances of the observations. In your case, you don't know the latter, because it's a combination of the "true" error variance and the measurement error variance. If you think the uncertainty in the dependent variable is large relative to the true error variance of your model, say as big or bigger, you could just apply weights derived from the uncertainties using the weights= option in lm and other regression functions. In fact, you could do that and compare the output to the unweighted estimates to see if there's any substantive difference.

Either way, I'd suggest using a heteroskedasticity-consistent covariance matrix estimator, for example, the one in the sandwich package. You might also find this answer to a somewhat-related question helpful if you are going to do hypothesis testing.