When are simpler models, with fewer free parameters, better? We can get insight into this problem using resampled fit trajectories, watching the fitted values for multiple possible samples comverge as sample size increases.
We will see that models with lots of free parameters start out with very widely-spread trajectories, indicating high sample variance. These trajectories also converge more slowly (as sample size increases) than models with few free parameters. That is not a huge surprise for most people, but it is helpful to see visually just how messy “convergence” can be. Because of this messiness, a biased model that omits key parameters can, in the short term, with smaller sample sizes, be more accurate than an unbiased model.
Then we take another step and graph what happens to the fit of two parameters at once. We will visually see how correlated pairs of variables are a lot harder to fit jointly on, because the search is squeezed into a warped space, than uncorrelated pairs of variables. One solution if our sample size is small and two variables are highly correlated, is simply to omit one.
# Use this seed to exactly replicate my graphs and values below.
set.seed(2)
# Remove it to see new resampling-- and whether my overall conclusions still
# hold.
But all of this is better understood with a real example. Let's consider a data set of Chicago high school drop-out rates from 1995. (See http://rpubs.com/jeanimal/19560)
dropout <- read.table("/tmp/Dropout_Source.csv", sep = ",", header = TRUE, na.strings = "missing")
# Shorten some column names for graph readability
names(dropout)[names(dropout) == "Low_Income_Students"] <- "Low_Income"
names(dropout)[names(dropout) == "Limited_English_Students"] <- "Limited_English"
names(dropout)[names(dropout) == "Percent_white_students"] <- "Percent_White"
Below is a scatterplot of the average student math score and the school's drop-out rate. We see that higher science scores are associated with a lower drop-out rate. Oddly, lower science scores are associated with both low and high drop-out rates– a low score is not very predictive. Notice that the best schools have a drop-out rate of less than 5% while the worst schools have drop-out rates above 30%.
# install.packages('ggplot2') # run this if you have never installed it
require(ggplot2)
## Loading required package: ggplot2
ggplot(dropout, aes(Math, Dropout_Rate)) + geom_point()
We will consider the other variables in the data set soon.
Let's sample from our existing sample– called “resampling”– to generate some alternate possible histories of the world.
sampleSize <- nrow(dropout)
numTrajectories <- 5
indexPermutationMatrix <- replicate(numTrajectories, sample(1:sampleSize))
For example, the first 8 points of the first resample are highlighted in red in the graph below.
trajectory_1_first8 <- indexPermutationMatrix[1:8, 1]
ggplot(dropout, aes(Math, Dropout_Rate)) + geom_point() + geom_point(data = dropout[trajectory_1_first8,
], colour = "red", size = 3) + geom_smooth(data = dropout[trajectory_1_first8,
], method = lm, se = FALSE)
Those points can be used when fitting a linear regression, shown in blue. The slope of that line, the coefficient of Math in the regression equation, is -0.17. That means increasing Math scores have a a negative effect o drop-out rates, as we expected. (The regression shows the p-value indicates this coefficient is not significant, but we'll consider significance later.)
model_1_first8 <- lm(as.formula("Dropout_Rate~Math"), data = dropout[trajectory_1_first8,
])
coef(model_1_first8)
## (Intercept) Math
## 45.7481 -0.1707
Now consider the same analysis for trajectory 2. Its first 8 points are shown below in green. This regression gets a Math coefficient of 0.057, indicating higher math scores are associated with higher drop-out rates, which is the wrong sign!
trajectory_2_first8 <- indexPermutationMatrix[1:8, 2]
ggplot(dropout, aes(Math, Dropout_Rate)) + geom_point() + geom_point(data = dropout[trajectory_2_first8,
], colour = "green", size = 3) + geom_smooth(data = dropout[trajectory_2_first8,
], method = lm, se = FALSE)
model_2_first8 <- lm(as.formula("Dropout_Rate~Math"), data = dropout[trajectory_2_first8,
])
coef(model_2_first8)
## (Intercept) Math
## 12.0670 0.0575
This is our first inkling of sampling variance, namely the idea that different samples will result in different coefficients, and that is a form of error or noise– although it also depends on the fitting method we are using, which in this case is multiple linear regression with one predictor, Math. We will learn that different fitting methods can also increase or decrease the sampling variance.
How do we fix this? The traditional answer is to increase the sample size. A sample of 8 is fairly small for a weak effect like Math scores, especially because low math scores are not predictive, only high math scores are. So let's increase trajectory 2's sample size to 15. Reassuringly, we now get a negative coefficient for math of -0.03.
trajectory_2_first15 <- indexPermutationMatrix[1:15, 2]
ggplot(dropout, aes(Math, Dropout_Rate)) + geom_point() + geom_point(data = dropout[trajectory_2_first15,
], colour = "green", size = 3) + geom_smooth(data = dropout[trajectory_2_first15,
], method = lm, se = FALSE)
model_2_first8 <- lm(as.formula("Dropout_Rate~Math"), data = dropout[trajectory_2_first15,
])
coef(model_2_first8)
## (Intercept) Math
## 36.1130 -0.1065
Now that our intuition for resampling has been built up, we can see what the 5 different resampling of coefficients look like for a range of sample sizes and seeing how they eventually converge.
First, let's define functions to automate the fits. (I know for loops are inefficient in R. Please, experts, help me rewrite these without for loops.)
# define fitter functions if not defined already
fitter <- function(formula, dataToFit, trajIndex, indexPermutation, startSampleSize,
endSampleSize) {
outputFrame <- data.frame()
for (sampleSize in startSampleSize:endSampleSize) {
sampledRowIndices <- indexPermutation[1:sampleSize]
model <- lm(formula, dataToFit[sampledRowIndices, ])
outputFrame <- rbind(outputFrame, cbind(trajectory = trajIndex, sampleSize = sampleSize,
t(coef(model))))
}
return(outputFrame)
}
manyFitter <- function(formula, dataToFit, numTraj, indexPermutationMatrix,
minSampleSize, maxSampleSize) {
outputFrame <- data.frame()
for (trajIndex in 1:numTraj) {
tempOut <- fitter(formula, dataToFit, trajIndex, indexPermutationMatrix[,
trajIndex], minSampleSize, maxSampleSize)
outputFrame <- rbind(outputFrame, tempOut)
}
return(outputFrame)
}
Now let's generate regression coefficients for each trajectory at each sample size and graph the results.
schOut_math <- manyFitter(as.formula("Dropout_Rate~Math"), dropout, numTrajectories,
indexPermutationMatrix, 8, 60)
g_math <- ggplot(schOut_math, aes(sampleSize, Math, color = as.factor(trajectory))) +
geom_line() + ylab("Math slope") + labs(color = "trajectory")
g_math
The start of the graph is at sample size 8. The red line is trajectory 1, beginning with a coefficient of -0.17. As the sample size is increased, the fits from trajectory 1 increase but overshoot and then go too low but eventually converge toward a number near -0.1.
But the real use of this graph is it gives an idea of the “spread” of the points, which is an indicator of sampling variance. At sample size 8, the spread is huge. At sample size 20, it's much better but still quite wide. It's not until sample size good thing is they converge by sample size gets to about 50 that the lines start to agree on a coefficient near -0.1.
In the final version of analyzing drop-out rates, we will of course use the whole sample. So why are we taking sub-samples? We are answering a meta-level question: What is the right way to analyze the data? What are the relative merits of a model with one free parameter, Math, relative to a model with many free parameters. And we now have the tools to answer that.
So far, we made fits of DropOut rates on Math. But there are many other variables in the data set. What if we did a regression on more of them? Here's a quick look using Math scores, Science scores, Reading scores, percent of students with limited English proficiency, percent of low income students, and percent of white students.
schOut_lots <- manyFitter(as.formula("Dropout_Rate~Math+Science+Reading+Limited_English+Low_Income+Percent_White"),
dropout, numTrajectories, indexPermutationMatrix, 8, 60)
g_lots <- ggplot(schOut_lots, aes(sampleSize, Math, color = as.factor(trajectory))) +
geom_line() + ylab("Math slope") + labs(color = "trajectory")
g_lots
At first, the graph might not look so bad. But look at the scales on the y axis! There is a much wider range of coefficients for Math scores than before.
Let's redo the graphs together the same y-axis. This time let's color all the trajectories from the same model as one color to focus on comparing across models. All five trajectories from the bigger model, called “lots”, are in red. Those from the model with just math are in a light blue.
schOut_both <- rbind(cbind(schOut_lots[, 1:4], frmla = "lots"), cbind(schOut_math,
frmla = "Math")) # combine the two dataframes into one
g_math_and_lots <- ggplot(schOut_both, aes(sampleSize, Math, group = interaction(frmla,
trajectory), color = as.factor(frmla))) + geom_line() + ylab("Math slope") +
labs(color = "frmla")
g_math_and_lots
Clearly, fitting on more variables greatly increased the sampling variance for estimating the impact of Math scores on drop-out rates.
Lesson: If your sample size was only 40, the “lots” model would NOT have have finished converging yet. You would have risked a garbage estimate for the impact of Math, including positive values, claiming Math increases drop-out rates!
Let's zoom in a bit because we can barely see the fits from the formula with just Math.
g_math_and_lots + coord_cartesian(ylim = c(-0.5, 0.5))
Now you can also see that the “lots” model at sample size 60 is converging to a smaller-sized negative slope of Math (closer to zero slope) than the model fit with just Math and no other varables. The lots model says that, taking other variables into account (e.g. low income students and limited English students), math scores do less to reduce dropout rates than claimed by the model fit just to Math.
Which is right? The model with just Math scores could be biased– that is the risk of a small model. It could be leaving out relevant, important variables. Or it could be that the “lots” model is being fooled by all those other variables.