You're on the right track, but usually this problem is interpreted as a form of dependence. Children's reading abilities will most likely be affected by their school, and so measurements from a particular school are not quite independent from each other.
Commonly for this kind of problem, you assume that the schools come from some larger population of schools and estimate the variance between them. This can be done in a mixed model, using a random effect for school. A random effect is like an extra error term, except that its variance is estimated between groups of observations (not individual observations). In essence you add an offset from the intercept for every school.
A simple example in R:
LMM <- lmer(reading_ability ~ siblings + age + books + (1 | school)) # random intercept
If you believe the effect of other variables is also dependent on the school (e.g., some schools have better access to books), you can estimate a random slope for that variable as well.
Continuing the example in R:
LMM <- lmer(reading_ability ~ siblings + age + (books | school)) # random intercept & slope
Lastly, schools are undoubtedly nested within cities. Cities may have different financial support for schools and so the entire reasoning used for schools can also be used for cities. You can model this hierarchically, either explicitly (using Bayesian hierarchical models) or using a nested random effect:
LMM <- lmer(reading_ability ~ siblings + age + (1 | city/school)) # nested effect
LMM <- lmer(reading_ability ~ siblings + age + (books | city/school)) # + slope for books
Actually your example is quite close to what most introductory text books on mixed models use as an example, so allow me to complete it: If you know the classes the children are in, you can make the same argument for nesting. Updating the example, your model would then look something like this:
LMM <- lmer(reading_ability ~ siblings + age + (1 | city/school/class)) # intercept only
LMM <- lmer(reading_ability ~ siblings + age + (both | city/school/class)) # intercept & slope
Have a look at other questions tagged mixed models, I've added the tag to your question!