Multiple Regression and Variable Selection — Ayaka Onohara

1What the method is for

What is it about a settlement that explains "the number of accents that changed"?

For nine settlements on Shōdoshima there exists a figure obtained by setting the description in Nakai [1998] against a survey carried out in 2012 and counting how many of the accent types of two-mora nouns had changed. High mora H, low mora L, mid mora M, intra-mora rise R and intra-mora fall F were each replaced by the numerical codes 11, 00, 0.5 0.5, 01 and 10; the absolute difference was taken between corresponding cells on the 1998 side and the 2012 side; and these were summed for each locality (source pp.49–50, Tables 4-6 to 4-8). That sum is the number of changes.

The question is: what is the settlement-to-settlement variation in this number of changes bound up with? The explanatory variables the dissertation assembled are population (small-area tabulations of the 2010 and 2000 censuses), households (the same source), population growth rate (population 2010 ÷ population 2000), three kinds of centrality (information, betweenness and degree centrality), and elevation (the 5 m mesh digital elevation model of the Fundamental Geospatial Data) (source pp.47–48, p.52). The centralities were computed by building an adjacency matrix recording the connections between settlements as 0/1 from bus-route and scheduled passenger shipping-route data (Table 4-5), and passing it to degree, betweenness and infocent in R's sna package (source pp.50–52).

Multiple regression becomes necessary here because these explanatory variables are strongly correlated with one another. A settlement with a large population also has many households; a settlement where roads converge scores high on both betweenness and degree centrality. Lining up nine simple correlations does not reveal which variable is a restatement of which. Multiple regression attempts an answer in the form of "the effect of one variable with the others held fixed" — but when the variables resemble one another so closely that the holding-fixed does not bite, the coefficients become unstable. Everything from section 3 onwards concerns how that problem is handled.

2The model and least squares

The coefficients that minimise the residual sum of squares can be written uniquely as the solution of the normal equations.

Suppose there are $n$ observations and $p$ explanatory variables. For the $i$-th observation,

THE MULTIPLE REGRESSION MODEL

\[ y_i \;=\; \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + \varepsilon_i \qquad (i = 1,\dots,n) \]

$y_i$ … the response variable (here, the number of accent changes in settlement $i$)
$x_{ij}$ … the value of the $j$-th explanatory variable
$\varepsilon_i$ … the error term. Least squares itself assumes no distribution, but at the stage where the $t$ test and the $F$ test are carried out one assumes that the $\varepsilon_i \sim \mathcal{N}(0,\sigma^2)$ are mutually independent

With the design matrix $X$ (a first column of ones, then each explanatory variable) and the coefficient vector $\beta$, write $y = X\beta + \varepsilon$. Least squares takes the residual sum of squares

RESIDUAL SUM OF SQUARES AND THE NORMAL EQUATIONS

\[ S(\beta) \;=\; \sum_{i=1}^{n}\Bigl(y_i - \beta_0 - \textstyle\sum_{j=1}^{p}\beta_j x_{ij}\Bigr)^{2} \;=\; \lVert y - X\beta \rVert^{2} \] \[ \frac{\partial S}{\partial \beta} \;=\; -2X^{\top}\bigl(y - X\beta\bigr) \;=\; 0 \quad\Longrightarrow\quad X^{\top}X\,\hat{\beta} \;=\; X^{\top}y \] \[ \hat{\beta} \;=\; \bigl(X^{\top}X\bigr)^{-1}X^{\top}y \qquad \mathrm{Var}\bigl(\hat{\beta}\bigr) \;=\; \sigma^{2}\bigl(X^{\top}X\bigr)^{-1} \]

The second expression is the normal equations. If $X^{\top}X$ is non-singular the solution is uniquely determined. What the next section turns on is not whether the inverse exists but how large the diagonal entries of that inverse are.

The dissertation regresses on standardised data

The data of Table 4-9 were standardised with R's scale before being passed to lm (source Appendix. 5). The intercept therefore comes out as 0.000, and each coefficient states how many standard deviations the number of changes moves when the explanatory variable moves by one standard deviation. The price of being able to compare magnitudes across variables is that interpretation in the original units (persons, metres and so on) is lost.

3The problem of explanatory variables correlated with one another

Multicollinearity is not a problem of "cannot be estimated" but of "the variance of the estimate inflates".

In Figure 4-1 of the dissertation (the correlation matrix and scatterplot matrix of the standardised data of Table 4-9, p.54), the four variables population (2010), households (2010), population (2000) and households (2000) show mutual correlations of 0.998 to 0.999. The correlations among the centralities are high as well: p.53 states that "the correlation between betweenness centrality and degree centrality is likewise high", and, on the two-mora-noun-plus-particle side, p.62 states that "the correlation between information centrality and degree centrality is also high". What happens when variables this similar enter the model together is best seen in the form of the variance of the coefficients.

Take the $j$-th explanatory variable $x_j$ and regress it on all the remaining explanatory variables (including the intercept), writing the residual as $\tilde{x}_j$. By the Frisch–Waugh–Lovell theorem, the multiple regression coefficient $\hat{\beta}_j$ coincides with the simple regression coefficient obtained by regressing $y$ on $\tilde{x}_j$ alone.

THE COEFFICIENT, AND ITS VARIANCE

\[ \hat{\beta}_j \;=\; \frac{\tilde{x}_j^{\top}y}{\tilde{x}_j^{\top}\tilde{x}_j} \qquad\Longrightarrow\qquad \mathrm{Var}\bigl(\hat{\beta}_j\bigr) \;=\; \frac{\sigma^{2}}{\tilde{x}_j^{\top}\tilde{x}_j} \;=\; \frac{\sigma^{2}}{\mathrm{RSS}_j} \]

$\mathrm{RSS}_j$ … the residual sum of squares from regressing $x_j$ on the other explanatory variables. The size of what is left over once the other variables have explained what they can is precisely what fixes the precision of the coefficient.

Writing $R_j^2$ for the coefficient of determination of that regression of $x_j$ on the other explanatory variables, by definition $\mathrm{RSS}_j = \mathrm{TSS}_j\,(1-R_j^2)$, where $\mathrm{TSS}_j = \sum_i (x_{ij}-\bar{x}_j)^2$. Substituting,

VIF AS THE INFLATION FACTOR OF THE VARIANCE

\[ \mathrm{Var}\bigl(\hat{\beta}_j\bigr) \;=\; \frac{\sigma^{2}}{\displaystyle\sum_{i}\bigl(x_{ij}-\bar{x}_j\bigr)^{2}} \;\times\; \underbrace{\frac{1}{1-R_j^{2}}}_{\textstyle \mathrm{VIF}_j} \] \[ \mathrm{VIF}_j \;=\; \frac{1}{1-R_j^{2}} \]

the leading fraction … the variance that would obtain if $x_j$ were uncorrelated with the other variables ($R_j^2=0$)
$\mathrm{VIF}_j$ … by what factor the actual variance is inflated relative to that ideal case. On the scale of the standard error, the factor is $\sqrt{\mathrm{VIF}_j}$
$R_j^2$ … the coefficient of determination from a regression that takes $x_j$ as the response and the other explanatory variables as predictors. The response $y$ is not used at all

The dissertation writes this very formula into a footnote

Footnote 16 on p.55 reads: "vif takes the explanatory variable in question as the response, obtains the R² from regressing it on the other explanatory variables, and computes 1/(1 − R²)". The same footnote adds that "VIF = 10 is the state in which more than 90% of the variation of one explanatory variable is explained away by the others. Even at VIF = 5, more than 80% is explained away." If $R_j^2 = 0.9$ then $\mathrm{VIF}_j = 10$, and if $R_j^2 = 0.8$ then $5$; this is a restatement of the formula above.

The computation uses the vif function of R's car package (source p.53, Appendix. 5).

A large VIF does not make the estimate itself biased

Multicollinearity does not break the unbiasedness of the least squares estimator. What it breaks is precision alone. One therefore cannot say "the coefficient is wrong because the VIF is large", only "this coefficient could move a great deal if the sample were drawn again". The practice of dropping a variable once its VIF exceeds 10 buys a smaller variance at the price of an operation that pushes the effect of the dropped variable into the coefficients of the variables retained. The dissertation carries out this deletion in two rounds, and the coefficients of the final model have to be read on that understanding.

4All possible subset selection

Fit all $2^p$ combinations and sort them by information criterion.

With $p$ explanatory variables, sorting each into "used" or "not used" yields $2^p$ combinations, or $2^p - 1$ once the empty set is excluded. All possible subset selection is the procedure of running lm on every one of these, computing the coefficient of determination, the adjusted coefficient of determination, AIC, BIC and so on, and sorting the results. What the dissertation used is the function All.possible.subset.selection from Tamaki Yano's R script papca.R, reproduced in full in Appendix. 4 (A-11 to A-12).

THE SCOPE OF THE EXHAUSTIVE SEARCH

\[ \#\{\text{candidate models}\} \;=\; 2^{p} - 1 \]

In the case at hand, Table 4-9 has nine explanatory variables, so $2^{9}-1 = 511$ combinations. At the second stage, after betweenness centrality and others had been dropped, $p=6$ gives $2^{6}-1=63$; at the third stage $p=4$ gives $2^{4}-1=15$. By default the script displays only the twenty models with the smallest AIC.

The sorting criterion is AIC. For a log-likelihood $\log L$ and a number of parameters $k$,

AIC, AND THE OTHER QUANTITIES THE SCRIPT REPORTS

\[ \mathrm{AIC} \;=\; -2\log L\bigl(\hat{\theta}\bigr) + 2k \qquad \mathrm{BIC} \;=\; -2\log L\bigl(\hat{\theta}\bigr) + k\log n \] \[ \mathrm{cAIC} \;=\; \mathrm{AIC} + \frac{2(k+1)(k+2)}{n-k-2} \]

The third is what the script in Appendix. 4 computes as result7[i] <- result3[i] + 2*(nd+1)*(nd+2)/(nrow(df)-nd-2): a small-sample correction in $k$ (nd in the code, R's ans$rank). Besides these the script outputs $R^2$, the adjusted $R^2$, and a quantity called Rh.

How this differs from stepwise selection

Footnote 14 on p.53 reads as follows. "One generally uses stepAIC, but stepAIC does not compute every combination of the explanatory variables; it removes variables one at a time and halts the process at an appropriate point. All.possible.subset.selection, by contrast, displays every combination of the explanatory variables in ascending order of AIC."

Stepwise selection depends on the path of the search. Once a variable has been dropped, no combination containing it is ever evaluated again. Exhaustive search has no such path dependence and cannot miss the best model in the sense of AIC. But not missing it and the selected model being correct are two different things. The more candidates there are, the greater the probability that some combination fits well by accident and rises to the top.

AIC is not a test

AIC is a comparison from the standpoint of prediction (an estimate of the Kullback–Leibler divergence); it carries neither a significance level nor a null hypothesis. One therefore cannot say "this model is significant because its AIC is smallest". The dissertation applies the $t$ test and $F$ test of summary(lm(...)) to the model chosen by AIC; the interpretive problem that raises is taken up in section 8.

5The actual data, and the process of narrowing the variables

Nine settlements × 2 rows, nine explanatory variables. From this, two variables remain.

Table 4-9 of the dissertation (two-mora nouns only, p.52) is reproduced as it stands. There are 18 rows, but only nine settlements. Footnote 13 explains: "where there were two survey points within a single settlement, or where variation in pronunciation meant that two phonetic patterns were recognised at one and the same point, the weight of that point in the multiple regression would become large. To prevent this, the data row for each survey point that had a single location and a single phonetic pattern was duplicated." The values of the explanatory variables agree exactly between rows a and b of the same settlement.

Analysis name	Changes	Pop. 2010	Households 2010	Pop. 2000	Households 2000	Pop. growth rate	Information centrality	Betweenness centrality	Degree centrality	Elevation
Ieura a	2	1246	562	1560	676	0.799	0.592	0	4	5.5
Ieura b	2	1246	562	1560	676	0.799	0.592	0	4	5.5
Kamano a	0	294	142	378	149	0.778	0.816	0	4	6.38
Kamano b	0	294	142	378	149	0.778	0.816	0	4	6.38
Karato a	2	648	290	854	328	0.759	0.804	24	6	110.06
Karato b	5	648	290	854	328	0.759	0.804	24	6	110.06
Kusakabe a	4	1300	626	1525	663	0.852	0.997	13.333	8	2.13
Kusakabe b	4	1300	626	1525	663	0.852	0.997	13.333	8	2.13
Tonoshō a	2	5106	2188	5665	2261	0.901	1.098	31.667	8	8.58
Tonoshō b	2	5106	2188	5665	2261	0.901	1.098	31.667	8	8.58
Fukuda a	0	749	338	967	378	0.775	0.722	1	4	24.09
Fukuda b	0	749	338	967	378	0.775	0.722	1	4	24.09
Yoshida a	0	115	48	148	56	0.777	0.923	10.333	6	93.42
Yoshida b	0	115	48	148	56	0.777	0.923	10.333	6	93.42
Kō a	3	71	42	120	60	0.592	0.592	0	4	29.33
Kō b	3	71	42	120	60	0.592	0.592	0	4	29.33
Tōhama a	0	102	46	154	58	0.662	0.736	1.667	4	32.64
Tōhama b	0	102	46	154	58	0.662	0.736	1.667	4	32.64

Source Table 4-9, "Data for the multiple regression (two-mora nouns only)" (p.52). By footnote 12 of the source, [changes] corresponds to the "sum" in Table 4-8, and [population growth rate] is population 2010 ÷ population 2000.

The steps the dissertation followed are these (source pp.53–58, Appendix. 5).

Standardisation　Table 4-9 is standardised with scale (data1 → standardised).
Inspecting the correlations　The correlation matrix and scatterplot matrix (Figure 4-1) are drawn, confirming that the correlations among population (2010), households (2010), population (2000) and households (2000) are all above 0.9.
First exhaustive search　All.possible.subset.selection is run on the nine explanatory variables (Table 4-10). Of the four population-related variables, only households (2010) — the one that appears most readily in the low-AIC combinations — is retained; population (2010), population (2000) and households (2000) are removed to form data1.2.
Second exhaustive search　The search is run again on data1.2 (six explanatory variables) (Table 4-11). The minimum AIC is households (2010) + information centrality + betweenness centrality + degree centrality + elevation.
VIF and sensitivity analysis　vif and sensitivity are run on this model.
Third exhaustive search　Households (2010) and betweenness centrality, whose VIFs exceeded 10, are removed to form data1.3, and the search is run once more (Table 4-14). The minimum AIC is information centrality + degree centrality.

$R^2$	Adj.$R^2$	Rh	AIC	cAIC	BIC	Formula
0.881	0.832	0.787	25.712	36.912	31.944	households 2010 + information centrality + betweenness centrality + degree centrality + elevation
0.881	0.817	0.759	27.686	43.686	34.809	households 2010 + population growth rate + information centrality + betweenness centrality + degree centrality + elevation
0.767	0.736	0.708	31.823	34.900	35.384	information centrality + degree centrality
0.773	0.724	0.680	33.397	38.397	37.849	information centrality + degree centrality + elevation
0.772	0.723	0.680	33.424	38.424	37.876	population growth rate + information centrality + degree centrality

Only the top five rows of source Table 4-11, "Model selection for Data1.2 (two-mora nouns only)" (p.55); the source lists twenty. That the AIC gap between the minimum-AIC row and the third row is about 6.1, while on cAIC the ranking reverses — 36.912 against 34.900 — can be read directly off the source's own table.

6Results

What remains is two variables: information centrality (negative) and degree centrality (positive).

Second stage (data1.2, two-mora nouns only). The regression equation of the model chosen by AIC is given on p.54 of the dissertation in the following form.

THE REGRESSION EQUATION OF SOURCE p.54 (STANDARDISED DATA)

\[ \text{changes} = 0.000 - 1.097\,\text{households 2010} - 1.590\,\text{information centrality} + 1.523\,\text{betweenness centrality} + 1.383\,\text{degree centrality} - 0.895\,\text{elevation} \]

Coefficient of determination 0.881, adjusted coefficient of determination 0.832, $F$ = 17.8 ($p$ = $3.514\times10^{-5}$, 5 and 12 degrees of freedom), residual standard error 0.410.

Variable	Estimate	Std. error	$t$	$p$	VIF
Intercept	0.000	0.097	0.000	1.000	—
Households 2010	−1.097	0.343	−3.196	0.008	11.890
Information centrality	−1.590	0.218	−7.283	0.000	4.813
Betweenness centrality	1.523	0.470	3.238	0.007	22.331
Degree centrality	1.383	0.297	4.653	0.001	8.924
Elevation	−0.895	0.268	−3.346	0.006	7.226

Estimates, standard errors, $t$ and $p$ values from source Table 4-12, "Results of the multiple regression on data1.2" (p.56). The VIFs are from the text on p.55 and from the vif(RegModel.1) output in Appendix. 5 (A-17). The dissertation states that "there is an indicator holding that multicollinearity is present once VIF reaches 10, and on that indicator degree centrality, information centrality and elevation are the appropriate variables".

Third stage (data1.3, two-mora nouns only). Dropping households (2010) and betweenness centrality, whose VIFs exceeded 10, and repeating the exhaustive search, the minimum AIC is a two-variable model: information centrality + degree centrality.

Variable	Estimate	Std. error	$t$	$p$	VIF
Intercept	0.000	0.121	0.000	1.000	—
Degree centrality	1.895	0.271	6.987	0.000	4.738
Information centrality	−1.588	0.271	−5.855	0.000	4.738

Source Table 4-15, "Results of the multiple regression on data1.3" (p.57); the VIFs are from p.58 and Appendix. 5 (A-19). Coefficient of determination 0.767, adjusted coefficient of determination 0.736, $F$ = 24.7 ($p$ = $1.79\times10^{-5}$, 2 and 15 degrees of freedom), residual standard error 0.514 (source p.58).
Note that the regression equation in the running text on p.58 reads "changes = 0.000 + 1.901 degree centrality − 1.595 information centrality", whose final digits do not agree with Table 4-15 or with the summary output in Appendix. 5 (1.90 / −1.59, $t$ = 6.99 / −5.85). The values of Table 4-15 are the ones given here.

The two-mora-noun-plus-particle case. An analysis was also run on the same nine settlements with the same explanatory variables, replacing only the response by "the number of changes including the forms followed by a particle" (source Table 4-18). At the second stage (data2.2) the minimum AIC was households (2010) + population growth rate + information centrality + betweenness centrality + degree centrality + elevation, with coefficient of determination 0.896, adjusted coefficient of determination 0.839, $F$ = 15.77 ($p$ = $8.00\times10^{-5}$) and residual standard error 0.401 (source p.63).

Variable	Estimate	Std. error	$t$	$p$	VIF
Intercept	0.000	0.095	0.000	1.000	—
Households 2010	−1.394	0.359	−3.880	0.003	13.634
Population growth rate	0.214	0.170	1.253	0.236	3.064
Information centrality	−1.586	0.228	−6.952	0.000	5.495
Betweenness centrality	1.291	0.473	2.732	0.020	23.593
Degree centrality	1.601	0.295	5.420	0.000	9.219
Elevation	−0.799	0.264	−3.025	0.012	7.370

Source Table 4-21, "Results of the multiple regression on data2.2" (p.64); the VIFs are from the text on p.63. Population growth rate is the only variable whose $p$ does not fall below 0.05.

Here too, dropping households (2010) and betweenness centrality for VIFs above 10 and repeating the exhaustive search (data2.3), the minimum AIC is again information centrality + degree centrality: coefficient of determination 0.727, adjusted coefficient of determination 0.69, $F$ = 20 ($p$ = $5.94\times10^{-5}$), residual standard error 0.556, with a VIF of 4.738 for both variables (source pp.66–67).

Variable	Estimate	Std. error	$t$	$p$
Intercept	0.000	0.131	0.000	1.000
Degree centrality	1.846	0.294	6.283	0.000
Information centrality	−1.552	0.294	−5.283	0.000

Source Table 4-24, "Results of the multiple regression on data2.3" (p.67).

The direction of the signs is the core of the dissertation's claim

The coefficient on degree centrality is positive, that on information centrality negative. Page 71 reads this as "settlements where roads converge underwent large accent change between 1998 and 2012" and "settlements that are easy to reach changed little in accent between 1998 and 2012", and connects it to Kindaichi's [2003] argument that "urban language has parts that resist change, while regional language has parts that change readily". At the same time the dissertation notes that "there is also a view that language changes more readily in cities, so one cannot simply equate the urban with ease of access".

The dissertation goes further, using the magnitude of the residuals (Figure 4-4, Table 4-17) to separate out which accent changes are influenced by extra-linguistic factors. Changes appearing only in settlements with small residuals (LLF→LHL, MLR→LLH, MLF→LHL) are set against changes appearing only in settlements with large residuals (MF→HL, MFL→HLL, LLR→LHL), and the former are assigned to contact-induced change, the latter to internal change (source pp.72–73).

7Sensitivity analysis

Perturb the explanatory variables and see how far the coefficients move.

Alongside VIF the dissertation uses a second diagnostic for multicollinearity: sensitivity analysis. The tool is the sensitivity function of R's accuracy package; p.53 states that "the function sensitivity in the package accuracy, which performs a Sensitivity Analysis that judges multicollinearity by applying perturbations to the explanatory variables, was used". The actual call in Appendix. 5 is data1.2.sens <- sensitivity(data1.2, lm, 変化数 ~ ., ptb.R = 2000), and the heading of the output is "Sensitivity of coefficients over 2000 perturbations:". The text on p.55 likewise says that "the check was carried out with 2000 perturbations".

The definition of Unstable

Footnote 17 on p.57 states that "a variable is Unstable when, after perturbing the explanatory variables and repeating this an arbitrary number of times, a test of how far each of the resulting estimates falls within the 95% confidence interval is rejected"; in that case an asterisk appears in the [Unstable] column of the output.

Variable	Perturbed estimate	Estimate	Std. error	2.50%	97.50%	Unstable
Intercept	0.000	0.000	0.100	−0.004	0.004	—
Households 2010	−1.100	−1.100	0.380	−1.100	−1.060	—
Population growth rate	−0.023	−0.023	0.180	−0.030	−0.015	—
Information centrality	−1.600	−1.600	0.240	−1.588	−1.570	—
Betweenness centrality	1.500	1.500	0.500	1.485	1.533	—
Degree centrality	1.400	1.400	0.320	1.378	1.403	—
Elevation	−0.890	−0.890	0.280	−0.903	−0.877	—

Source Table 4-13, "Results of the Sensitivity Analysis on data1.2 (1)" (p.56). The "Unstable" column is blank throughout, and p.57 states that "the result obtained was that no variable was unstable (Unstable)". The same output is also printed in Appendix. 5 (A-18) with digits = 2, and there are places where the final digit of the 2.5% / 97.5% values disagrees with Table 4-13 (for example, the 2.5% figure for households 2010 is −1.100 in Table 4-13 and −1.099 in Appendix. 5). Comparable output for data2.2, data1.3 and data2.3 appears in source Tables 4-22, 4-16 and 4-25, and in none of them is a verdict of instability returned.

The model handed to the sensitivity analysis is not the model chosen by AIC

The call in Appendix. 5 is sensitivity(data1.2, lm, 変化数 ~ ., ptb.R = 2000). Since 変化数 ~ . means "take every column remaining in the data frame as an explanatory variable", what is fitted here is the model containing all six variables of data1.2, not the five-variable model chosen by AIC (which excludes population growth rate). That is why Table 4-13 has a row for population growth rate. The same holds for data1.3: Table 4-16 has rows for population growth rate and elevation, whereas the model chosen by AIC is the two-variable model information centrality + degree centrality.

In other words, the object the sensitivity analysis judged "stable" and the object whose coefficients and $p$ values are reported are not the same model. The dissertation does not mention this difference.

The two diagnostics disagree

VIF judged households 2010 = 11.890 and betweenness centrality = 22.331 to be "above 10", yet the sensitivity analysis found no variable Unstable. The text on p.57 reads "households H22 and betweenness centrality, judged unstable by vif (>10), were removed from data1.2", so the VIF verdict was given priority in dropping variables. Which diagnostic one adopts changes the final model, and the dissertation does not state the reason for the choice.

What could not be determined

Within the dissertation and Appendices 4 to 6, neither the distribution nor the magnitude of the "perturbation" that sensitivity applies to the explanatory variables is recorded. Only the number of repetitions, ptb.R = 2000, is specified; the rest relies on the default perturbation settings. The accuracy package is moreover already described, in footnote 15 on p.53, as "no longer being developed and unobtainable from CRAN".

8Limits

There is a problem in the procedure itself: selecting by exhaustive search and then testing.

Testing after selecting makes the $p$ values too generous

All possible subset selection fits $2^p-1$ models on one and the same sample and picks the one that fits best. When the $t$ test and $F$ test are then carried out on that selected model using the same sample, the reported $p$ values lose their nominal meaning. Both the $t$ and the $F$ distribution are derived on the premise that the model was fixed before the data were seen, and the operation of looking at the data and then picking the best breaks that premise. The $p$ values actually reported are smaller than they should be, and the confidence intervals narrower (the problem of post-selection inference).

At the second stage, 63 models were compared from $p=6$; at the first stage, 511 from $p=9$. Attaching a $p$ value of $3.514\times10^{-5}$ to the $F$ value of 17.8 of the model that was best among 511 means carrying over unchanged the calculation appropriate to having looked at a single model. The dissertation does not address this point.

There are many variables relative to the sample size, and the rows are duplicated

Table 4-9 has 18 rows, but there are nine settlements, and the values of the explanatory variables agree exactly between rows a and b of the same settlement. Footnote 13 explains this duplication as a device for equalising the weight against the case where "there were two survey points within a single settlement". The response (the number of changes) can differ between a and b for a given settlement, but on the explanatory side there are only nine distinct sets of values.

Despite this, the regression is run on 18 rows of data, and the residual degrees of freedom are computed as 12 (second stage) or 15 (third stage). Since there are only nine independent settlements, these degrees of freedom are very likely too large. The bias runs towards standard errors that are too small, $t$ values too large and $p$ values too small. The dissertation does not address this point.

In addition, at the first stage there are nine explanatory variables and nine independent settlements. Searching 511 combinations under that ratio leaves ample room to pick up an accidental fit.

The observations are spatially arranged, and the independence of the errors is doubtful

The nine settlements are laid out on Shōdoshima. The explanatory variables are themselves centralities built from the road and shipping-route connections between settlements (Table 4-5). Settlements close to one another are likely to resemble each other both in the movement of people and in dialect, and in that case the error terms of the regression are not mutually independent. Positive spatial autocorrelation in the errors leads to underestimated standard errors.

The dissertation plots the absolute values of the residuals as circles sized over the settlement locations (Figure 4-5), so a spatial sketch is provided. It does not, however, test the residuals for spatial autocorrelation, nor build the spatial structure into the model. The dissertation does not state this as a limitation. For the position that puts spatial correlation of the errors into the model explicitly, see ICAR spatial random effects.

The data are not treated as nested

The 18 rows have a nested structure: nine settlements × 2 rows. The two rows of a settlement share the same explanatory variables and derive from the same body of speakers. Treating each row of such data as an independent observation ignores the within-settlement correlation. The straightforward remedy is to treat settlement as a random effect and estimate the within- and between-settlement variance separately, a framework covered under linear mixed models and AIC. The dissertation proceeds with a fixed-effects-only multiple regression and does not mention this option.

The response is a non-negative count

The number of changes is a value of 0 or more, falling in the range 0 to 5 in Table 4-9 and 0 to 10 in Table 4-18. Because mid mora M is replaced by 0.5 during the computation, a cell-level difference can take the value 0.5, but the per-locality sums listed in Tables 4-9 and 4-18 are all integers. Either way, there is a floor at 0, and the variance may change with the level. A linear model with normal errors is treating this as an approximation. That the fitted values can go negative is a further caution for extrapolating from the model. The dissertation does not address this point.

The two analyses are not independent replications

That information centrality + degree centrality was selected in both the "two-mora nouns only" and the "two-mora nouns plus particle" analyses is used on p.71 as support for the conclusion. But the two use the same nine settlements and the same explanatory variables, and both responses are built from the same survey of the same words (water, sound, mountain, needle, autumn). This is not replication on an independent sample, so the agreement cannot be taken to double the credibility. The dissertation does not address this point.

This is not a causal claim

That a regression coefficient expresses the association "a settlement one standard deviation higher in degree centrality has a number of changes 1.9 standard deviations higher" is one thing; "building more roads will change the accent" is another. Variables that could confound the relation — the industrial structure of the settlement, actual patterns of commuting to school and work, the age composition of the speakers, and differences in the accent system already present in 1998 — are not in the model. The dissertation itself enters a reservation where it connects the results to Kindaichi's [2003] argument: "one cannot simply equate the urban with ease of access".

On the VIF threshold itself

"VIF > 10 means multicollinearity" is a rule of thumb, not a test. As footnote 16 of the dissertation writes, it says no more than $R_j^2 > 0.9$. Set the threshold at 5 and everything at the second stage except information centrality (4.813) would be dropped. Since where one cuts changes the final model, the threshold ought to be reported as part of the result.

9Sources

Every figure on this page is quoted from Chapter 4 and the end matter Appendix of Ayaka Onohara, A study of language change and its relation to extra-linguistic factors from a mathematical approach (Graduate School of Culture and Information Science, Doshisha University, doctoral dissertation, academic year 2013). Page numbers in the body text refer to the continuous pagination of the dissertation itself; Appendix page numbers refer to the end matter, prefixed "A-".

Onohara, Ayaka (2013) A study of language change and its relation to extra-linguistic factors from a mathematical approach. Doctoral dissertation, Graduate School of Culture and Information Science, Doshisha University.
　Chapter 4, section 4, "Method — centrality and multiple regression" pp.49–52 ／ section 5, "Results" pp.53–70 ／ section 6, "Discussion" pp.71–73.
　Table 4-9 (analysis data, p.52), Tables 4-10 and 4-11 (model selection, p.55), Tables 4-12 and 4-13 (coefficients and sensitivity analysis, p.56), Tables 4-14 and 4-15 (the same, p.57), Table 4-18 (data with particle, p.62), Table 4-21 (coefficients, p.64), Table 4-22 (sensitivity analysis, p.65), Tables 4-24 and 4-25 (the same, p.67).
　Appendix. 4, "Case 2: the All.possible.subset.selection script" A-11 to A-12 ／ Appendix. 5, "Case 2: R script for the multiple regression (two-mora nouns)" A-13 onwards ／ Appendix. 6, "Case 2: R script for the multiple regression (two-mora nouns plus particle)".

References for the methods and tools cited in the dissertation

Stephenson, K. & Zelen, M. (1989) Rethinking Centrality: Methods and Applications. Social Networks, 11, 1–37.　Cited on p.52 as the source for information centrality.
Suzuki, Tsutomu (2009) Centrality. Data Science with R: Network Analysis (vol. 8), Kyoritsu Shuppan.　Cited on pp.51–52 as the source for the accounts of degree, betweenness and information centrality.
Butts, C. T. (2010) Package ‘sna’.　The R package used to compute the centralities.
The vif function of the R package car (source p.53).
The sensitivity function of the R package accuracy (source p.53, footnote 15). At the time of writing it is noted as no longer developed and unobtainable from CRAN.
The function All.possible.subset.selection in Tamaki Yano's R script papca.R (source p.53, footnote 14). The comment at the head of the script reads “originaly by S.Aoki / modified by T.Yano”.

Provenance of the data

Accent: Nakai, Yukihiko (1998) "The accent of the eastern islands of Kagawa Prefecture", 61, 1–26, together with the survey of March 2012 (JSPS Grant-in-Aid, Kiban (A) 2324202, "A detailed account of change in dialect distribution", principal investigator Takuichiro Onishi; led by Professor Shinsuke Kishie, Japanese Linguistics Laboratory, Tokushima University). Source p.47, footnote 10.
Population and households: Statistics Bureau, Ministry of Internal Affairs and Communications, "2010 Population Census, small-area tabulation, Table 2: population by sex and number of households — by town, block and district", and "2000 Population Census (small area): total population by sex and total number of households". Source p.48.
Connections between settlements: Ministry of Land, Infrastructure, Transport and Tourism, "National Land Numerical Information: bus routes" (SHAPE format, world geodetic system, Kagawa Prefecture, 2011) and "National Land Numerical Information: scheduled passenger shipping routes" (2012). Source p.48.
Elevation: Geospatial Information Authority of Japan, MLIT, "Fundamental Geospatial Data" digital elevation model (JPGIS format, 5 m mesh). Source p.48.

Which perspective this method belongs to

Multiple regression treats the true coefficient $\beta$ as a fixed unknown constant. It is the estimator $\hat{\beta}$ that is the random quantity, and its fluctuation is evaluated as sampling variation: how far would the coefficient move if the sample were drawn again from the same settlements? The standard error $\sqrt{\sigma^2 (X^{\top}X)^{-1}_{jj}}$, the $t$ value and the $p$ value are all statements about hit rates under that repetition. What VIF measures is likewise nothing other than the factor by which the coefficient scatters within that repetition. This coincides exactly with Perspective 2 (frequentist) of the seven perspectives in equations.

The position that places a distribution on the coefficients instead and speaks in terms of priors and posteriors (Bayesian estimation), the position that handles nested structure as random effects (linear mixed models and AIC), and the position that builds spatial dependence itself into the model (ICAR spatial random effects) apply a different yardstick to the very same data.

This page sets out a method Onohara used in her own research, written for readers in the same field. Every figure is quoted from the body text and end matter Appendix of the doctoral dissertation; no estimates or approximations have been supplied. Where values within the source do not agree, this has been stated.
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.