Estimate linear regression coefficients and standard errors using sub-samples of dataset I am working with a very large dataset (250 million records) and I want to use linear regression to estimate how some variables are related to the outcome variable.  I have a lot of categorical variables that I am controlling for using one-hot encoding, so when I create the model matrix the dimensions are about 250,000,000 by 4,000.
After trying a number of different approaches, I think it will be too difficult to run the regression using the entire dataset, even on a large cloud computer.  So now I am interested in estimating my parameters of interest using sub-samples.
Using a toy dataset, it seems that getting many subsamples of the data and taking the mean coefficient estimates yields a value very close to the true mean.  However, the estimates of standard errors obviously depends a lot on the sample size, and so I am not sure how to go about estimating this parameter.
For example, with the dataset diamonds from the ggplot2 package in R:
library(tidyverse)

set.seed(100)

n_sample <- 1000

truemod <- lm(price ~ depth, data=diamonds)

#With subsamples of 90% of the full dataset
val90 <- NULL
se90 <- NULL
for (i in 1:n_sample){
  m <- lm(price ~ depth, data=diamonds[sample(1:nrow(diamonds), nrow(diamonds)*0.9), ])
  val90 <- c(val90, m$coefficients['depth'])
  se90 <- c(se90, summary(m)$coefficients['depth', 'Std. Error'])
}

#With subsamples of 50% of the full dataset
val50 <- NULL
se50 <- NULL
for (i in 1:n_sample){
  m <- lm(price ~ depth, data=diamonds[sample(1:nrow(diamonds), nrow(diamonds)*0.5), ])
  val50 <- c(val50, m$coefficients['depth'])
  se50 <- c(se50, summary(m)$coefficients['depth', 'Std. Error'])
}

#With subsamples of 10% of the full dataset
val10 <- NULL
se10 <- NULL
for (i in 1:n_sample){
  m <- lm(price ~ depth, data=diamonds[sample(1:nrow(diamonds), nrow(diamonds)*0.1), ])
  val10 <- c(val10, m$coefficients['depth'])
  se10 <- c(se10, summary(m)$coefficients['depth', 'Std. Error'])
}

plt <- bind_rows(data.frame(value=c(val90, se90), sample_percent='90',
                            parameter = rep(c('coefficient', 'standard error'), each=n_sample)),
                 data.frame(value=c(val50, se50), sample_percent='50',
                            parameter = rep(c('coefficient', 'standard error'), each=n_sample)),
                 data.frame(value=c(val10, se10), sample_percent='10',
                            parameter = rep(c('coefficient', 'standard error'), each=n_sample)))

true <- data.frame(value=c(truemod$coefficients['depth'], 
                       summary(truemod)$coefficients['depth', 'Std. Error']),
                   parameter = c('coefficient', 'standard error'))
                   
ggplot() + 
  geom_density(data=plt, aes(x=value, fill=sample_percent), alpha=0.5) + 
  geom_vline(data=true, aes(xintercept=value, color=parameter)) + 
  facet_wrap(parameter ~ ., scales="free")
ggsave('~/example.png')




So my main question is: how can I use sub-sampling to get an estimate of the true standard errors?  By "true standard errors" I mean the standard errors I would get if I ran a linear regression with the entire dataset.
I'm also unsure of whether it would make more sense to use sampling with replacement or without replacement, as well as whether this approach could be appropriately called a bootstrap approach.
 A: The variance of an ordinary least squares estimate (that is, the square of its standard error) is inversely proportional to the sample size.
In practice, there is an additional form of uncertainty in this general rule: the standard error itself has to be estimated because the variance of the error terms in the model is unknown.  However, provided you obtain a sufficiently large sample to estimate this error variance reasonably accurately, you may safely extrapolate.
Consider, then, fitting the model with one very small subsample of around 20,000 observations (randomly selected, of course, but in such a way that all values of all categories occur in the subsample) to see what kinds of standard errors you are getting.  (This number is about the minimum needed to do a regression with 4,000 variables.)  Do that again with a larger subsample, perhaps around the 100,000 shown here, to check that the standard errors tend to decrease by about a factor of $\sqrt{5}$ and to improve your estimates of them.  If these standard errors aren't good enough for your purposes, extrapolate from there.  E.g., if you want them to be half of what you're seeing with 100K observations, you will need to quadruple the subsample size.  Time your work: the computation should scale almost linearly with the subsample size, giving you a good basis to anticipate how long you will have to wait to get results for the larger subsamples.

To illustrate, I estimated coefficients for a categorical variable with 4000 levels using subsamples of sizes 8000 through 100,000 taken from a dataset of a quarter million observations.  (It took a few minutes ;-).)  The inverse proportion rule implies that on a log-log plot, the squares of their standard errors ought to be scattered around a line of slope $-1.$

All $4000$ estimation variances are plotted (distinguished by color) for each of twelve sample sizes.  The gray line has a slope of $-1.$
You can see the trend follows what theory tells us to expect. For this example, may therefore feel comfortable extrapolating it to a sample size of 250,000,000, where the typical estimation variance ought to be close to $0.000017,$ corresponding to a standard error of $0.004.$
As another example, I performed a series of fits with a 40-category variable (almost instantaneously) with subsample sizes from 80 out to a million:

The slope of $-1$ as well as the regular, steady shrinking of the range of standard errors are characteristic.
A: A possible frame to answer your question could be based on the sampling theory linear regression estimator of the finite population variance. You will first, however, need an estimate for N of the superpopulation from which you have amassed 250 million records.
Then apply the usual sampling theory regression estimates and associated variance (which is subject to a finite sampling correction factor, f = (N-n)/N, see reference here, to further assess the variance associated with subpopulation sampling).
Here is also a related work The Finite-Population Linear Regression Estimator and Estimators of its Variance—An Empirical Study, to quote:

The usual variance estimator for the linear regression estimator of a finite population total is examined under some prediction (superpopulation) models. Its bias is compared with that of the least squares variance estimator. Also described are three bias-robust alternatives, one of which is the jackknife variance estimator. The theoretical results are supported by an empirical study in which simple and restricted random samples as well as some purposive (nonrandom) samples are drawn from six real populations.

