Bootstrapping regression

When the linear regression assumption is violated, options are limited to generalize the findings:

1. If the residuals show problems with heteroscedasticity or non-normality you can try transforming the raw data or try weighted regression.
2. If you have a linear assumption violation, you can try logistic regression.
3. You can also try a robust bootstrapping regression.

Bootstrapping regression is very similar to bootstrapping correlation.

# function for boot
bootReg <- function(formula, data, index) {
   d <- data[index,]
   model <- lm(formula, data = d)
   return (coef(model))
}
# create bootReg object
bootRegObj <- boot(statistics = bootReg, formula = y ~ x1 + x2 + x3, data = data, R=2000)
# get confidence interval for each predictor
boot.ci(bootRegObj, type = "bca", index = 1)
boot.ci(bootRegObj, type = "bca", index = 2)
boot.ci(bootRegObj, type = "bca", index = 3)

Bootstrapping Correlation

Bootstrapping comes in handy when there is doubt that the usual distributional assumptions and asymptotic results are valid and accurate. We can use it to compute estimated standard errors, confidence intervals and hypothesis testing.

The basic idea of Bootstrapping is as the following:

1. You treat your sample as population and repeatedly draw new samples from it with replacement. All original observations have equal probability of being drawn into the new sample.
2. This is repeated n times. In each iteration, some observations from your original sample are drawn multiple times while some observations may not be drawn at all. After n iterations, you have n stored bootstrap estimates of the statistic(s) of interest (e.g., mean, or correlation coefficient).
3. summary statistics such as the mean, median and the standard deviation of the n bootstrap-estimates are calculated.

The following is an example of application of bootstrapping in correlation computation.

Suppose our data violates the assumption of normal distribution, and we decide to apply bootstrapping analysis of kendall’s tau on this set of data.

library(boot)
# function for boot object
bootTau <- function(data, indice) {
     cor(data$NameOfColumn2[indice], data$NameOfColumn3[indice], use="complete.obs", method="kendall")
}
# create a boot object. Parameters: data, function, repeat times
boot_kendall <- boot(data, bootTau, 2000)
# display summary of the created boot object
boot_kendall
# get 95% confidence interval for the boot object
boot.ci(boot_kendall, conf= 0.95)

Kendall’s tau

Kendall’s tau is another non-parametric correlation, and it hsould be used rather than Spearman’s coefficient when you have a small data set with a large number of tied ranks.

# correlation
cor(data, method="kendall")
# confidence interval
cor.test(data, method="kendall")

Spearman’s Correlation

When you want to test correlation between variables, but find out that the data violates parametric assumption such as non-normally distributed data, you can try Spearman’s correlation.

The procedure of carrying out Spearman’s correlation is almost the same to Pearson’s correlation, but you do need to specify the functions’ parameters.

# correlation
cor(data, method="Spearman")
# significance test
matrixData <- as.matrix(data)
rcorr(matrixData)
# confidence interval
cor.test(data, method="Spearman")

Person’s Correlation

Assumption of Pearson’s r:

• Data for the two variables are interval.
• The distribution of two variables are normally distributed; there could be an exception that one variable is categorical, but with only two categories.

Pearson’s Correlation in R:

If you want to get the Person’s correlation coefficients, cor() function will be enough.

# correlation matrix
cor(data)
# correlation matrix
cor(data[, c("Name of Column 1", "Name of Column 2", "Name of Column 3")])

Since correlation coefficients are effect sizes, so we can interpret these values without really needing to worry about p-values. However, if you want to get the p-values of each coefficient, you will need Hmisc package, and then you can use rcorr() function.

library(Hmisc)
# rcorr() only works on matrix
matrixData <- as.matrix(data[, c("Name of Column 1", "Name of Column 2", "Name of Column 3")])
# rcorr(matrixData)

If you are interested in looking at the confidence interval of correlation coefficient, cor.test() would be the function serving the purpose.

# cor.test() can only work on one pair of variables
cor.test(data$Column1, data$Column2)

ggplot and regression

This post serves as a sample of how to use ggplot to plot key graphs for regression.

The model that is produced by lm() is a type of data set, which has variables in it. One of those variables is predicted values (or called fitted values). The fitted value can be accessed by modelNmae$fitted.values; it is a good idea to save these fitted values in the original data frame:

data$fitted <- model$fitted.values

Adjusted predicted value:

The computer calculates a new model without a particular case and then uses this new model to predict the value of the outcome variable for the case that was excluded. If a case does not exert a large influence over the model, then we would expect the adjusted predicted value to be similar to the predicted value when the case is included. Put it similarly, if a model is stable then the predicted value of a case should be the same regardless of whether or not that case was used to calculate the model.

Studentized residual:

Studentized residual = (The adjusted predicted value – the original observed value)/standard error

This residual can be compared across different regression analyses because it is measured in standard unit. It is very useful to assess the influence of a case on the ability of the model to predict that case. However, how a case influences the model as a whole can not be told by Studentized residual. Cook’s distance measures the overall influence of a case on the model.

Suppose we have studentized residuals saved in our original data frame, then we can plot a histogram of the studentized residuals.

# histogram
histogram <- ggplot(data, aes(studentized.residuals)) + opts(legend.position = "none") + geom_histogram(aes(y=..density..), color = "black", fill = "white") + labs(x="Studentized Residual", y="Density")
# density curve
histogram + stat_function(fun = dnorm, args = list(mean = mean(data$studentized.residuals), na.rm = TRUE), sd = sd(data$studentized.residuals, na.rm = TRUE)), color = "red", size = 1)

We can create a Q-Q plot.

qqplot.resid <- qplot(sample = data$studentized.residuals, stat = "qq") + labs(x = "Theoretical values", y="Observed values")

We can also plot a scatterplot of studentized residuals against predicted values:

scatter <- ggplot(data, aes(fitted, studentized.residuals))
scatter + geom_point() + geom_smooth(method="lm", color="Blue") + labs(x="Fitted value", y="Studentized Residual")