Estimation of population parameters on the basis of multiple samples regressions

Question

I am working with a population of approximately 1 million subjects for whom I have repeated measures of expenditure (expend). My objective is to evaluate the link between 2 variables (let's say A and B) and the amount of expenditure. Given the longitudinal nature of the data and the presence of confounding factors, I cannot afford to estimate these associations directly. I therefore undertook to carry out mixed modelling using R and the lme4 library. Here is the code for the part that interests us.

glmer(expend ~ A*B + (1|ID), family = Gamma(link = "log"), data = data)

Given the population size, it is obviously not feasible to carry out this regression on all the data. Therefore, I undertook to iteratively repeat this regression model using random samples drawn from the population with replacement.

In the end, I have a sample of 5,000 models for which I obviously have estimators of the parameters and their variance.

Using only the distribution of the estimators obtained for the fixed effects allows me, I imagine, to estimate the standard error associated with them, but I'm more interested in estimating their distribution in the initial population and therefore obtaining their standard deviation in this population.

To do this I followed the procedure below with the following notations :

$ N $ number of samples / regressions

$ n_{Ai} $ number of observations with characteristic A in sample i and $ n_{\overline{A}i} $ the number of observations without characteritic A

$ \beta_{Ai} $ estimator of the parameter A of regression i

1 : Center the parameters obtained for regression : $ \beta_{Ai}' = \beta_{Ai} - \frac{\sum{\beta_{Ai}}}{n} $

2 : Multiply by the mean squared number of observation implied in each regression : $ \beta_{Ai}'' = \beta_{Ai}' \times \frac{\sum{(\sqrt{n_{Ai}}/2}+\sqrt{n_{\overline{A}i}}/2)}{n} $

3 : Recenter the parameters : $ \beta_{Ai}''' = \beta_{Ai}'' + \frac{\sum{\beta_{Ai}}}{n} $

4 : Estimate the standard deviation of the parameter in the population on the basis of the distribution of $\beta_{Ai}'''$

To test this approach, I carried out simulations using distributions with known parameters. It turns out that the variance obtained by the procedure I have just described always slightly overestimates the variance of the source distribution.

Any help you can give me to understand this problem better would be most welcome.

Answer 1

Not quite an answer to the question you pose, but I don't think this is the right approach.

You write glmer in your question, which makes me think you're using the lme4 package. It is not really suited for very large data.

While it will take a long, long time to run, you can use parallelization and a discretization trick to strongly reduce the memory-footprint of the model.$^\dagger$

An implementation for this is available in the bam function of mgcv. This does not fit a mixed model, but a generalized additive model. However, you can use it much the same way by fitting a smooth term with the bs = "re" (random effect) option. Your model would look something like this:

library(mgcv)
GAM <- bam(y ~ A * B + s(ID, bs = "re"), family = Gamma(link = "log"), data = data,
           nthreads = 16,  discrete = TRUE)

(Change nthreads to the number of logical processors you want to use.)

Whether it is appropriate for your end goals I don't know, but another, much faster option would be to forego the GLMM and use an ordinary LMM instead on the logarithm of your outcome. This model has a different outcome, and a different set of assumptions, but this would decrease the computation time and memory-footprint considerably.

Finally, if you have a million patients, there's nothing wrong with just analyzing a subset of the total data set. Not a small subset, but halving for example would still leave you with a sample size most researchers can only dream of.

There's a great blogpost by Michael Clark on this topic that goes in much greater depth.

---

$^\dagger$: Wood, S. N., Li, Z., Shaddick, G., & Augustin, N. H. (2017). Generalized Additive Models for Gigadata: Modeling the U.K. Black Smoke Network Daily Data. Journal of the American Statistical Association, 112(519), 1199–1210. https://doi.org/10.1080/01621459.2016.1195744

Answer 2

Besides the options proposed by Frans Rodenburg (+1) you could also consider to use gee for gamma distributed data. Such model will run much faster than a random intercept model estimated with lmer. Here is an example showing how to specify such model with R's geepack.

Further, a bit in line with Frans R's answer, if you really want to stick to lmer, draw a number of "smaller" samples and show that their fixed effects and std. errors and p values hardly differ. Then report the results of one of these samples. However, I would prefer gee or one of Frans R's ideas.