Model Complications

Sociology 312/412/512, University of Oregon

Aaron Gullickson

The Linear Model Revisited

Review of linear model

${\hat{y}}_{i} = β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + \dots + β_{p} x_{i p}$

${\hat{y}}_{i}$ is the predicted value of the outcome/dependent variable which must be quantitative.
The $x$ values are different independent/explanatory variables which may be quantitative or categorical (by using dummy coding).
The $β$ values are the slopes/effects/coefficients that, holding all other independent variables constant, tell us the:
- The predicted change in $y$ for a one unit increase in $x$ if $x$ is quantitative
- The mean difference in $y$ between the indicated category and the reference category if $x$ is categorical.

Interpret these values

coef(lm(nominations~nsports+I(parent_income-45)+gender, data=popularity))

          (Intercept)               nsports I(parent_income - 45) 
           4.29062324            0.50750823            0.01071398 
           genderMale 
          -0.62846269

The residual term

$ϵ_{i} = y_{i} - {\hat{y}}_{i}$

The residual/error term $ϵ_{i}$ gives us the difference between the actual value of the outcome variable for a given observation and the value predicted by the model.

Lets use algebra to rewrite this equation with $y_{i}$ on the left-hand side:

$y_{i} = {\hat{y}}_{i} + ϵ_{i}$

If we plug in our linear model formula for ${\hat{y}}_{i}$ , we can get the full model formula:

$y_{i} = \underset{{\hat{y}}_{i}}{\underset{⏟}{β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + \dots + β_{p} x_{i p}}} + ϵ_{i}$

Linear model as a partition of the variance in $y$

$y_{i} = \underset{structural}{\underset{⏟}{β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + \dots + β_{p} x_{i p}}} + \underset{stochastic}{\underset{⏟}{ϵ_{i}}}$

The structural part is the predicted value from our model which is typically a linear function of the independent variables.
The stochastic component is the leftover residual or error component, that is not accounted for by the model.

Depending on disciplinary norms, there are different conceptual ways to view this basic relationship:

Description: observed = summary + residual
Prediction: observed = predicted + error
Causation: observed = true process + disturbance

Calculating marginal effects

The marginal effect of $x$ on $y$ is simply the predicted change in $y$ for a one unit increase in $x$ from its current value, holding all else constant.

In a basic linear model, the marginal effect is just given by the slope itself.

$y_{i} = β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + ϵ_{i}$

$β_{1}$ is the marginal effect of $x_{1}$
$β_{2}$ is the marginal effect of $x_{2}$

Technically, the marginal effect is the derivative of $y$ with respect to a given $x$ . This gives us the tangent line for the curve at any value of $x$ .

$\frac{\partial y}{\partial x_{1}} = β_{1}$ $\frac{\partial y}{\partial x_{2}} = β_{2}$

😱 I am not expecting you to know the calculus behind this, but it may help some people.

Marginal effects with interaction terms

Marginal effects can get more complicated when we move to more complex model. Consider this model with an interaction term:

$y_{i} = β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + β_{3} (x_{i 1}) (x_{i 2}) + ϵ_{i}$

The marginal effects are now given by:

$\frac{\partial y}{\partial x_{1}} = β_{1} + β_{3} x_{i 2}$ $\frac{\partial y}{\partial x_{2}} = β_{2} + β_{3} x_{i 1}$

This is a very math-y way of saying: the main effects depend on the effect of the other variable.

model <- lm(nominations~nsports*gender, data=popularity)
as.data.frame(coef(model))

                   coef(model)
(Intercept)         4.21374011
nsports             0.59169144
genderMale         -0.54255556
nsports:genderMale -0.08760095

The marginal effect for number of sports played is:

$0.49 - 0.09 (m a l e_{i})$

The marginal effect for gender is:

$- 0.54 - 0.09 (n s p o r t s_{i})$

How are linear model parameters estimated?

⚠️ Heavy math ahead!

We have simple formulas for the slope and intercept for a bivariat model.
With multiple independent variables, a simple formula will not suffice. To estimate model parameters with multiple independent variables we need to use some matrix algebra.

The matrix algebra approach to linear models

We can use matrix algebra to represent our linear regression model equation using one-dimensional vectors and two-dimensional matrices.

$y = X β + ϵ$

$y$ is a vector of known values of the independent variable of length $n$ .
$X$ is a matrix of known values of the independent variables of dimensions $n$ by $p + 1$ .
$β$ is a vector of to-be-estimated values of intercepts and slopes of length $p + 1$ .
$ϵ$ is a vector of residuals of length $n$ that will be equal to $y - X β$ .

⚠️ Minimize sum of squared residuals

Our goal is to choose a $β$ vector that minimizes the sum of squared residuals, $S S R$ , which is just given by the $ϵ$ vector squared and summed up. We can rewrite the matrix algebra formula to isolate $e$ on one side:

$\begin{aligned} y & = & X β + ϵ \\ ϵ & = & y - X β \end{aligned}$

In matrix form, $S S R$ can be represented as a function of $β$ , like so:

$\begin{aligned} S S R (β) & = & (y - X β)^{'} (y - X β) \\ = & y^{'} y - 2 y^{'} X β + β^{'} X^{'} X β \end{aligned}$

⚠️ Minimize sum of squared residuals

$\begin{aligned} S S R (β) & = & (y - X β)^{'} (y - X β) \\ = & y^{'} y - 2 y^{'} X β + β^{'} X^{'} X β \end{aligned}$

We want to choose a $β$ to plug into this function that provides the smallest possible value (the minimum).

It turns out that we can get this value by using calculus to get the derivative with respect to $β$ and solving for zero:

$0 = - 2 X^{'} y + 2 X^{'} X β$

Applying some matrix algebra will give us the estimator of $β$ :

$β = (X^{'} X)^{- 1} X^{'} y$

⚠️ Estimating standard errors

If we treat $σ^{2}$ as the variance of the error term $ϵ$ , then we can also use matrix algebra to calculate the covariance matrix:

$σ^{2} (X^{'} X)^{- 1}$

The values of this matrix give us information about the correlation between different independent variables. Most importantly, the square root of the diagonal values of this matrix are the standard errors for the estimated values of $β$ .

In practice, we don’t have $σ^{2}$ , but we can estimate from the fitted values of $y$ by:

$s^{2} = \frac{\sum (y_{i} - {\hat{y}}_{i})^{2}}{n - p - 1}$ We can then use these estimated standard errors to calculate t-statistics and p-values, confidence intervals, and so on.

Modeling Non-Linearity

Title	IMDB Rating	Box Office (millions)	Smoothed Box Office (median)
Nobel Son	6.2	$0.5	-
Resident Evil: Extinction	6.2	$50.6	-
Rush Hour 3	6.2	$140.1	$50.6
Spider-Man 3	6.2	$336.5	-
The Last Mimzy	6.2	$21.5	-

The natural log

Although ggplot uses log with a base 10, we usually use the natural log transformation in practice. Both transformations have the same effect on the relationship, but the natural log provides results that are easier to interpret.

In R, it is easy to take the natural log of a number by just using the log command. Any positive number can be logged.

log(5)

[1] 1.609438

The natural log of 5 is 1.609, but what does this mean?

The natural log of any number is the power that you would have to raise the constant $e$ to in order to get the original number back. In R, we can calculate $e^{x}$ with exp(x):

exp(log(5))

[1] 5

What does it mean for the model?

To fit the model, we can just use the log transformation directly in the formula of our lm command:

model <- lm(log(box_office)~rating_imdb, data=movies)
coef(model)

(Intercept) rating_imdb 
 0.01557736  0.37140161

$\log ({\hat{r e t u r n s}}_{i}) = 0.016 + 0.371 (r a t i n g_{i})$

How do we interpret the slope?

A one point increase in IMDB rating is associated with a 0.371 increase in … the log of box office returns? 😕

Converting to the original scale

To get back to the original scale of box office returns for the dependent variable, we need to exponentiate both side side of the regression equation by $e$ :

$e^{\log ({\hat{r e t u r n s}}_{i})} = e^{0.016 + 0.371 (r a t i n g_{i})}$

On the left hand side, I will get back predicted box office returns by the definition of logs. On the right hand side, I can apply some of the mathematical properties of logarithms and powers.

${\hat{r e t u r n s}}_{i} = (e^{0.016}) (e^{0.371})^{r a t i n g_{i}}$

exp(0.016)

[1] 1.016129

exp(0.371)

[1] 1.449183

I now have: ${\hat{r e t u r n s}}_{i} = (1.016) (1.449)^{r a t i n g_{i}}$

A multiplicative relationship

${\hat{r e t u r n s}}_{i} = (1.016) (1.449)^{r a t i n g_{i}}$

Lets calculate predicted box office returns for IMDB Ratings of 0, 1, and 2.

${\hat{r e t u r n s}}_{i} = (1.016) (1.449)^{0} = 1.016$

${\hat{r e t u r n s}}_{i} = (1.016) (1.449)^{1} = 1.449 (1.449)$

${\hat{r e t u r n s}}_{i} = (1.016) (1.449)^{2} = 1.449 (1.449) (1.449)$

For each one unit increase in the independent variable, you multiply the previous predicted value by 1.449 to get the new predicted value. Therefore the predicted value increases by 44.9%.
The model predicts that movies with a zero IMDB rating make 1.02 million dollars, on average. Every one point increase in the IMDB rating is associated with a 44.9% increase in box office returns, on average.

General form and interpretation

$\log ({\hat{y}}_{i}) = b_{0} + b_{1} (x_{i})$

You must apply the exp command to your intercept and slopes in order to interpret them.
The model predicts that the mean of $y$ when $x$ is zero will be $e^{b_{0}}$ .
The model predicts that each one unit increase in $x$ is associated with a multiplicative change of $e^{b_{1}}$ in $y$ . It is often easiest to express this in percentage terms.
An absolute change in $x$ is associated with a relative change in $y$ .

Interpret these numbers

coef(lm(log(box_office)~relevel(maturity_rating, "PG"), data=movies))

                        (Intercept)     relevel(maturity_rating, "PG")G 
                          3.2822510                           0.7482935 
relevel(maturity_rating, "PG")PG-13     relevel(maturity_rating, "PG")R 
                         -0.2725515                          -1.7874802

$\log ({\hat{r e t u r n s}}_{i}) = 3.22 + 0.50 (G_{i}) - 0.24 (P G 13_{i}) - 1.93 (R_{i})$

$e^{3.28} = 26.58$ : PG-rated movies make 26.58 million dollars on average.
$e^{0.75} = 2.12$ : G-rated movies make 112% more than PG-rated movies, on average.
$e^{- 0.27} = 0.76$ : PG-13 rated movies make 76% as much as (or 24% less than) PG-rated movies, on average.
$e^{- 1.79} = 0.17$ : R-rated movies make 17% as much as (or 83% less than) PG-rated movies, on average.

Approximating small effects

When $x$ is small (say $x < 0.2$ ), then $e^{x} \approx 1 + x$ .

We can use this fact to roughly approximate coefficients/slopes as percent change when they are small.

model <- lm(log(box_office)~runtime, data=movies)
coef(model)

(Intercept)     runtime 
 -1.7528129   0.0383472

exp(coef(model))

(Intercept)     runtime 
  0.1732858   1.0390919

If we do the full exponentiating, we can see that a one minute increase in runtime is associated with a 3.9% increase in box office returns.
The actual percentage increase is very close to what we got for the slope of the non-exponentiated slope (0.038). So, you can often get a ballpark estimate without having to exponentiate.

Logging the independent variable

Lets return to the relationship between GDP per capita and life expectancy that fit well as a linear relationship when we logged GDP per capita. Lets run the model:

model <- lm(lifeExp~log(gdpPercap), 
            data=subset(gapminder, 
                        year==2007))
round(coef(model), 5)

   (Intercept) log(gdpPercap) 
       4.94961        7.20280

How do we interpret these results? This case requires something different than the case where we logged the dependent variable.

Our basic model for life expectancy by GDP per capita is:

${\hat{y}}_{i} = 4.9 + 7.2 \log x_{i}$ What is the predicted value of life expectancy at $1 GDP?

${\hat{y}}_{i} = 4.9 + 7.2 \log 1 = 4.9 + 7.2 * 0 = 4.9$

What happens when we increase GDP per capita by 1% (from 1 to 1.01)?

${\hat{y}}_{i} = 4.9 + 7.2 \log 1.01 = 4.9 + 7.2 * 0.01 = 4.9 + 0.072$

Predicted life expectancy increases by 0.072 years.

The elasticity model

model <- lm(log(wages)~log(age), data=earnings)
coef(model)

(Intercept)    log(age) 
  1.1558386   0.5055849

This is actually the easiest model to interpret. We can interpret the slope directly as the percent change in $y$ for a one percent increase in $x$ .
Calculating the percent change in one variable by percent change in another is what economists call an elasticity so this model is often called an elasticity model.
The model predicts that a one percent increase in age is associated with a 0.51% increase in wages, on average.

Which variable logged	Non-linear shape	Change in x	Change in y	Interpret $β_{1}$
Independent variable	diminishing returns	relative	absolute	$β_{1} / 100$
Dependent variable	exponential	absolute	relative	$e^{β_{1}}$
Both variables	both types	relative	relative	$β_{1}$

Polynomial models

Polynomial expression

A polynomial expression is one that adds together terms involving multiple powers of a single variable.

if we include a squared value of $x$ we get the classic formula for a parabola:

$y = a + b x + c x^{2}$

A polynomial model

We can fit such a parabola in a linear model by including a new variable that is simply the square of the original variable:

${\hat{y}}_{i} = β_{0} + β_{1} x_{i} + β_{2} x_{i}^{2}$

model <- lm(wages~I(age-40)+I((age-40)^2), 
            data=earnings)
coef(model)

    (Intercept)     I(age - 40) I((age - 40)^2) 
    26.92134254      0.31711167     -0.01783204

Our model is:

${\hat{y}}_{i} = 26.9 + 0.3171 (x_{i} - 40) - 0.0178 (x_{i} - 40)^{2}$

How do we interpret the results?

Calculate the marginal effect

With polynomial terms, the marginal effect of $x$ is more complicated because our $x$ shows up in more than one place in the equation:

$y_{i} = β_{0} + β_{1} x_{i} + β_{2} x_{i}^{2} + ϵ_{i}$

By calculus, we can mathemagically get the marginal effect for this model as:

$\frac{\partial y}{\partial x} = β_{1} + 2 * β_{2} x$

So, for our case, the marginal effect of age on wages is given by:

$0.3171 + 2 * - 0.0178 (x - 40) = 0.3171 - 0.0356 (x - 40)$

At age 40 (the zero value), a one year increase in age is associated with a salary increase of $0.32, on average.
For every year over 40, this increase is smaller by $0.0356. For every age younger than 40, this increase is larger by $0.0356.

Finding the inflection point

If the effect of age on wages goes down by $0.0356 for every year over 40, then at some point the positive effect of age on wages will become negative.

We can figure out the value of age at this inflection point by setting the effect to zero and solving for x. In general, this will give us:

$β_{1} / (- 2 * β_{2})$

In our case, we get:

$0.3171 / (- 2 * - 0.0178) = 8.91$

So the model predicts that the effect of age on wages will shift from positive to negative at age 48.91.

Spline models

The basic idea of a spline model is to allow the slope of the relationship between $x$ and $y$ to be different at different cutpoints or “hinge” values of $x$ .
These cutpoints create different linear segments where the effect of x on y is different.
We will look at the case of one hinge value which gives us an overall slope that looks like a “broken arrow.”

Creating the spline variable

The relationship between age and wages suggests that the relationship shifts considerably around age 35. To model this we create a spline variable like so:

$s p l i n e_{i} = {\begin{cases} a g e - 35 & if age>35 \\ 0 & otherwise \end{cases}$

We can then add this variable to our model:

earnings$age.spline <- ifelse(earnings$age<35, 0, earnings$age-35)
model <- lm(wages~age+age.spline, data=earnings)
coef(model)

(Intercept)         age  age.spline 
 -6.0393027   0.9472390  -0.9549087

How do we interpret the results?

Interpreting results in spline model

Up to age 35, the spline term is zero, so our model is given by:

${\hat{w a g e s}}_{i} = - 6.04 + 0.9472 (a g e_{i})$

The model predicts that for individuals age 35 and under, a one year increase in age is associated with a $0.9472 increase in wages, on average.

After age 35, the spline variable increases by one every time age increases by one, so the marginal effect of age is given by:

$0.9472 - 09549 = - 0.0077$

The model predicts that for individuals over age 35, a one year increase in age is associated with a $0.0077 reduction in wages, on average.

Linear model as data-generating process

$y_{i} = \underset{{\hat{y}}_{i}}{\underset{⏟}{β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + \dots + β_{p} x_{i p}}} + ϵ_{i}$

To get values of $y_{i}$ , you feed in values of $x_{i}$ to the structural component and get back out a predicted $\hat{y_{i}}$ value.
To get the stochastic (random) part, you then reach into some distribution to grab a random value of $ϵ_{i}$ that you add to your predicted value to get an actual value of $y_{i}$ .
What distribution are you reaching into when you grab $ϵ_{i}$ ?
- Independence: The number you pull out each time doesn’t depend on other numbers that you pull out. This is most commonly violated by autocorrelation or the clustering of repeated observations.
- Identical distribution: You reach into the same distribution for all observations. This is most commonly violated by heteroscedasticity which means non-constant variance in the residuals.
Together these assumptions give us the IID assumption of the linear model: the error terms are independent and identically distributed.

Heteroscedasticity

Heteroscedasticity means that the variance of the residuals is not constant but depends on the values of $x_{i}$ , and therefore, implicitly, ${\hat{y}}_{i}$ .
A classic example of heteroscedasticity is when the variance of the residuals increases with larger values of ${\hat{y}}_{i}$ giving you a cone shape in a residual by fitted value plot.

Figure 10: The variance of the residuals increases with the predicted value

Using Weighted Least Squares

The weighted least squares technique uses a weighting matrix $W$ in its calculation of regression slopes like so:

$β = (X^{'} W^{- 1} X)^{- 1} X^{'} W^{- 1} y$ $S E_{β} = \sqrt{σ^{2} (X^{'} W^{- 1} X)^{- 1}}$

The exact form of this weighting matrix depends on the nature of the iid violation, but in general it is used to represent the covariance between residuals.

Values in the diagonal cells adjust for heteroscedasticity.
Values in other cells adjust for autocorrelation.

The reality of survey sampling

The simple random sample

In a simple random sample (SRS) of size $n$ , every possible combination of $n$ observations from the population has an equally likely chance of being drawn.
Eevery statistic from an SRS should be representative of the population, except for random sampling bias.

Reality

In practice, large-scale surveys never use SRS for pragmatic and design reasons.
For correct statistical inference, we typically need to make adjustments for sample design.
The primary ways in which sample design affects estimation are clustering, stratification, and weighting.

Weighting

Numerous factors can lead to a sample being unrepresentative such as sampling strata with different probabilities, differential non-response rates, and a lack of fit between sample frame and population.
Sampling weights allow researchers to correct summary statistics from the sample so that they are representative of the population.
The sampling weight for an observation should be $1 / p_{i}$ where $p_{i}$ is the probability of being sampled. The sampling weight indicates the number of observations in the population that observation $i$ in the sample represents.
Calculating sampling weights can be quite complex. In some cases, researchers may know $p_{i}$ from the study design. In other cases, researchers may create post-stratification weights by comparing the sample to some other data source (e.g. census, school records) for a set of demographic characteristics and applying weights to make the sample align with the other data source.
When sampling weights are present in a dataset, they must be used to generate statistics representative of the population.
Variation in sampling weights will increase sampling variability above and beyond that expected for an SRS.

type	income	prob
Star-bellied	3	0.6666667
Star-bellied	7	0.6666667
Star-bellied	6	0.6666667
Regular	2	0.2857143
Regular	1	0.2857143
Regular	4	0.2857143
Regular	5	0.2857143
Regular	0	0.2857143
Regular	2	0.2857143
Regular	1	0.2857143

type	income	prob	weight
Star-bellied	6	0.6666667	1.5
Star-bellied	7	0.6666667	1.5
Regular	1	0.2857143	3.5
Regular	4	0.2857143	3.5

Design issue	Representative?	Design effect (change to SE)
Clustering	Yes, if PSUs are proportionally drawn. Otherwise, weighting necessary.	Increases with difference between clusters.
Stratification	Yes, if strata are sampled with same probability. Otherwise weighting necessary.	Decreases with homogeneity within strata
Weights	Only if weights are applied.	Increases with the variance of the weights.

	unweighted OLS	weighted OLS	robust SE	svy: weights	svy: weights+cluster
Intercept	3.997***	4.027***	4.027***	4.027***	4.027***
	(0.072)	(0.073)	(0.080)	(0.080)	(0.115)
Number of sports played	0.504***	0.533***	0.533***	0.533***	0.533***
	(0.043)	(0.043)	(0.063)	(0.063)	(0.073)
Num.Obs.	4397	4397	4397	4397	4397
R2	0.031	0.033		0.033	0.033
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

	(1)	(2)	(3)
Intercept	3.997***	3.851***	3.391***
	(0.072)	(0.079)	(0.120)
Number of sports	0.504***	0.508***	0.454***
	(0.043)	(0.043)	(0.049)
Drinker		0.707***	0.663***
		(0.155)	(0.180)
Smoker		0.194	0.383*
		(0.162)	(0.188)
Parental income (1000s USD)			0.012***
			(0.002)
Num.Obs.	4397	4343	3332
R2	0.031	0.037	0.048
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

	(1)	(2)	(3)
Intercept	4.051***	3.876***	3.391***
	(0.085)	(0.091)	(0.120)
Number of sports	0.490***	0.494***	0.454***
	(0.049)	(0.049)	(0.049)
Drinker		0.717***	0.663***
		(0.181)	(0.180)
Smoker		0.372*	0.383*
		(0.189)	(0.188)
Parental income (1000s USD)			0.012***
			(0.002)
Num.Obs.	3332	3332	3332
R2	0.029	0.037	0.048
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Sample	$r$ (nominations and income)	SD (income)
Valid cases	0.132	33.7
Valid cases +mean imputed	0.117	29.5
Valid cases +random imputed	0.105	33

	Drop missing	Missing dummy
Intercept	3.931***	3.931***
	(0.107)	(0.106)
Parental income (1000s USD)	0.015***	0.015***
	(0.002)	(0.002)
Parent income imputed		-0.187
		(0.131)
Num.Obs.	3370	4397
R2	0.017	0.014
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

	deletion	mean	mean + dummy	random	regression	random regression	chained equations
Intercept	3.391***	3.360***	3.399***	3.465***	3.367***	3.403***	3.341***
	(0.120)	(0.111)	(0.115)	(0.105)	(0.110)	(0.108)	(0.103)
Number of sports	0.454***	0.476***	0.474***	0.483***	0.461***	0.464***	0.463***
	(0.049)	(0.043)	(0.043)	(0.043)	(0.045)	(0.045)	(0.043)
Drinker	0.663***	0.666***	0.668***	0.666***	0.647***	0.654***	0.614***
	(0.180)	(0.155)	(0.155)	(0.155)	(0.163)	(0.163)	(0.153)
Smoker	0.383*	0.201	0.203	0.202	0.257	0.257	0.221
	(0.188)	(0.161)	(0.161)	(0.161)	(0.170)	(0.170)	(0.159)
Parental income (1000s USD)	0.012***	0.012***	0.012***	0.009***	0.013***	0.012***	0.013***
	(0.002)	(0.002)	(0.002)	(0.002)	(0.002)	(0.002)	(0.002)
Parent income imputed			-0.162
			(0.130)
Num.Obs.	3332	4343	4343	4343	4100	4100	4397
R2	0.048	0.046	0.046	0.044	0.048	0.047	0.049
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Multiple imputation process

Use imputation process with random component to impute missing values and repeat this process to produce $m$ separate complete datasets. Each of these datasets will be somewhat different due to the randomization of imputation. Usually $m = 5$ is sufficient.
Run $m$ separate parallel models on each imputed complete dataset. As a result, you will have $m$ sets of regression coefficients and standard errors.
Pool the regression coefficients across datasets by taking the mean across all $m$ datasets.
Pool standard errors by taking the mean across all $m$ datasets plus the between model standard deviation in coefficients. The formula for the overall standard error is:

$S E_{β} = \sqrt{W + (B + \frac{B}{m})}$ Where $W$ is the squared mean standard error across all $m$ datasets, and $B$ is the variance in coefficient estimates calculated across all $m$ models.

b.pool	se.pool	t.pool	pvalue.pool
3.3762564	0.10682681	31.604954	0.000000e+00
0.0117380	0.00175684	6.681312	2.368128e-11
0.1827329	0.15963841	1.144668	2.523468e-01
0.6447113	0.15501409	4.159050	3.195743e-05
0.4664538	0.04300696	10.846009	0.000000e+00

	(1)	(2)	(3)	(4)	(5)	(6)	(7)
intercept	4.142***	4.135***	10.238***	4.121***	-24.699***	3.078	-29.597***
	(0.172)	(0.231)	(2.052)	(0.226)	(4.941)	(3.272)	(5.309)
poverty rate	0.039***			0.038***	0.127***		0.124***
	(0.007)			(0.010)	(0.017)		(0.017)
unemployment rate		0.091***		0.005		0.100**	0.080*
		(0.024)		(0.033)		(0.036)	(0.033)
median income (logged)			-0.493*		2.516***	0.091	2.912***
			(0.191)		(0.431)	(0.282)	(0.458)
Num.Obs.	341	341	341	341	341	341	341
R2	0.080	0.041	0.019	0.080	0.164	0.041	0.178
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Variance inflation factors

The variance inflation factor (VIF) is the multiplicative factor by which the variance in the estimation of any one coefficient from the current model is increased due to multicollinearity relative to a model in which the variables are uncorrelated. The square root of the VIF is roughly the expected factor increase in the standard error. It can be shown that the VIF for the $i$ th variable is given by:

$V I F_{i} = \frac{1}{1 - R_{i}^{2}}$

Where $R_{i}^{2}$ is the r-squared value when the given independent variable is predicted by all of the other independent variables in the model. For example, we could calculate the VIF for poverty in the full model by:

1/(1-summary(lm(poverty~unemployed+lincome, data = nyc))$r.squared)

[1] 5.984485

The square root of this VIF is 2.45, indicating that the standard error for poverty is more than doubled due to multicollinearity in the full model.

Cronbach’s alpha

Cronbach’s $α$ is a statistic developed in psychology to test the degree to which different variables measure the same underlying concept.

Cronbach’s $α$ is often thought of as a test of the internal reliability of a scale.
You can think about Cronbach’s $α$ as a summary measure of the correlation matrix we saw earlier. It goes from 0 to 1, with 0 indicating no shared correlation, and 1 indicating perfect correlation between all items.

The pysch package includes and alpha command that will calculate $α$ .

psych::alpha(nyc[,c("poverty.z","unemployed.z","lincome.z")])$total

 raw_alpha std.alpha   G6(smc) average_r      S/N         ase          mean
 0.9159293 0.9159293 0.9013009  0.784091 10.89474 0.008313234 -1.373942e-16
        sd  median_r
 0.9252355 0.7419472

The $α$ of 0.916 indicates a very high level of shared covariance between the three variables.

Factor analysis

Another approach to creating a scale is factor analysis. Factor analysis is a method to extract the underlying latent variables (the factors) from a set of observed variables. We will start with an example where we assume a single underlying latent factor.

Lets say we have three variables $z_{1}$ , $z_{2}$ , and $z_{3}$ all measured as z-scores. We can construct a system of equations that relate each of these variables to a common shared factor $F$ (also on a standard scale with a mean of 0 and SD of 1) and three unique factors, $Q_{1}$ , $Q_{2}$ , and $Q_{3}$ .

$z_{i 1} = b_{1} F_{i} + u_{1} Q_{i 1}$ $z_{i 2} = b_{2} F_{i} + u_{2} Q_{i 2}$ $z_{i 3} = b_{3} F_{i} + u_{3} Q_{i 3}$

The $F_{i}$ are called the factor scores and are the values of the common factor for each observation. The $b$ values are called the factor loadings and give the correlation between each observed variable and the common factor. The unique factor components are not actually estimated in factor analysis but serve as the “residual” components.

Factor analysis with multiple factors

The previous example only used one factor and produced results that were very similar to a simple summated scale. However, it is also possible to use factor analysis to identify more than one common factor shared among a set of variables.

The formulas for factor analysis above can be generalized to a set of $J$ observed variables and $m$ factors by a set of $J$ equations. For the $j$ th observed variable:

$z_{j i} = b_{j 1} F_{1 i} + b_{j 2} F_{2 i} + \dots + b_{j m} F_{j i} + u_{j} Q_{j i}$

There will be a set of $J$ factor loadings for each of the $m$ factors. The key question with this technique is what are the appropriate number of factors? This is typically determined by an analysis of:

how much total variation is explained by a given number of factors.
whether the factor loadings for a given number of factors make theoretical sense.

Model Complications Sociology 312/412/512, University of Oregon Aaron Gullickson