1What the method is for

When observations are nested within groups, what breaks in a regression that assumes independence?

In the Maibara–Sekigahara survey the response is the rate of vocabulary loss in each settlement, and settlements belong to the catchment area of a lower-secondary school district. Settlements in the same district share a commuting zone; the same cohort exchanges language at the same school. The premise that settlements within a district are mutually independent observations therefore does not hold as it stands.

What happens if a regression assuming independence is applied anyway? In the framework of Perspective 2, the variability of an estimate is set by the inverse of the information \(I(\theta)\), and the confidence interval follows from it. When observations carry overlapping information but are counted as independent, there appears to be more information than there is, and \(I(\theta)\) is overstated. The consequence is that standard errors come out too small, confidence intervals too narrow, and tests too ready to declare significance. The frequentist promise — the hit rate of a procedure under repetition — is broken at that point.

There are two remedies. One is to estimate a separate intercept for each district as a fixed effect; the other is to represent the spread of the district intercepts with a single variance parameter. The linear mixed model is the second of these.

2What a random effect puts into the model

Not a separate value for each group, but the width of the spread of those values.

WRITTEN AS A FIXED EFFECT
\[ y_i \;=\; \beta_0 \;+\; \sum_{d=2}^{D} \gamma_d\,\mathbf{1}\!\left[d(i)=d\right] \;+\; \mathbf{x}_i^{\top}\boldsymbol{\beta} \;+\; \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}\!\left(0,\sigma_\varepsilon^2\right) \]

\(\gamma_d\) … the intercept for district \(d\). Excluding the reference district, \(D-1\) parameters are added
\(d(i)\) … the district to which settlement \(i\) belongs. Here \(D=5\)

WRITTEN AS A RANDOM EFFECT
\[ y_i \;=\; \beta_0 \;+\; u_{d(i)} \;+\; \mathbf{x}_i^{\top}\boldsymbol{\beta} \;+\; \varepsilon_i, \qquad u_d \sim \mathcal{N}\!\left(0,\sigma_u^2\right),\quad \varepsilon_i \sim \mathcal{N}\!\left(0,\sigma_\varepsilon^2\right) \]

\(u_d\) … the departure of district \(d\)'s intercept. It is not itself a parameter but a random variable drawn from a normal distribution
\(\sigma_u^2\) … this one quantity is what gets estimated. It does not grow whether there are five districts or fifty

The difference lies in what is estimated. The fixed-effect version adds a parameter per district, and each district is determined by its own observations alone. The random-effect version estimates a single width \(\sigma_u^2\), and the individual \(u_d\) are pulled towards the overall mean under that width. This pull, stronger for districts with fewer observations, is called shrinkage (partial pooling).

WHERE INDEPENDENCE IS USED — THE MARGINAL COVARIANCE
\[ \mathrm{Cov}\!\left(y_i,\,y_j\right) \;=\; \begin{cases} \sigma_u^2 & d(i)=d(j),\ i \neq j \\[4pt] 0 & d(i) \neq d(j) \end{cases} \qquad\Longrightarrow\qquad \mathrm{ICC} \;=\; \frac{\sigma_u^2}{\sigma_u^2+\sigma_\varepsilon^2} \]

The errors \(\varepsilon_i\) are mutually independent, but observations in the same district share \(u_d\) and are therefore correlated marginally. The size of that correlation is the intraclass correlation. The independence assumption has not vanished; it has been pushed down into the error layer.

The point is being able to name the line on which independence is used

The independence assumption of an ordinary regression is the claim that observations are unrelated to one another. A mixed model weakens it into independence conditional on district. The assumption has not disappeared; where it sits has changed. Pointing out, in each derivation, the line on which independence is used is a habit cultivated in Statistics II so that this kind of re-reading can be done for oneself. Correlation across districts — neighbouring districts resembling one another, say — is still set to 0 in this expression. Handling that requires a different structure (Section 5).

3Correspondence with Perspective 2 — what maximum likelihood maximises

The random effects are integrated out, and what remains is maximised.

Because \(u_d\) is a random variable rather than a parameter, it cannot be left in the likelihood. It is integrated out. What remains is a function of the fixed effects \(\boldsymbol{\beta}\) and the variances \(\sigma_u^2,\sigma_\varepsilon^2\), and that is maximised.

MARGINAL LIKELIHOOD — INTEGRATING OUT THE RANDOM EFFECTS
\[ L\!\left(\boldsymbol{\beta},\sigma_u^2,\sigma_\varepsilon^2\right) \;=\; \int p\!\left(\mathbf{y}\mid \mathbf{u}\right)\, p\!\left(\mathbf{u}\right)\, d\mathbf{u} \qquad\Longrightarrow\qquad \mathbf{y} \sim \mathcal{N}\!\left(X\boldsymbol{\beta},\; ZGZ^{\top}+R\right) \]

\(Z\) … the design matrix indicating which observation belongs to which district
\(G=\sigma_u^2 I\), \(R=\sigma_\varepsilon^2 I\) … the variances of the random effects and of the error
Only in the linear, normal case does this integral have a closed form, giving a multivariate normal distribution. For a non-normal response such as beta regression it does not, and numerical integration or MCMC is required.

From here the procedure is the standard one of Perspective 2. Differentiate the log-likelihood, set it to zero, solve numerically. The variability of the estimate follows from the information, and standard errors, confidence intervals and tests are assembled from it. The parameter is a fixed unknown value; what varies at random is the sample.

REML — correcting the bias in the variance estimate, at the cost of what may be compared

Estimating variance components by maximum likelihood produces a downward bias, because the fixed effects are estimated at the same time. Restricted maximum likelihood (REML) corrects for this by using a likelihood from which the influence of the fixed effects has been removed. The REML likelihood, however, depends on the design matrix of the fixed effects, so neither the REML likelihood nor the AIC can be compared between models with different fixed-effect structures (Pinheiro & Bates 2000). Using AIC to select a set of predictors requires refitting by maximum likelihood.

4What changes when the same structure is estimated in a Bayesian framework

The structure is the same. What differs is whether \(\sigma_u\) receives a point or a distribution.

In this research the same Maibara–Sekigahara data were later re-analysed with a hierarchical Bayesian model. That is what the Bayesian estimation page covers, and the (1|district) in its model formula is the same object as \(u_{d(i)}\) here. The structure of placing district as a random effect did not change. What changed is the estimation framework and the distribution placed on the response.

 Linear mixed model (this page)Hierarchical Bayes (Bayesian estimation)
Hierarchical structuredistrict as a random effectdistrict as a random effect (the same)
Distribution of the responsenormalbeta, with a logit link
Treatment of \(\sigma_u\)returns a single point estimatereturns a posterior distribution
Uncertainty in the coefficientsconfidence interval (hit rate of a procedure)credible interval (probability for the parameter)
Quantity that can be reportedp-value\(P(\beta>0)\)
Tool for selectionAIC, BICapproximate leave-one-out cross-validation

Two frameworks applied to the same data and the same hierarchical structure. On the change in the distribution of the response, see the beta regression page.

With few clusters, this difference starts to matter

There are five districts and twenty-five settlements. \(\sigma_u^2\) is therefore estimated from the spread of five values, which is very little information. In this situation maximum likelihood can return \(\hat{\sigma}_u = 0\), a boundary solution: since a variance cannot be negative, the maximum of the likelihood falls at the edge of the parameter space. An estimate of 0 does not mean that no district differences exist — it means only that the data could not determine the width — but the number printed is the same 0 either way.

In a Bayesian treatment a prior is placed on \(\sigma_u\), so the posterior does not collapse onto a single point at 0. It is returned as an interval, and one can read that it touches 0 while still having width. Indeed, for the spatial variance in the same research the ICAR page reports sdcar as 0.04 with a 95% credible interval of [0.00, 0.16] — a lower bound touching 0, yet displayed as an interval. Under maximum likelihood this would be a single number, 0, and whether it merely touches or genuinely has width could not be distinguished.

The figures from the two analyses must not be conflated

The AIC values and coefficients in the table in Section 7 are the output of maximum likelihood estimation under an assumed normal distribution. The \(\beta = +0.231\) and \(P(\beta>0)=0.993\) reported on the Bayesian estimation page are a different quantity, obtained from a posterior distribution under an assumed beta distribution. The response distribution, the link function and the principle of estimation all differ, so the magnitudes of the coefficients cannot be compared directly. They are outputs of different models, even though the research and the data are the same.

5Handling spatial autocorrelation within a frequentist frame — the spatial error model

Resemblance between neighbours, which districts cannot capture, parameterised as a correlation in the errors.

The marginal covariance in Section 2 set correlation across districts to 0. But settlements laid out across a landscape resemble one another when they are adjacent, whatever district they belong to. Spatial models address this structure. In the analysis of eighteen points in the upper Ibi River valley, a spatial error model was used.

SPATIAL ERROR MODEL — THE ERROR ITSELF IS SPATIALLY CORRELATED
\[ \mathbf{y} \;=\; X\boldsymbol{\beta} + \mathbf{u}, \qquad \mathbf{u} \;=\; \lambda W \mathbf{u} + \boldsymbol{\varepsilon} \qquad\Longrightarrow\qquad \mathbf{u} \;=\; \left(I-\lambda W\right)^{-1}\boldsymbol{\varepsilon} \]

\(W\) … the neighbourhood matrix, specifying which points count as adjacent
\(\lambda\) … the strength of spatial autocorrelation. At \(\lambda=0\) the model reduces to an ordinary regression
Unlike a random effect, the spatial structure is written as a correlation in the errors. Estimation is by maximum likelihood.

Contrast with ICAR — the same problem, written in a different layer

ICAR places a spatial random effect \(\phi_i\) at each point and gives \(\phi\) a prior encoding that neighbours take similar values. The spatial error model does not name \(\phi\) at all; it expresses the covariance structure of the errors through the single parameter \(\lambda\). The former writes the spatial structure as a prior within a Bayesian frame, the latter as a correlation structure within a frequentist one. The problem addressed is the same; what differs is the layer in which the structure is placed.

The definition of a neighbour is a construct here too

The neighbour list in this analysis is built not from straight-line distance but from walking-cost distance between settlements. For settlements in mountain valleys, two places close in a straight line may be far apart in practice if a ridge has to be crossed. A walking-cost function (Aiken 1977; Langmuir 1984) is used to correct for this. Which function to use, and where to cut off adjacency, are nonetheless decisions made by the analyst. The point made on the ICAR page — that what counts as a neighbour is not determined by the data — applies unchanged on the frequentist side.

6Model selection by AIC — a bridge to Perspective 5

Complexity is subtracted from fit. A quantity from information theory rides on top of a frequentist estimate.

AKAIKE INFORMATION CRITERION
\[ \mathrm{AIC} \;=\; -2\,\ell\!\left(\hat{\theta}\right) + 2k \]

\(\ell(\hat{\theta})\) … the maximised log-likelihood, the quantity maximised in Section 3
\(k\) … the number of estimated parameters
The first term is goodness of fit, the second a penalty on complexity (Akaike 1974).

AIC is a quantity of Perspective 5 (information theory). It estimates, from the data at hand, the wasted code length incurred when the true distribution is represented by the model — the Kullback–Leibler divergence. Fit by maximum likelihood, select by an information-theoretic criterion. The two perspectives are not in conflict; they are combined in a division of labour where estimation is frequentist and selection is information-theoretic. Where normalised compression distance, also under Perspective 5, measures nearness by how much can be compressed, AIC compares models by code length.

In Maibara–Sekigahara, selection was by exhaustive search over all subsets

The predictors are numerous: longitude, latitude, area, elevation, land use (built-up land, forest, paddy field), population, population density, population growth rate, ageing rate and the proportion under fifteen. In the GISA2016 analysis, MuMIn::dredge searched all subsets exhaustively, the lowest-AIC models were extracted, and models with \(\Delta\mathrm{AIC} < 4\) were averaged. The same procedure was carried out for gain and loss of the overall set and of the Kansai, Tokyo and Mino types separately. The hazards of this search are discussed in Section 8.

In the upper Ibi River valley, AIC decided whether a spatial term was needed

The AIC of the spatial error model is reported as \(-52.203\), against \(-50.466\) for a linear model with the same predictors and no spatial term. The difference of roughly 1.7 favours including the spatial term. AIC is a relative criterion, however, and what a difference of that size means is a separate judgement.

7How the method is used in this research

Two regions, two responses, two treatments of space.

(A) Maibara–Sekigahara — rates of vocabulary loss and gain against environmental and social factors

  1. Constructing the response For settlements in Maibara (Shiga Prefecture) and Sekigahara (Gifu Prefecture), the mean rates of loss and gain of forms from the cohort aged 60 and over (142 respondents) to the lower-secondary school cohort (544 respondents) are computed per settlement from a survey of 106 words across 40 items. Loss and gain are built separately for each of the three types: Kansai, Mino–Owari and Tokyo.
  2. Predictors Land use (the proportions of built-up land, forest and paddy field) is taken from the Land Use Subdivision Mesh of the National Land Numerical Information (Ministry of Land, Infrastructure, Transport and Tourism, FY2006); population, population density, population growth rate, ageing rate, proportion under fifteen and area from the small-area tabulations of the 2010 Population Census on e-Stat (Statistics Bureau, Ministry of Internal Affairs and Communications).
  3. Standardisation Predictors are standardised with scale so that the magnitudes of the coefficients can be compared.
  4. Model A linear mixed model is fitted with the lower-secondary school district as a random effect.
  5. Selection All subsets are searched exhaustively and the model with the lowest AIC is selected for each response.
Response (lowest-AIC model)AICBIClogLikdeviancedf.resid
Kansai type, loss70.0576.53−30.0360.0522
Tokyo type, loss63.3074.98−22.6545.3018
Mino type, loss75.2784.35−30.6461.2720
Kansai type, gain76.0381.22−34.0268.0323
Tokyo type, gain66.1275.19−28.0652.1220
Mino type, gain69.8280.19−26.9153.8219

Fit statistics for the lowest-AIC model for each response, from Table 1 of the GISA2016 abstract. Coefficients for eleven variables appear in the same table; only the findings the text states explicitly are recorded below.

The abstract states three findings. For every type, settlements with larger populations show lower rates of both loss and gain. The factor specific to loss is the population growth rate, and the factor specific to gain is area: where population is growing, words are more readily lost, and where a settlement covers a wider area, new words are more readily acquired. Further, for Tokyo-type and Mino–Owari-type items, both loss and gain are higher towards the larger longitudes, that is, on the Sekigahara side.

The abstract says 27 settlements, but the model appears to have taken 25

The df.resid in the table above can be read as the sample size minus the number of fixed effects, the intercept included. Kansai-type gain has one predictor plus an intercept and a df.resid of 23; Tokyo-type loss has six predictors plus an intercept and 18. Both are consistent with \(n=25\). The abstract records the study area as 27 settlements, but it is coherent to take the 25 settlements for which loss and gain rates were obtained as the ones entered into the model. The working record also states n=25, and the later Bayesian re-analysis uses 25 neighbourhood associations.

(B) The upper Ibi River valley — a composite variable from phylogenetic estimation against settlement conditions

  1. Study area Eighteen points in the upper Ibi River valley, north-west Gifu Prefecture: eight in Tokuyama, two in Fujihashi, three in Sakauchi and five in Kuze.
  2. Language data The 47 basic-vocabulary items reported by Yamada (1978, 1980).
  3. Phylogenetic estimation A phylogeny is estimated from the basic vocabulary. Branch support is assessed by 1000 bootstrap replicates, and the split separating Tokuyama from the rest is reported at 68.3%, that separating Tokuyama and Sakauchi from the rest at 95.2%.
  4. Building the composite variable Multidimensional scaling is applied to the distances produced by the phylogenetic estimation, and the first axis (MDS1), which maximises the variance of inter-settlement distance, is taken as the response.
  5. Predictors Number of households, population (1875–1920), paddy-field area, upland-field area, kokudaka (assessed rice yield) and elevation, all drawn from the Ibi County gazetteer and the Kadokawa dictionary of Japanese place names.
  6. Spatial error model A neighbour list is built from walking-cost distance and a spatial error model is fitted by maximum likelihood.
TermEstimateStd. errorzp
Intercept0.108437710.039873662.71950.0065374
Households0.003521940.001069483.29310.0009908
Population−0.000563700.00017643−3.19510.0013979
Paddy field−0.001341130.00067133−1.99770.0457472
Upland field−0.001156080.00037578−3.07650.0020944
Kokudaka0.000293810.000180601.62690.1037645
Elevation−0.000431520.00013314−3.24110.0011909

Partial regression coefficients of the spatial error model, quoted from the errorsarlm output printed in the presentation material. 18 observations, 9 estimated parameters. \(\lambda = -1.2505\); likelihood-ratio statistic 3.7368 (p = 0.053227); Wald statistic 24.445 (p = 7.648e−07); AIC −52.203 (−50.466 for the linear model).

The material reads this result as a relation in which, as MDS1 increases, households and kokudaka increase while paddy field, upland field and elevation decrease. For kokudaka the p-value of 0.104 somewhat exceeds 0.05, which the material describes as a level that may be regarded as a tendency towards significance.

\(\lambda\) is negative

The estimated \(\lambda\) is \(-1.2505\), a negative value. What a spatial model ordinarily anticipates is positive spatial autocorrelation — neighbours resembling one another — and a negative \(\lambda\) means the opposite, a structure in which neighbours differ. For settlements strung along a mountain valley, the material does not record what this sign reflects. The p-value of the likelihood-ratio test is also 0.053, so the need for the spatial term itself is not clearly supported. Eighteen observations against nine parameters is a further condition that counsels caution in interpreting this test.

8Limits — what cannot be claimed with this method

"The variance of the random effect is zero" cannot be put to an ordinary test

It is tempting to examine \(H_0:\sigma_u^2=0\) with a likelihood-ratio test, but this null hypothesis lies on the boundary of the parameter space. A variance cannot be negative, so the usual regularity condition — that the parameter is an interior point under the null — is not met. The likelihood-ratio statistic then does not follow a chi-squared distribution with one degree of freedom. Self & Liang (1987) give a general framework, and Stram & Lee (1994) show that for the longitudinal mixed model it is a 50:50 mixture of \(\chi^2_0\) and \(\chi^2_1\). Applying the standard chi-squared distribution errs on the conservative side, that is, in the direction of finding significance less readily.

AIC is a relative criterion and requires the same data and the same response

The absolute value of an AIC means nothing in itself. The 70.05 and the \(-52.203\) above say nothing on their own. What can be compared is only the difference between models fitted to the same data with the same response. If the handling of missing values differs and the number of rows changes, the comparison fails at that point. The same caution is taken up in Session 11 of Statistics I.

In addition, as stated in Section 3, the AIC of a model fitted by REML cannot be compared across models with different fixed effects. Using it to select a set of predictors requires fitting by maximum likelihood. Neither of the two sources consulted here records which was used.

An exhaustive search over all subsets grows more hazardous as the sample shrinks

With eleven predictors there are more than two thousand subsets. Computing an AIC for every one and taking the minimum carries a substantial risk of overfitting and false positives at a sample size of \(n=25\). The model selected is the combination that fitted these data best, which is not necessarily the best model in the population. Moreover, coefficients and p-values obtained after a selection step cannot be read as though no selection had taken place.

This is not an external criticism: the research's own working record states that "exhaustive search at n=25 carries a large risk of overfitting and false positives". The same record also lists the assumption of a normal distribution for proportion data as a misspecification of the distribution, and sets the move to beta regression as a task to be undertaken.

A normal distribution cannot accommodate a response that is a proportion

Loss and gain rates are proportions confined to \([0,1]\). A linear mixed model assuming a normal distribution has no structural means of preventing predicted values from leaving that interval. Even where the fit is good within the range of the sample, the specification is wrong. Replacing the response distribution with a beta distribution to address this is what beta regression does, and the later analyses in the same research take that form. What this page covers is the analysis that preceded it.

The definition of a neighbour, and the small number of clusters

The neighbourhood matrix of the spatial error model and the clusters of the mixed model (school districts) are alike frameworks chosen by the analyst: which walking-cost function, where to cut the distance, what unit to treat as a group. None of these is determined by the data. The same point is recorded on the ICAR page as the fact that, within this same body of research, changing the definition of adjacency produced different indications as to whether spatial autocorrelation was present.

The same holds for the number of clusters. There are five districts and either eighteen or twenty-five points. Estimating a variance component at this scale yields unstable estimates and tests with low power. In such a sample, "no district difference was detected" and "there is no district difference" are hard to tell apart.

9Sources and data

Two analyses are covered. (A) rests on a three-generation dialect vocabulary survey conducted across settlements from Maibara in Shiga Prefecture to Sekigahara in Gifu Prefecture. The original survey was carried out in FY2015 by Habiro Net (Kashiwabara, Maibara, Shiga) under the Maibara collaborative project proposal scheme, with 809 respondents in total. (B) covers eighteen points in the upper Ibi River valley in Gifu Prefecture, using the basic vocabulary reported by Yamada (1978, 1980).

References on the method

  • Akaike, H. (1974) A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. doi:10.1109/TAC.1974.1100705
  • Self, S. G. & Liang, K.-Y. (1987) Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. doi:10.1080/01621459.1987.10478472
  • Stram, D. O. & Lee, J. W. (1994) Variance components testing in the longitudinal mixed effects model. Biometrics, 50(4), 1171–1177.
  • Pinheiro, J. C. & Bates, D. M. (2000) Mixed-Effects Models in S and S-PLUS. Springer.

References on the research

  • Bromham, L., Hua, X., Fitzpatrick, T. G. & Greenhill, S. J. (2015) Rate of language evolution is affected by population size. Proceedings of the National Academy of Sciences, 112(7), 2097–2102. doi:10.1073/pnas.1419704112
  • Aoki, K., Lehmann, L. & Feldman, M. W. (2011) Rates of cultural change and patterns of cultural accumulation in stochastic models of social transmission. Theoretical Population Biology, 79(4), 192–202.
  • Trudgill, P. (2011) Sociolinguistic Typology: Social Determinants of Linguistic Complexity. Oxford University Press.
  • Tatsuya Yamada (1978) On the basic vocabulary of the upper Ibi River valley. Nagoya City University, Bulletin of the Faculty of Liberal Arts, Studies in Humanities and Social Sciences, 22, 33–66.
  • Tatsuya Yamada (1980) On the basic vocabulary of the Imajō and Ikeda dialects of Fukui Prefecture and the Kinomoto dialect of Shiga Prefecture, in relation to the Tokuyama dialect. Nagoya City University, Bulletin of the Faculty of Liberal Arts, Studies in Humanities and Social Sciences, 24, 23–34.
  • Aiken, R. (1977) Wilderness areas in Scotland. Unpublished Ph.D. thesis, University of Aberdeen. / Langmuir, E. (1984) Mountaincraft and leadership. The Scottish Sports Council/MLTB, Cordee. (Cited by (B) as the source of the walking-cost function. The body text of (B) writes "Aitken" while its reference list writes "Aiken".)

Which perspective this method belongs to

The linear mixed model treats the parameter as a fixed unknown value and maximises the marginal likelihood obtained by integrating out the random effects. What it returns is a point estimate and a confidence interval, and the probability attaches to the interval. It belongs to Perspective 2 (frequentist) of the seven perspectives in equations. AIC, used as the tool for selecting a model, is a quantity of Perspective 5 (information theory), where it sits alongside normalised compression distance. That the method runs on a combination of frequentist estimation and information-theoretic selection is what fixes its position. For the same data solved again under Perspective 3 (Bayesian), see Bayesian estimation and ICAR spatial random effects.

This page sets out a method Onohara uses in her own research, written for readers in the same field. Every figure is quoted from values printed in the presentation materials and from the working records of the analysis; no estimates or approximations have been supplied. Bibliographic details that could not be confirmed are not recorded.
For how the methods sit relative to one another see the seven perspectives in equations; for which course session falls under which perspective see Methods & Mathematics.