Can I improve linear model coefficient estimates using group information without working it into model? I am fitting a linear model in order to predict future observations. The training data consists of about 1000 observations. Each observation comes from one of 10 individuals, which means I have about 100 observations per individual. The future data to which the model will be applied will come from different individuals than those in the training set. The individuals are expected to vary systematically (they could have different slopes, intercepts, etc).
I have an understanding of how one would usually account for this individual-level variation in a linear model. Dummy variables, interaction terms, partial pooling, etc. However, these options only work for prediction if I have the same individuals in the prediction data as I do in the training data. This is not the case in my situation; the individuals at prediction time are not the same individuals at training time, so the model cannot contain any terms which represent the individual.
My question: if my sole concern is in getting the best coefficient estimates possible, is there anything more I can do beyond fitting the linear model as usual with no accounting for the individuals?
I came across some info about clustered standard errors, but the entire purpose of that seems to be in correcting for underestimation of error-term variance. It doesn't seem like introducing clustered standard errors will at all impact the coefficient estimates themselves (what I really care about), though please correct me if I'm wrong.
Some relevant points:

*

*Thankfully, there is strong reason to believe that the individuals in the training set are a representative sample of individuals that will be seen in the wild.

*I have no information that could account for individual-level variation. In the training data, I only have indicators for individuals. There is no way to model the individuals' slopes/intercepts as depending on some underlying property of the individuals.

Visual explanation:
library(ggplot2)

set.seed(123)
n_individuals <- 5
obs_per_individual <- 100

df_params <- data.frame(individual = factor(1:n_individuals),
                        slope = rnorm(n = n_individuals, mean = 1, sd = 0.2),
                        intercept = rnorm(n = n_individuals, mean = 0, sd = 0.2))

df_list <- list()
for (i in 1:n_individuals) {
  df_iter <- data.frame(individual = rep(df_params$individual[i], obs_per_individual),
                    x = runif(n = obs_per_individual))
  df_iter$y <- df_params$intercept[i] + df_iter$x * df_params$slope[i] + rnorm(n = obs_per_individual, mean = 0, sd = 0.1)
  df_list[[i]] <- df_iter
}
df <- bind_rows(df_list)

ggplot(df, aes(x = x, y = y)) +
  geom_point(aes(color = individual), size = 0.5) +
  geom_smooth(aes(color = individual), method = "lm", fullrange = TRUE, se = FALSE, size = 0.5) +
  geom_smooth(mapping = aes(x = x, y = y), method = "lm", fullrange = TRUE, se = FALSE, linetype = "dashed") +
  scale_color_discrete(guide = NULL)




The dashed line is the overall regression line. Colors represent individuals. Is there any way to improve upon the regression in a way that does not rely on knowing the individuals in advance at prediction time?

When estimating a mixed effects model we have to distinguish between marginal (sometimes called "unconditional") and conditional estimates. Conditioning refers to whether these estimates are subject-specific or not. In a way, all models other than conditional GLMMs produce some variation of a marginal estimators. Please note that while I mention "subjects" here, this extends to an arbitrary level on nesting.
The practical implications of this distinction is that when we are making estimates about new "unseen" individuals, e.g. subjects $s_{n+1}$ and $s_{n+2}$, we can only predict their marginal estimates. That is because we assume that any difference between them and our training sample will be due to random variation which is by definition "0" in expectation. (For example in a single-level mixed effects model with $y_{i,j} = X_{ij}\beta + Z_{ij}u_i + \epsilon_{ij}$, our working assumption for the random effect $u_i$ is that $u_i \sim N(0,\sigma^2 D)$.) As such having conditional estimates for $s_{n+1}$ and $s_{n+2}$ is inapplicable; if we assume there are any systemic differences between these new subjects and our training sample, those differences will be realised in the fixed effects ("explanatory variables") rather than the random effects ("nuisance parameters") of our model.
Two relatively mature reference on the matter are: A Comparison of Cluster-Specific and Population-Averaged Approaches for Analyzing Correlated Binary Data (1991) by Neuhaus et al. and Marginalized multilevel models and likelihood inference (2000) by Heagerty and Zeger; CV.SE has a related thread on Difference between marginal and conditional treatment effect? Relating to regression vs. propensity score methods that will be handy to check too.  West et al.'s "Linear Mixed Models: A Practical Guide Using Statistical Software" (see Sect. 2.3.2 on the "The Marginal Model Implied by an LMM" particularly) is also a good and easy-to-follow resource on the matter.
