1
$\begingroup$

I want to show that the variance of a dependent variable y is a function of the interaction of two independent variables x and z. The value of y may depend on x and z, but I am only interested in how the variance of y is a function of the interaction of x and z. How do I quickly convince my audience that var(y) = f(x*z)?

The best I can do is a quantile regression, and in the example below I use R’s quantreg package. The example below shows that the interaction of x and z widens the interquartile range of y. However, I would prefer to

  1. not rely on the audience’s knowledge of quantile regressions and
  2. show this relation with a plot (hopefully not 3D).

I find plots most convincing, but I do not know how to visualize the relation between the interaction of two variables and the variance of a third variable.

# toy data
library(quantreg)
#> Warning: package 'quantreg' was built under R version 3.5.1
#> Loading required package: SparseM
#> 
#> Attaching package: 'SparseM'
#> The following object is masked from 'package:base':
#> 
#>     backsolve
library(tidyverse)
df <- tibble(
    x = runif(n = 1e4, min = 0.1, max = 5),
    z = runif(n = 1e4, min = 0.1, max = 5),
    y = 5 + rnorm(n = 1e4, mean = 0, sd = 0.1*x*z)
)

# no relation between y and x|z at mean or median
lm1 <- lm(y ~ x + x:z + z, data = df)
summary(lm1)
#> 
#> Call:
#> lm(formula = y ~ x + x:z + z, data = df)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5.0102 -0.2608 -0.0045  0.2742  5.2933 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  4.98083    0.03630 137.202   <2e-16 ***
#> x            0.01371    0.01247   1.100   0.2715    
#> z            0.01816    0.01235   1.470   0.1416    
#> x:z         -0.01133    0.00423  -2.679   0.0074 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.8458 on 9996 degrees of freedom
#> Multiple R-squared:  0.001729,   Adjusted R-squared:  0.00143 
#> F-statistic: 5.772 on 3 and 9996 DF,  p-value: 0.0006123

# but interquartile range of y increases in x*z
rq1 <- rq(y ~ x + x:z + z, tau = (1:9)/10, data = df)
summary(rq1)
#> Warning in summary.rq(xi, U = U, ...): 11 non-positive fis
#> Warning in summary.rq(xi, U = U, ...): 15 non-positive fis

#> Warning in summary.rq(xi, U = U, ...): 15 non-positive fis
#> Warning in summary.rq(xi, U = U, ...): 4 non-positive fis
#> Warning in summary.rq(xi, U = U, ...): 24 non-positive fis
#> Warning in summary.rq(xi, U = U, ...): 70 non-positive fis
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.1
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99615    0.00267 1868.37060    0.00000
#> x              0.00522    0.00373    1.40100    0.16124
#> z              0.00163    0.00277    0.58882    0.55600
#> x:z           -0.13175    0.00350  -37.62846    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.2
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99480    0.00118 4225.96500    0.00000
#> x              0.00727    0.00252    2.88320    0.00394
#> z              0.00329    0.00231    1.42482    0.15424
#> x:z           -0.08901    0.00299  -29.78778    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.3
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99433    0.00203 2457.00881    0.00000
#> x              0.00560    0.00068    8.18339    0.00000
#> z              0.00342    0.00292    1.17092    0.24166
#> x:z           -0.05684    0.00236  -24.06055    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.4
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99796    0.00339 1472.55930    0.00000
#> x              0.00222    0.00291    0.76178    0.44621
#> z              0.00142    0.00288    0.49415    0.62121
#> x:z           -0.02759    0.00286   -9.64994    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.5
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99570    0.00139 3597.40873    0.00000
#> x              0.00403    0.00188    2.14321    0.03212
#> z              0.00527    0.00270    1.95175    0.05100
#> x:z           -0.00435    0.00262   -1.65853    0.09724
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.6
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99490    0.00331 1507.24530    0.00000
#> x              0.00349    0.00278    1.25449    0.20969
#> z              0.00696    0.00254    2.74252    0.00611
#> x:z            0.02066    0.00285    7.24476    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.7
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99694    0.00127 3941.41330    0.00000
#> x              0.00328    0.00261    1.25537    0.20937
#> z              0.00384    0.00088    4.37494    0.00001
#> x:z            0.04930    0.00255   19.34221    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.8
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99574    0.00338 1476.73533    0.00000
#> x              0.00448    0.00326    1.37449    0.16932
#> z              0.00460    0.00295    1.55719    0.11946
#> x:z            0.07849    0.00304   25.80507    0.00000
#> 
#> Call: rq(formula = y ~ x + x:z + z, tau = (1:9)/10, data = df)
#> 
#> tau: [1] 0.9
#> 
#> Coefficients:
#>             Value      Std. Error t value    Pr(>|t|)  
#> (Intercept)    4.99854    0.00098 5113.53938    0.00000
#> x              0.00230    0.00212    1.08323    0.27873
#> z              0.00407    0.00067    6.10957    0.00000
#> x:z            0.12184    0.00324   37.59828    0.00000
plot(rq1)

Update based on the comment from @Repmat about the Bruesch-Pagan test

df <- df %>% 
    mutate(y_mean = mean(y),
                 resid = y - y_mean,
                 resid_sq = resid^2,
                 x_z = x*z)
plot(df$x_z, df$resid)

lm2 <- lm(resid_sq ~ x + z + x:z, data = df)
summary(lm2)
#> 
#> Call:
#> lm(formula = resid_sq ~ x + z + x:z, data = df)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.8838 -0.5245 -0.0225  0.2044 23.8736 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  0.202596   0.073145   2.770  0.00562 ** 
#> x           -0.236146   0.025121  -9.400  < 2e-16 ***
#> z           -0.213769   0.024883  -8.591  < 2e-16 ***
#> x:z          0.253388   0.008522  29.734  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.704 on 9996 degrees of freedom
#> Multiple R-squared:  0.2538, Adjusted R-squared:  0.2536 
#> F-statistic:  1133 on 3 and 9996 DF,  p-value: < 2.2e-16

The coefficient on the interaction between x and z is significantly positive, which is what I want. However, how should I interpret the significantly negative coefficient estimates on x and z? In this context, can I omit the x and z level terms?

Created on 2018-07-05 by the reprex package (v0.2.0).

$\endgroup$
  • 2
    $\begingroup$ How about something like this: 1. For each $y_i$ calculate: $(y_i - \bar{y})^2$ call this $\tilde{y}$ 2. Run the (linear) regression:$\tilde{y}$ on $1, x , z$ This directly mimics the approach of the Breush-Pegan test for heteroscedasticity in linear models. $\endgroup$ – Repmat Jul 5 '18 at 19:20
  • $\begingroup$ @Repmat makes sense. Thanks. In this context, can I omit the x and z level terms? Please see my update above. $\endgroup$ – Richard Herron Jul 5 '18 at 20:39
  • $\begingroup$ It sounds like you're trying to estimate the conditional variance $Var(Y \mid XZ)$. If that's the case, it's not clear to me that you'd want to subtract the mean of $Y$ when implementing @Repmat's suggestion. I think it would make more sense to subtract the conditional mean $E[Y \mid XZ]$, which you could estimate using regression, as in the Breusch-Pagan test. $\endgroup$ – user20160 Jul 7 '18 at 11:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.