There won't be any advantage to bootstrapping the SE in your example because you have a very large sample size. The distribution of means of that sample size is going to be normal, not skewed, because of the central limit theorem [CLT] (try hist(skewLeftbootData)). Even if it were skewed the SE is going to be so small because of N that the SE is not going to be appreciably skewed anyway. Then you're using for proof the backward calculation of an SD based on the SE of the bootstrap distribution calculated through conventional means. Even if the bootstrap distribution were skewed you've just tossed out one of the reasons you might do bootstrap in this case.
Bootstrapping would be more compelling if you had substantially smaller sample (say 12) and calculated your SE as the middle 67% of the bootstrapped data by cutoffs of the sorted bootstrap distribution. Then you would see that that is a different estimate than an SE calculated from the conventional SD. You also wouldn't then calculate a bootstrapped SD based on the cut offs.
I think I have managed to rewrite the Bootstrap 0.632+ from R to Java using the Weka Java API. The original R function can be found in the bootpred method inside the bootstrap package (link).
As you can see from the source code, I have used the corrected $\hat{R'}$ and $\hat{Err}^{(1)'}$ with the final (corrected) equation (32) from the original article.
However, despite correcting for abnormal values, I sometimes still get negative error rate, which is of course impossible and therefore invalid. Also, I have noticed that the difference between the 0.632 and the 0.632+ method is minimal, if any.
If someone finds any errors in my source code, I would be really grateful if you could point them out.
public class Bootstrap632plus extends AbstractPerformance {
private final int repeats;
private double Err632;
private double resub;
public Bootstrap632plus(Instances instances, int repeats) {
super(instances);
this.repeats = repeats;
}
public double getErr632() {
return Err632;
}
public double getResub() {
return resub;
}
@Override
public double getErrorRate(final MachineLearningAlgorithm machineLearningAlgorithm, final Random seed) throws Exception {
// First component
double err = predictionError(machineLearningAlgorithm);
this.resub = err;
// Error rates
List<Double> errorRates = Collections.synchronizedList(new ArrayList<>());
// GAMA related stuff
final int numClasses = instances.numClasses();
AtomicIntegerArray p_l = new AtomicIntegerArray(numClasses);
AtomicIntegerArray q_l = new AtomicIntegerArray(numClasses);
// Bootstrap iterations
seed.ints(repeats).parallel().forEach(randomSeed -> {
// Get error rate
Evaluation evaluation = bootstrapIteration(machineLearningAlgorithm, randomSeed);
errorRates.add(evaluation.errorRate());
/*
GAMA VARIABLE
Confusion matrix:
- first dimension (rows): real distribution for first class
- second dimension (columns): predicted distribution for first class
p_l = observed proportions of responses where y_i equals l
- sum by l-th row (first dimension)
q_l = observer proportions of predicted responses where y_i equals l
- sum by l-th column (second dimension)
GAMA = SUM_by_l(p_l * (1 - q_l))
*/
double[][] confusionMatrix = evaluation.confusionMatrix();
for(int l = 0; l < numClasses; l++) {
int p_tmp = 0, q_tmp = 0;
for(int n = 0; n < numClasses; n++) {
// Sum for l-th class
p_tmp += confusionMatrix[l][n];
q_tmp += confusionMatrix[n][l];
}
// Add data for l-th class
p_l.addAndGet(l, p_tmp);
q_l.addAndGet(l, q_tmp);
}
});
// Second component
double Err1 = errorRates.stream().mapToDouble(i -> i).average().orElse(0);
// Plain 0.632 bootstrap
Err632 = .368*err + .632*Err1;
// GAMA
final double observations = instances.size() * repeats;
double gama = 0;
for(int l = 0; l < numClasses; l++) {
// Normalize numbers -> divide by number of all observations (repeats * dataset size)
gama += ((double)p_l.get(l) / observations) * (1 - ((double)q_l.get(l) / observations));
}
// Relative overfitting rate (R)
double R = (Err1 - err) / (gama - err);
// Modified variables (according to original journal article)
double Err1_ = Double.min(Err1, gama);
double R_ = R;
// R can fall out of [0, 1] -> set it to 0
if(!(Err1 > err && gama > err)) {
R_ = 0;
}
// The 0.632+ bootstrap (as used in original article)
double Err632plus = Err632 + (Err1_ - err) * (.368 * .632 * R_) / (1 - .368 * R_);
return Err632plus;
}
/**
* Prediction error: first component of the 0.632+ bootstrap.
* Train the classifier on the whole dataset and then also test it on the whole dataset.
*
* @param machineLearningAlgorithm Specified machine learning algorithm
* @return prediction error [0, 1]
* @throws Exception exception
*/
private double predictionError(final MachineLearningAlgorithm machineLearningAlgorithm) throws Exception {
// Train
Classifier classifier = ClassifierFactory.instantiate(machineLearningAlgorithm);
classifier.buildClassifier(instances);
// Test
Evaluation evaluation = new Evaluation(instances);
evaluation.evaluateModel(classifier, instances);
// Return error rate
return evaluation.errorRate();
}
/**
* One iteration of the Leave-one-out Bootstrap Cross-Validation.
* @return
* @throws Exception
*/
private Evaluation bootstrapIteration(final MachineLearningAlgorithm machineLearningAlgorithm, final int randomSeed) {
try {
final int SIZE = instances.size();
final Random r = new Random(randomSeed);
// Custom sampling (100%, with replacement)
List<Instance> TRAIN = new ArrayList<>(SIZE); // Empty list (add one-by-one)
List<Instance> TEST = new ArrayList<>(instances); // Full (remove one-by-one)
for(int i = 0; i < SIZE; i++) {
// Random select instance
Instance instance = instances.get(r.nextInt(SIZE));
// Add to TRAIN, remove from TEST
TRAIN.add(instance);
TEST.remove(instance);
}
// Train
Instances trainSet = new Instances(instances, TRAIN.size());
trainSet.addAll(TRAIN);
Classifier classifier = ClassifierFactory.instantiate(machineLearningAlgorithm);
classifier.buildClassifier(trainSet);
// Test set
Instances testSet = new Instances(instances, TEST.size());
testSet.addAll(TEST);
// Test
Evaluation evaluation = new Evaluation(instances);
evaluation.evaluateModel(classifier, testSet);
// Return the evaluation (for further processing)
return evaluation;
} catch(Exception e) {
throw new RuntimeException(e);
}
}
}
Best Answer
First, you need to write a function that calculates the "percent attenuation" for a data set and a "Risk Factor." It would: (a) fit both "Model 1" and "Model 1 + Risk Factor" to the data set; (b) extract the $\beta_{SES}$ regression coefficient for
SESfrom both models; (c) calculate the "percent attenuation" following your equation.Then you take multiple bootstrap samples of the original data set, apply that function to each of them, and keep track of all the corresponding "percent attenuation" values.
Finally, for a simple estimate of the 95% CI, put the "percent attenuation" values in order, and select the 2.5th and 97.5th percentiles. For example, if you did 999 bootstrap samples you could take the 25th and 975th from the ordered list.
This is straightforward in R, using tools in the
bootpackage. Theboot.ci()function in that package provides for more refined CI estimates that might be helpful if there is bias or skew in your "percent attenuation" measure. I don't use Stata. (Software-specific programming details are off-topic on this site anyway unless there is associated statistical content.)There is at least one detail with respect to the cited paper that you might need to consider. The paper used multiple imputation to deal with some types of missing variable values. If that's the case, then the $\beta$ values should not have been from a single regression but should have been averages among separate fits on each of the imputed data sets. See Stef van Buuren's multiple imputation book. If you want to repeat what was done in that paper (and the paper handled the multiple imputation correctly), then your function to calculate "percent attenuation" should do the same.
A couple of warnings
First, that measure of "percent attenuation" is not additive across multiple "Risk Factors." You can see that in Tables 3 and 4 of the paper. The attenuation when all risk factors are combined is less than the sum of the attenuations associated with each of the individual risk factors. That's to be expected when risk factors are correlated among each other and with
SES. So be very careful when interpreting those values.Second, it is hard to draw causal inferences from a study like this. If some of the "risk factors" are themselves caused by
SESand then mediate the effect ofSESon mortality, then the "attenuation" estimate derived by adding them to the model might be an "overadjustment for mediators" as Hernán and Robins explain in Section 18.2 of What If?. IfSESis driving those risk factors, those "attenuation" values don't necessarily mean thatSESis correspondingly less important. They would then be related to how much of theSESeffect the risk factors mediate.