2.1 What is SEM?

Structural equation modeling is not one thing. Rather, it is a family of techniques for testing multivariate hypotheses about the relationships among observed and (often) latent variables. In general, the goal of SEM is to understand structural patterns in the associations among variables. For example, can twenty items on a depression inventory be thought to measure an underlying (hypothetical) construct?

Critically, SEM allows researchers to test hypotheses about observed variables, as well as latent variables that are not directly measured. This opens the door to develop the construct validity of a latent variable (e.g., depression), or to test hypotheses about how constructs are related to each other (e.g., the relationship between anxiety and depression).

SEM has two important components: structure and measurement.

2.2 Structural model

The structural model represents the relationships among constructs, whether observed or latent. This is closest to standard regression techniques, where believe that one variable causes (or at least is associated with) another.

For example, age often predicts level of performance on working memory tests. Thus, \(Age \rightarrow WorkingMemory\). A developmental psychologist might be quick to note that it may not be age that is causal per se. Rather, the proximate cause might be increasing functional integration of brain networks, or some other outcome of aging. We might think of this in terms of age having an indirect effect on working memory via one or more intermediate (mediating) variables:

#make up data
df <- matrix(rnorm(1000), nrow=100, ncol=4, 
             dimnames = list(NULL, c("Age", "WorkingMemory", "Education", "NetworkIntegration"))
)

#lavaan model syntax
m <- '
WorkingMemory ~ Age
WorkingMemory ~ NetworkIntegration
WorkingMemory ~ Education
Education ~ Age
NetworkIntegration ~ Age
'
mm <- sem(m, data=df)
suppressWarnings(semPaths(mm, style='lisrel', what="diagram", layout = "spring", 
         nCharNodes=20, sizeMan=25, sizeMan2=10, label.scale=FALSE))

Did you notice we just fit our first SEM?! Pretty exciting… I mean, it was technically a path analysis because all variables were observed, but we’ll save that discussion for another day.

2.3 Measurement model

The measurement model specifies how observed variables are related to each other to form underlying latent variables.

That is, observed variables are considered indicators (i.e., indirect measures) of a latent variable. And latent variables typically represent hypothetical constructs in SEM. For example, narcissism is a construct that describes the tendency to be grandiose, vengeful, self-enhancing, and potentially vulnerable to rejection. Research on narcissism (e.g., Pincus & Lukowitsky, 2010) suggests that it likely has multiple facets, such as grandiosity versus vulnerability. So, perhaps it is a hypothetical construct composed of separable latent subcomponents.

The role of a measurement model is to test whether the underlying latent structure of the indicators matches the researcher’s theory. For example, a psychometric item intended to capture grandiosity may not correlate with other grandiosity items. Thus, it would likely be a poor indicator of that latent construct. We’ll spend more time on the notation next week, but this is the basic idea:

grViz("
digraph narc {

  # a 'graph' statement
  graph [overlap = true, fontsize = 12]

  # nodes for observed and latent
  node [shape = circle, fontname = Helvetica]
  'Pathological\nNarcissism'; 'Narcissistic\nGrandiosity'; 'Narcissistic\nVulnerability';
  
  node [shape = box, fontname = Helvetica]
  i1; i2; i3; i4; i5; i6

  # edges
  'Pathological\nNarcissism'->'Narcissistic\nGrandiosity' 'Pathological\nNarcissism'->'Narcissistic\nVulnerability'
  'Narcissistic\nGrandiosity'->'i1' 'Narcissistic\nGrandiosity'->'i2' 'Narcissistic\nGrandiosity'->'i3'
  'Narcissistic\nVulnerability'->'i4' 'Narcissistic\nVulnerability'->'i5' 'Narcissistic\nVulnerability'->'i6'
}")

6.1 Factors affecting correlations

• Linearity: less linear, lower \(r\)
• Variability: smaller \(s\), smaller \(r\)
• Distribution of scores: e.g., opposing skews, smaller \(r\)
• Reliability: lower measurement precision in either variable, lower \(r\)

7.1 Regression diagnostics

This is clearly not a straightforward relationship! The distribution of crime has a long positive tail. We can visualize the residuals of this model, \(\textrm{Price} - \hat{\textrm{Price}}\) to get a sense of things:

ggplot(BostonSmall, aes(x=crim, y=cmedv)) + 
  geom_point() + 
  geom_line(aes(y=pred)) + 
  geom_segment(aes(xend=crim, yend=pred), color="blue", alpha=0.2) +
  theme_bw()

Or the association between the predicted data and residuals:

BostonSmall <- BostonSmall %>% add_residuals(m)
ggplot(BostonSmall, aes(x=pred, y=resid)) + 
  geom_point() + stat_smooth(method="loess")

Or the distribution of residuals (which should be normal):

ggplot(BostonSmall, aes(x=resid)) + geom_histogram(bins=15)

The residuals do not look white and there is evidence of heteroscedasticity. The distribution of crime is very difficult to work with, and we won’t solve it perfectly now.

ggplot(BostonSmall, aes(x=crim)) + geom_histogram(bins=15)

But a log2 transformation of crime moves things in the right direction in terms of normalizing the variance while retaining interpretability (one unit change = 2x crime). Remember that regression doesn’t assume normality of predictors, but instead of the residuals.

ggplot(BostonSmall, aes(x=log2(crim))) + geom_histogram(bins=15)

Refit the model using the log2 transformation. This is moving in the right direction and will have to do for now.

BostonSmall$log_crim <- log2(BostonSmall$crim)
m2 <- lm(cmedv ~ log_crim, BostonSmall)
BostonSmall <- BostonSmall %>% add_predictions(m2)
BostonSmall <- BostonSmall %>% add_residuals(m2)
ggplot(BostonSmall, aes(x=log_crim, y=cmedv)) + 
  geom_point() + 
  geom_line(aes(y=pred)) + 
  geom_segment(aes(xend=log_crim, yend=pred), color="blue", alpha=0.2) + 
  theme_bw()

7.2 Orthogonality of residuals with predictor(s)

Let’s return to the relationships among correlation, covariance, and regression. As Kline notes (p. 27), the residuals of a regression should be orthogonal to the predictor. Indeed, the least squares criterion requires this: \(\underset{\mathbf{X} \in \mathbb{R}}{\textrm{min}} \sum(Y - \hat{Y})^2\).

cor.test(residuals(m2), BostonSmall$log_crim)

## 
##  Pearson's product-moment correlation
## 
## data:  residuals(m2) and BostonSmall$log_crim
## t = 1e-15, df = 500, p-value = 1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0872  0.0872
## sample estimates:
##      cor 
## 4.34e-17

7.3 Swapping criterion and predictor

Unlike correlation or covariance, which are symmetric and undirected, swapping the criterion and predictor in single-predictor regression changes the slope estimate of their association. If we treat crime as the criterion and housing price as the predictor, the slope now represents change in housing prices as a function of change in crime (under log2, a one unit change indicates a 2x increase in crime rate). To wit:

\[ \textrm{log}_2(\hat{\textrm{Crime)}} = B_0 + B_1 \textrm{Price} \]

m3 <- lm(log_crim ~ cmedv, BostonSmall)
summary(m3)

## 
## Call:
## lm(formula = log_crim ~ cmedv, data = BostonSmall)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.951 -2.186 -0.655  2.521  8.600 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.3738     0.3273    7.25  1.6e-12 ***
## cmedv        -0.1553     0.0135  -11.54  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.78 on 504 degrees of freedom
## Multiple R-squared:  0.209,  Adjusted R-squared:  0.208 
## F-statistic:  133 on 1 and 504 DF,  p-value: <2e-16

We can, however, still recover the correlation based on the unstandardized regression coefficient and the standard deviations of the variables:

\[r_{\textrm{Price,log2(Crime)}} = B_1 \frac{s_{\textrm{log2(Crime)}}}{s_{\textrm{Price}}}\]

(coef(m3)["cmedv"] * sd(BostonSmall$cmedv)) / sd(BostonSmall$log_crim)

##  cmedv 
## -0.457

cor.test(~log_crim + cmedv, BostonSmall)

## 
##  Pearson's product-moment correlation
## 
## data:  log_crim and cmedv
## t = -10, df = 500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.524 -0.386
## sample estimates:
##    cor 
## -0.457

7.4 Standardized parameter estimates

This brings us to standardization in regression, which attempts to re-express the parameter estimates in units that are not specific to the scale of each variable. More specifically, standardized parameter estimates represent the association in standard deviation units. In the case of crime and housing prices, we could ask how many standard deviations lower are housing prices given a one SD increase in crime.

Standardization is often achieved by z-scoring all variables (or at least dividing each by its own SD) prior to calculating the regression:

BostonSmall_std <- BostonSmall %>% select(-pred, -resid) %>% mutate_all(scale)
m4 <- lm(cmedv ~ log_crim, BostonSmall_std)
summary(m4)

## 
## Call:
## lm(formula = cmedv ~ log_crim, data = BostonSmall_std)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.885 -0.563 -0.270  0.290  3.627 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.58e-18   3.96e-02     0.0        1    
## log_crim    -4.57e-01   3.96e-02   -11.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.89 on 504 degrees of freedom
## Multiple R-squared:  0.209,  Adjusted R-squared:  0.208 
## F-statistic:  133 on 1 and 504 DF,  p-value: <2e-16

Thus, in the context of single-predictor regression, the standardized regression coefficient is the correlation.

\[r_{xy} := \beta_1 := B_1 \frac{s_y}{s_x} \]

coef(m2)["log_crim"] * sd(BostonSmall$log_crim) / sd(BostonSmall$cmedv)

## log_crim 
##   -0.457

Standardization can also be achieved after model estimation by dividing by the relevant standard deviations. I’ll skip the details for now, but check out the lm.beta package.

lm.beta(m2)

## 
## Call:
## lm(formula = cmedv ~ log_crim, data = BostonSmall)
## 
## Standardized Coefficients::
## (Intercept)    log_crim 
##       0.000      -0.457

7.5 A few other regression equivalencies

As Kline notes, there are a few other handy properties here:

1. The unsigned correlation of predictor and criterion is equal to the unsigned correlation of the criterion with the fitted/predicted values

\[ |r_{XY}| = r_{Y\hat{Y}} \]

More generally, as in multiple regression, this should be expressed:

\[ R_{Y \cdot X_{1...p}} = r_{Y\hat{Y}} \]

cor.test(fitted(m2), BostonSmall$cmedv)

## 
##  Pearson's product-moment correlation
## 
## data:  fitted(m2) and BostonSmall$cmedv
## t = 10, df = 500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.386 0.524
## sample estimates:
##   cor 
## 0.457

2. Total variance in \(Y\) can be decomposed into the predicted part and residual part

\[ s_{Y}^2 = s_{\hat{Y}}^2 + s_{Y - \hat{Y}}^2 \]

var(BostonSmall$cmedv)

## [1] 84.3

var(fitted(m2)) + var(resid(m2))

## [1] 84.3

3. The proportion of shared variance between \(X\) and \(Y\) is the proportion of \(Y\) variance explained relative to its total variance

\[ r_{XY}^2 = \frac{s_{\hat{Y}}^2}{s_{Y}^2} \]

rsquare(m2, BostonSmall)

## [1] 0.209

var(fitted(m2)) / var(BostonSmall$cmedv)

## [1] 0.209

8.1 Standardized parameters in multiple regression

To extend standardized parameter estimates to the multiple regression context, we have to account for the shared variance among predictors.

\[ \beta_\textrm{Nox} = \frac{r_\textrm{Nox,Price} - r_\textrm{Crime,Price} r_\textrm{Nox,Crime}}{1 - r_\textrm{Nox,Crime}^2} \]

So, we adjust the association Nox with Price by the product of the correlations of Crime with Price (other predictor with criterion) and the correlation of Nox and Crime (predictor correlation). The magnitude of the standardized partial coefficient is reduced by adjusting for the correlation among predictors in the denominator. In the case of orthogonal (uncorrelated) predictors, note that the standardized beta simplifies to the bivariate correlation, as in single-predictor regression:

\[ \begin{aligned} r_\textrm{Nox,Crime} :&= 0 \\ \beta_\textrm{Nox} & = r_\textrm{Nox,Price} \end{aligned} \]

Unstandardized betas can be computed from the standardized betas as above:

\[ B_\textrm{Nox} = \beta_\textrm{Nox}\left(\frac{s_\textrm{Price}}{s_\textrm{Nox}}\right) \]

8.2 Partial and part correlations

Also note that we can consider partial relationships in multiple regression in correlational (standardized) terms. But we have to decide how we wish to think about shared variance among predictors. There are two options:

8.2.3 Graphical expression

Venn diagram of shared variance

• \(a\) is the part correlation of \(X_1\) with \(Y\), adjusting for \(X_1-X_2\) association.

• \(\frac{a}{a+d}\) is the partial correlation of \(X_1\) with \(Y\), controlling for \(X_2\).

8.3 Decomposing variance explained in regression (standardized coefficients, part, or partial correlations)

There is not a single best measure of association between a predictor and criterion in multiple regression. Rather, the challenge is that regression coefficients, part correlations, and partial correlations all represent the association of a predictor with a criterion conditional on the other predictors in the model.

The good news is that standardized regression coefficients are interpretable because they are on the same scale and have a known interpretation (relationship between \(Y\) and \(X_1\) at a given level of \(X_2 ... X_p\)). The bad news is that they do not sufficiently weight the direct effect of a predictor on the criterion, which is often of substantive interest (Johson and Lebreton, 2004).

There are additional measures of the relative importance of predictors in a regression model that go beyond what we will cover in class. I would suggest reading this paper and using the relaimpo package to examine relative importance of predictors in regression.

One intuitive, but flawed measure, of partial association is the product of the standardized coefficient with the correlation of the predictor and the criterion (Hoffman, 1960):

\[ r^2_{Y,X_1\cdot (X_2...X_p)} = \beta_{X_1}r_{Y,X_1} \]

By using the marginal correlation (i.e., the bivariate relationship, averaging over all other predictors), the sum of these estimates across all predictors equals the model \(R^2\).

If you’re interested in the relative importance of a predictor in a multiple regression model, I would suggest using the Lindeman, Merenda and Gold (‘lmg’) method of computing relative importance over standardized coefficients, Hoffman’s method, part, or partial correlation. It’s described in the paper above. Here’s a quick example:

m_relimport <- lm(cmedv ~ log_crim + nox, BostonSmall)
stdbetas <- lm.beta(m_relimport)$standardized.coefficients #standardized coefficients

stdbetas["log_crim"] * cor(BostonSmall$cmedv, BostonSmall$log_crim)

## log_crim 
##    0.144

stdbetas["nox"] * cor(BostonSmall$cmedv, BostonSmall$nox)

##   nox 
## 0.078

#here is the broader set of relatience importance measures to consider (esp. 'lmg')
calc.relimp(m_relimport, type = c("lmg", "first", "last", "betasq", "pratt", "genizi", "car"))

## Warning in rev(variances[[p]]) - variances[[p + 1]]: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Response variable: cmedv 
## Total response variance: 84.3 
## Analysis based on 506 observations 
## 
## 2 Regressors: 
## log_crim nox 
## Proportion of variance explained by model: 22.2%
## Metrics are not normalized (rela=FALSE). 
## 
## Relative importance metrics: 
## 
##             lmg   last first betasq pratt genizi    car
## log_crim 0.1232 0.0373 0.209 0.0986 0.144 0.1232 0.1310
## nox      0.0984 0.0125 0.184 0.0330 0.078 0.0984 0.0906
## 
## Average coefficients for different model sizes: 
## 
##              1X     2Xs
## log_crim  -1.35  -0.925
## nox      -34.02 -14.391

As a reminder, the complexity of partial relationships in multiple regression largely comes from the conceptual (and statistical) difficulty of handling the relationship among predictors. If the predictors are uncorrelated with each other, then everything you know from single-predictor regression holds true – for example, the equivalence of bivariate correlation and standardized coefficients.

A useful way to think about this is that if the predictors are truly orthogonal to each other, then if I regress \(\mathbf{y}\) on \(\mathbf{x}_1\) alone, none of the variance accounted for by \(\mathbf{x}_2\) has been removed. Indeed, one could take the residuals of the \(Y ~ X_1\) regression into a new model:

\[ \begin{align} r_{X_1,X_2} :&= 0 \\ \hat{Y} &= B_0 + B_1 X_1 + \varepsilon \\ R_{X_1} &= Y - \hat{Y} \\ \hat{R}_{X_1} &= B_0 + B_2 X_2 + \varepsilon \end{align} \]

The coefficients for \(B_1\) and \(B_2\) will be identical when estimated jointly:

\[ \hat{Y} = B_0 + B_1 X_1 + B_2 X_2 + \varepsilon \]

mu <- rep(0,3)
Sigma <- tribble(
  ~y, ~x1, ~x2,
  1, 0.4, 0.4,
  0.4, 1, 0,
  0.4, 0, 1
)

print(Sigma)

## # A tibble: 3 x 3
##       y    x1    x2
##   <dbl> <dbl> <dbl>
## 1   1     0.4   0.4
## 2   0.4   1     0  
## 3   0.4   0     1

sim_nocorr <- data.frame(MASS::mvrnorm(n=5000, mu=mu, Sigma=Sigma, empirical=TRUE))
names(sim_nocorr) <- c("y", "x1", "x2")
cor(sim_nocorr) #verify that x1 and x2 correlate with y, but not each other

##      y       x1       x2
## y  1.0 4.00e-01 4.00e-01
## x1 0.4 1.00e+00 1.96e-16
## x2 0.4 1.96e-16 1.00e+00

lm.beta(lm(y ~ x1 + x2, sim_nocorr))

## 
## Call:
## lm(formula = y ~ x1 + x2, data = sim_nocorr)
## 
## Standardized Coefficients::
## (Intercept)          x1          x2 
##         0.0         0.4         0.4

lm.beta(lm(y ~ x1, sim_nocorr))

## 
## Call:
## lm(formula = y ~ x1, data = sim_nocorr)
## 
## Standardized Coefficients::
## (Intercept)          x1 
##         0.0         0.4

lm.beta(lm(y ~ x2, sim_nocorr))

## 
## Call:
## lm(formula = y ~ x2, data = sim_nocorr)
## 
## Standardized Coefficients::
## (Intercept)          x2 
##         0.0         0.4

with(sim_nocorr, spcor.test(y, x1, x2))

with(sim_nocorr, spcor.test(y, x2, x1))

with(sim_nocorr, cor(y, x1))

## [1] 0.4

with(sim_nocorr, cor(y, x2))

## [1] 0.4

#demonstrate equivalence of sequential versus join estimation under orthogonal predictors
m1 <- lm(y ~ x1, sim_nocorr)
summary(m1)

## 
## Call:
## lm(formula = y ~ x1, data = sim_nocorr)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.832 -0.620  0.014  0.637  3.076 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.48e-17   1.30e-02     0.0        1    
## x1          4.00e-01   1.30e-02    30.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.917 on 4998 degrees of freedom
## Multiple R-squared:  0.16,   Adjusted R-squared:  0.16 
## F-statistic:  952 on 1 and 4998 DF,  p-value: <2e-16

sim_nocorr$resid_x1 <- resid(m1)
m2 <- lm(resid_x1 ~ x2 - 1, sim_nocorr) #remove intercept since that was pulled out in m1
summary(m2)

## 
## Call:
## lm(formula = resid_x1 ~ x2 - 1, data = sim_nocorr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0906 -0.5686  0.0151  0.5740  2.8157 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## x2   0.4000     0.0117    34.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.825 on 4999 degrees of freedom
## Multiple R-squared:  0.19,   Adjusted R-squared:  0.19 
## F-statistic: 1.18e+03 on 1 and 4999 DF,  p-value: <2e-16

m3 <- lm(y ~ x1 + x2, sim_nocorr)
summary(m3)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = sim_nocorr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0906 -0.5686  0.0151  0.5740  2.8157 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.80e-17   1.17e-02     0.0        1    
## x1           4.00e-01   1.17e-02    34.3   <2e-16 ***
## x2           4.00e-01   1.17e-02    34.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.825 on 4997 degrees of freedom
## Multiple R-squared:  0.32,   Adjusted R-squared:  0.32 
## F-statistic: 1.18e+03 on 2 and 4997 DF,  p-value: <2e-16

Sigma <- tribble(
  ~y, ~x1, ~x2,
  1, 0.4, 0.4,
  0.4, 1, 0.6,
  0.4, 0.6, 1
)
print(Sigma)

## # A tibble: 3 x 3
##       y    x1    x2
##   <dbl> <dbl> <dbl>
## 1   1     0.4   0.4
## 2   0.4   1     0.6
## 3   0.4   0.6   1

sim_withcorr <- data.frame(MASS::mvrnorm(n=5000, mu=mu, Sigma=Sigma, empirical=TRUE))
names(sim_withcorr) <- c("y", "x1", "x2")

#now, semipartial and standardized beta are not equivalent
#likewise, bivariate correlation does not match 
lm.beta(lm(y ~ x1 + x2, sim_withcorr))

## 
## Call:
## lm(formula = y ~ x1 + x2, data = sim_withcorr)
## 
## Standardized Coefficients::
## (Intercept)          x1          x2 
##        0.00        0.25        0.25

lm.beta(lm(y ~ x1, sim_withcorr))

## 
## Call:
## lm(formula = y ~ x1, data = sim_withcorr)
## 
## Standardized Coefficients::
## (Intercept)          x1 
##         0.0         0.4

lm.beta(lm(y ~ x2, sim_withcorr))

## 
## Call:
## lm(formula = y ~ x2, data = sim_withcorr)
## 
## Standardized Coefficients::
## (Intercept)          x2 
##         0.0         0.4

with(sim_withcorr, spcor.test(y, x1, x2))

with(sim_withcorr, spcor.test(y, x2, x1))

with(sim_withcorr, cor(y, x1))

## [1] 0.4

with(sim_withcorr, cor(y, x2))

## [1] 0.4

#demonstrate *non*-equivalence of sequential versus join estimation under correlated predictors
m1 <- lm(y ~ x1, sim_withcorr)
summary(m1)

## 
## Call:
## lm(formula = y ~ x1, data = sim_withcorr)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.795 -0.638 -0.016  0.624  3.325 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.45e-16   1.30e-02     0.0        1    
## x1           4.00e-01   1.30e-02    30.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.917 on 4998 degrees of freedom
## Multiple R-squared:  0.16,   Adjusted R-squared:  0.16 
## F-statistic:  952 on 1 and 4998 DF,  p-value: <2e-16

sim_withcorr$resid_x1 <- resid(m1)
m2 <- lm(resid_x1 ~ x2 - 1, sim_withcorr) #remove intercept since that was pulled out in m1
summary(m2)

## 
## Call:
## lm(formula = resid_x1 ~ x2 - 1, data = sim_withcorr)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.839 -0.626 -0.011  0.617  3.756 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## x2   0.1600     0.0128    12.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.902 on 4999 degrees of freedom
## Multiple R-squared:  0.0305, Adjusted R-squared:  0.0303 
## F-statistic:  157 on 1 and 4999 DF,  p-value: <2e-16

m3 <- lm(y ~ x1 + x2, sim_withcorr)
summary(m3)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = sim_withcorr)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.831 -0.618 -0.016  0.606  3.410 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.59e-16   1.27e-02     0.0        1    
## x1           2.50e-01   1.58e-02    15.8   <2e-16 ***
## x2           2.50e-01   1.58e-02    15.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.895 on 4997 degrees of freedom
## Multiple R-squared:   0.2,   Adjusted R-squared:   0.2 
## F-statistic:  625 on 2 and 4997 DF,  p-value: <2e-16

10.1 Classical suppression

Including a variable that is uncorrelated with the criterion, but correlated with a useful predictor, can inflate the parameter estimate of the useful predictor. For example:

\[ \begin{align} r_{YX} &= 0.5 \\ r_{YZ} &= 0 \\ r_{XZ} &= 0.4 \\ \hat{Y} &= B_0 + B_1 X + B_2 Z + \varepsilon \end{align} \]

In this cirumstance, as we discussed above, the partial regression coefficient adjusts for the associations among predictors.

\[ \beta_\textrm{X} = \frac{r_\textrm{XY} - r_\textrm{YZ} r_\textrm{XZ}}{1 - r_\textrm{XZ}^2} \]

Thus, when there is no association between the criterion \(Y\) and the predictor \(Z\), we get

\[ \beta_\textrm{X} = \frac{r_\textrm{XY}}{1 - r_\textrm{XZ}^2} \] Relative to the bivariate correlation \(r_{XY}\), the denominator will be smaller under classical suppression, leading \(\beta_X\) to be larger.

10.2 Net suppression

Definition: We observe positive bivariate associations among 2+ predictors and a criterion, but one of the predictors has a negative regression coefficient.

\[ \begin{align} r_{YX} &= 0.3 \\ r_{YZ} &= 0.5 \\ r_{XZ} &= 0.4 \\ \hat{Y} &= B_0 + B_1 X + B_2 Z + \varepsilon \end{align} \]

In this case, \(X\) is the predictor having a weaker association with \(Y\) and will tend to take have a negative parameter estimate in the multiple regression, largely reducing error variance in \(Z\), but contributing little to predicting \(Y\). This form of suppression is often identified when the standardized beta for a predictor is opposite in sign to the bivariate correlation of that predictor with the criterion.

10.3 Cooperative suppression

Definition: two predictors correlate negatively with each other, but positively with the criterion. This will lead to a high multiple \(R\) because the predictors will account for more of the criterion when they are considered together than in bivariate correlation. Thus, this form of suppression is often identified when the standardized beta for a predictor exceed its bivariate correlation with the criterion.

\[ \begin{align} r_{YX} &= 0.5 \\ r_{YZ} &= 0.5 \\ r_{XZ} &= -0.4 \\ \hat{Y} &= B_0 + B_1 X + B_2 Z + \varepsilon \end{align} \]

More detail here

11.1 Consider a correlation matrix of items related to subjective well-being, detachment, and narcissism

lower <- '
   1
 .53     1
 .46   .48     1
 .34   .24   .20     1
 .07  -.06  -.02   .50     1
 .24   .13   .15   .53   .60     1
-.02  -.22  -.12   .17   .09   .10     1
-.08  -.27  -.14   .16   .07   .06   .75     1
-.11  -.30  -.13   .12   .03   .00   .72   .80    1
'

narc.cmat <- getCov(lower, 
                    names = c("Restrict_Affect", "Withdrawn", "Avoid_Intimacy",
                              "Grandiose", "Att_Seeking", "Manipulative",
                              "Happy", "Ideal", "Change_Nothing"))

corrplot(narc.cmat, type="lower", diag=FALSE, 
         addCoef.col = "black", method="color", number.cex=1)

11.2 What can SEM bring to this problem?

A crucial point about SEM is that the ‘goodness’ of a model can be thought of in terms of how well it predicts (captures) the full covariance matrix of the data using fewer parameters. So in a covariance matrix, we the number of free parameters, \(p\), is given by:

\[p = \frac{k(k+1)}{2}\]

That is, we have to estimate each cell in the covariance matrix. In the case of dataset above with 9 variables, we have 45 parameters:

\[ p = \frac{9 \cdot 10}{2} = 45\]

This reflects 36 covariances and 9 variances (diagonal). SEM seeks to capture the full covariance matrix \(\mathbf{S_{XX}}\) using a model-implied covariance matrix \(\hat{\mathbf{\Sigma}}\). If this can be done with fewer than 45 parameters, it suggests that there is some structure underlying the covariance matrix that is more parsimonious than a purely unstructured account. This relates to the degrees of freedom in the model, which we will dive into when discussing identification later in the course.

11.3 Testing a model of covariance structure

We can test how well a candidate structural equation model captures the 9 x 9 covariance matrix of these data using fewer than 45 parameters. For example, what if we inferred that there latent traits underlie the three narcissism-related items (Grandiose, Att_Seeking, and Manipulative), three detachment-related items (Restrict_Affect, Withdrawn, Avoid_Intimacy), and three subjective wellbeing items (Happy, Ideal, Change_Nothing)? Furthermore, what if we tested the hypothesis that individuals who are more detached and have greater narcissism report lower subjective well-being?

The first consideration – the presence of latent traits underlying observed responses – represents a measurement question. The second – the associations among latent variables – reflects a structural question. Here is what we would get from testing such a SEM in lavaan. Don’t be concerned with understanding the details at this point!

model_string <- '
#measurement model
Detach =~ Restrict_Affect + Withdrawn + Avoid_Intimacy
Narc =~ Grandiose + Att_Seeking + Manipulative
Subj_Well =~ Happy + Ideal + Change_Nothing

#structural model
Subj_Well ~ Detach + Narc
'

fit_model <- sem(model_string, sample.cov=narc.cmat, sample.nobs=250)
summary(fit_model, standardized=TRUE)

## lavaan 0.6-3 ended normally after 29 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         21
## 
##   Number of observations                           250
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                      73.316
##   Degrees of freedom                                24
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   Detach =~                                                             
##     Restrict_Affct    1.000                               0.700    0.701
##     Withdrawn         1.089    0.135    8.083    0.000    0.762    0.764
##     Avoid_Intimacy    0.904    0.116    7.787    0.000    0.633    0.634
##   Narc =~                                                               
##     Grandiose         1.000                               0.696    0.698
##     Att_Seeking       1.029    0.114    9.033    0.000    0.717    0.718
##     Manipulative      1.145    0.126    9.117    0.000    0.798    0.799
##   Subj_Well =~                                                          
##     Happy             1.000                               0.820    0.821
##     Ideal             1.112    0.066   16.840    0.000    0.911    0.913
##     Change_Nothing    1.067    0.065   16.323    0.000    0.875    0.876
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   Subj_Well ~                                                           
##     Detach           -0.389    0.098   -3.975    0.000   -0.332   -0.332
##     Narc              0.247    0.092    2.680    0.007    0.210    0.210
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   Detach ~~                                                             
##     Narc              0.131    0.043    3.040    0.002    0.268    0.268
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .Restrict_Affct    0.506    0.069    7.341    0.000    0.506    0.508
##    .Withdrawn         0.415    0.071    5.834    0.000    0.415    0.417
##    .Avoid_Intimacy    0.595    0.069    8.631    0.000    0.595    0.598
##    .Grandiose         0.511    0.063    8.117    0.000    0.511    0.513
##    .Att_Seeking       0.482    0.063    7.678    0.000    0.482    0.484
##    .Manipulative      0.360    0.064    5.582    0.000    0.360    0.361
##    .Happy             0.324    0.037    8.741    0.000    0.324    0.326
##    .Ideal             0.166    0.032    5.184    0.000    0.166    0.166
##    .Change_Nothing    0.231    0.033    6.929    0.000    0.231    0.232
##     Detach            0.490    0.093    5.285    0.000    1.000    1.000
##     Narc              0.485    0.088    5.518    0.000    1.000    1.000
##    .Subj_Well         0.593    0.080    7.428    0.000    0.883    0.883

And this is the graphical expression (notation will be covered later!).

semPaths(fit_model, what='std', nCharNodes = 8, sizeMan=11, sizeMan2=4, sizeLat = 13, sizeLat2=7,
         fade=FALSE, fixedStyle = 1, label.cex=1.2, curvePivot=TRUE,
         edge.color='black', esize=5, edge.label.cex=1)

In this model, we’ve represented the full covariance matrix using 22 parameters. Not a bad start, assuming that it captures all 45 covariance values!