4
$\begingroup$

I have a dataset dat where each row represents a soil sample, with independent variables chemical measurements a, b and c, and a dependent variable soil erosion rate value.

Here is a printout of the data:

      a    b     c value
1  1.68 1.84  5.33  1.00
2  1.85 1.75  5.52  1.42
3  1.64 1.48  5.12  1.08
4  2.13 1.68  5.77  1.07
5  2.07 1.68  8.65  2.18
6  1.84 1.74  9.06  2.16
7  2.10 1.97  8.72  1.89
8  2.15 1.91  8.48  2.00
9  2.12 1.20  0.70  1.51
10 1.94 1.19  0.39  2.45
11 2.30 1.38  0.32  1.58
12 2.08 1.35  0.63  1.32
13 2.05 1.65  0.97  1.37
14 1.73 1.65  1.06  1.71
15 2.40 1.86  1.07  1.30
16 2.16 1.95  0.94  1.38
17 2.14 1.18 10.61  3.69
18 2.33 1.38 10.70  3.33
19 2.31 1.28 10.96  3.09
20 2.26 1.08 10.42  2.84
21 2.24 1.74 19.69  3.49
22 2.21 1.45 19.64  4.04
23 2.07 1.96 19.79  3.51
24 1.66 2.07 19.51  3.39
25 1.67 0.83  1.95  1.17
26 1.54 0.69  1.93  1.52
27 1.52 1.02  1.76  0.92
28 1.25 0.88  2.11  0.97
29 2.11 1.99  4.55  2.25
30 2.29 1.50  4.46  1.41
31 2.48 1.65  4.63  2.19
32 2.18 1.51  4.43  2.16

And here is reproducible code to create it in R:

dat = structure(list(a = c(1.68, 1.85, 1.64, 2.13, 2.07, 1.84, 2.1,  
    2.15, 2.12, 1.94, 2.3, 2.08, 2.05, 1.73, 2.4, 2.16, 2.14, 2.33,  
    2.31, 2.26, 2.24, 2.21, 2.07, 1.66, 1.67, 1.54, 1.52, 1.25, 2.11,  
    2.29, 2.48, 2.18), b = c(1.84, 1.75, 1.48, 1.68, 1.68, 1.74,  
    1.97, 1.91, 1.2, 1.19, 1.38, 1.35, 1.65, 1.65, 1.86, 1.95, 1.18,  
    1.38, 1.28, 1.08, 1.74, 1.45, 1.96, 2.07, 0.83, 0.69, 1.02, 0.88,  
    1.99, 1.5, 1.65, 1.51), c = c(5.33, 5.52, 5.12, 5.77, 8.65, 9.06,  
    8.72, 8.48, 0.7, 0.39, 0.32, 0.63, 0.97, 1.06, 1.07, 0.94, 10.61,  
    10.7, 10.96, 10.42, 19.69, 19.64, 19.79, 19.51, 1.95, 1.93, 1.76,  
    2.11, 4.55, 4.46, 4.63, 4.43), value = c(1, 1.42, 1.08, 1.07,  
    2.18, 2.16, 1.89, 2, 1.51, 2.45, 1.58, 1.32, 1.37, 1.71, 1.3,  
    1.38, 3.69, 3.33, 3.09, 2.84, 3.49, 4.04, 3.51, 3.39, 1.17, 1.52,  
    0.92, 0.97, 2.25, 1.41, 2.19, 2.16)), row.names = c(NA, -32L), class = "data.frame")

The variables a, b and c seem to be uncorrelated which I thought means it's OK to include them in the same model:

cor(dat[, c("a", "b", "c")])                                                         
##           a         b         c
## a 1.0000000 0.3323974 0.1891983
## b 0.3323974 1.0000000 0.3353348
## c 0.1891983 0.3353348 1.0000000

My goal is to evaluate the effects of a, b and c on value. According to lm, it seems like all three effects are significant:

fit = lm(value ~ a + b + c, dat) 
summary(fit)                                                                         
## Call:
## lm(formula = value ~ a + b + c, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.88732 -0.22240 -0.05573  0.12671  1.07648 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.08922    0.56288   0.159  0.87519    
## a            1.01036    0.28126   3.592  0.00124 ** 
## b           -0.60997    0.23447  -2.601  0.01467 *  
## c            0.12839    0.01361   9.431 3.45e-10 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.4333 on 28 degrees of freedom
## Multiple R-squared:  0.7995,    Adjusted R-squared:  0.7781 
## F-statistic: 37.23 on 3 and 28 DF,  p-value: 6.611e-10

However, when evaluating the effect of b on its own, the effect is not significant:

fit = lm(value ~ b, dat) 
summary(fit)                                                                         
## Call:
## lm(formula = value ~ b, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1638 -0.7137 -0.2697  0.6349  2.0208 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   1.4815     0.7027   2.108   0.0435 *
## b             0.3708     0.4510   0.822   0.4174  
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.9246 on 30 degrees of freedom
## Multiple R-squared:  0.02204,   Adjusted R-squared:  -0.01056 
## F-statistic: 0.6761 on 1 and 30 DF,  p-value: 0.4174

A scatterplot of b vs. value also shows there is no association between the two:

plot(value ~ b, dat)
abline(lm(value ~ b, dat))

scatterplot of b vs. value

Moreover, the scatterplot suggests a positive effect (if any), while the model value~a+b+c suggests a (significant) negative effect of b!

My questions are:

  1. What should be the conclusion regarding the effect of b on value?
  2. What other tests or diagnostics should be done to investigate this type of situation?

Nota bene: This thread previously had a different title and was marked as a duplicate. Those prior, similar questions focus on how this situation might arise. My question on the practical aspect of reporting the results. Namely, in the described situation, would it be correct to conclude that "b negatively affects value", and if so, what other diagnostics or tests (if any) is it necessary to report.

$\endgroup$
  • 1
    $\begingroup$ Hi, Michael. I think in your correlation coefficient matrix. 'b' is weakly correlated to 'a' and 'c'. It's not uncorrelated. $\endgroup$ – Tbone Jun 5 '20 at 20:40
  • $\begingroup$ @Tbone Thanks, I used the correlation matrix mainly to conclude whether to include all three of the independent variables in the same analysis, using the rule of thumb that all possible pair have weak (<0.7) correlations. $\endgroup$ – Michael Dorman Jun 6 '20 at 8:22
5
$\begingroup$

Let me address some preliminaries first, before addressing the main question. For reference, it will help to read these threads:

  1. Significant predictors become non-significant in multiple logistic regression
  2. How can adding a 2nd IV make the 1st IV significant?
  3. Is there a difference between 'controlling for' and 'ignoring' other variables in multiple regression?
  4. Estimating $b_1x_1+b_2x_2$ instead of $b_1x_1+b_2x_2+b_3x_3$
  5. Intuition behind the names 'partial' and 'marginal' correlations
  6. Basic Simpson's paradox

As those threads discuss, the reason a variable can move from non-significant to significant as more variables are included is that the additional variables are associated with the response and so reduce the error variance, increasing the power of the test of the focal variable. To come to understand how this phenomenon is playing out in your data, you can use the method @whuber demonstrates in his answer to the second linked thread above.

Likewise, the reason the sign flips is because the variable is correlated with the other variables that are being added. I gather seeing the sign flip was perplexing, because you had checked for correlations and found they were small to moderate (you thought they "seem to be uncorrelated"). Nonetheless, the phenomenon is due to those correlations. It's important to recognize that those correlations, and the scatterplot of value x b, are marginal projections (see linked thread #5 above), and relationships can exist between more than two variables oriented in such a way within the full multidimensional space such that they cannot be seen in any of the 2D projections. The best way to see this is to use visualizations that move beyond 'flatland' via conditioning, pseudo-3D representations, motion, or interactivity (e.g., brushing). These can all be done conveniently in R (e.g., via ?coplot, ?lattice, ?scatterplot3d, ?rgl, or ?rggobi). I'm not sure how easy they are in other software, though. A simple hack is to examine a scatterplot matrix and use color, or different symbols, to distinguish different values of your focal variable.

cols = ifelse(dat$b<median(dat$b), "red", "blue")  # low b values are red, high blue
windows()
  pairs(dat[,c(4,1:3)], col=cols, pch=16)
windows()
  plot(value~c, dat, col=cols, pch=16)
  abline(a=1, b=.15, col="gray")

enter image description here

Now it is possible to see if a value for b is low or high. The plots with b aren't the ones to look at, instead, look at how the colors play out in the plots that are collapsing over b. The key plot is the plot of value vs c in the top right corner (row 1, column 4). To get a better look, we can make a larger version of that scatterplot alone:

enter image description here

Here you can see that there is a strong correlation between value and c, with a band of points moving from lower left to upper right. Running parallel within that band, we see stripes of mostly red points (lower values of b) along the band's top and mostly blue points (higher values of b) along the band's bottom (separated by a line I arbitrarily drew through the data). As a result, the relationship between b and value, after controlling for c, is downward sloping. That gets you your negative coefficient sign in the multiple regression model.


Having now attempted to establish the phenomena at play here, and to understand how they are occurring in these data, let's move to the main question: What interpretations can be made, given what you've found?

When we see conflicting results (e.g., non-significant vs. significant, or positive vs negative) when modeled different ways (alone or controlling for covariates), or having made different judgment calls (e.g., which measurement of a construct, transformations, outliers, etc.), it is common to ask, 'which is right?' There isn't necessarily an answer to that. In many cases, transparency is the best course of action: report the analysis planned a-priori, but then discuss other analyses / results as 'sensitivity analyses'.

In this case, I gather the question is about causality ("the effect of b on value"). There are three things to bear in mind:

  1. This appears to have been an observational study. Causal inferences are valid primarily due to the design of the study. That is, you can infer causality because the study is a true experiment due to randomization and independent manipulation of the treatment, or due to exploiting naturally occurring exogeneity.
  2. These effects are due to the correlations amongst variables, and there will likewise be correlations with other variables that are not included in the study. (Variable b is confounded with a and c, which are 'measured confounders', and is certainly confounded with any number of other variables that are unmeasured in this study.)
  3. Each of these models / correlations is a marginal association, and each of the plots is a marginal projection.

Point three is really important to understand. It may help to read linked threads #s 5 and 3 above. It is entirely possible that both models are correct: there is no (or just a very small positive) association with value when ignoring all other possible variables, and there is a clear negative association with value when ignoring all other variables except a and c (which are being controlled for). If you were to gather data on additional variables and assess other marginal associations (controlling for d and e, controlling for a and d, for c and e, for all four covariates, etc.), you could get completely different answers and they could also be right.

I don't necessarily think other tests or diagnostics are needed. The exploratory data visualizations listed above and in @whuber's linked answer can help you with understanding the results of the tests you've already run, though.

Ultimately, if you want to know if b causes value, you need to run a true experiment. Find plots, and divide them randomly into subplots. Independently manipulate the levels of the variables of interest (b, but possibly all three) and treat the subplots. Then wait whatever period of time is appropriate in this context and assess the resulting values. Good experimental design will make a, b, and c orthogonal, and randomization will make all background variables uncorrelated at the population level. That will allow for valid causal inferences.

$\endgroup$
  • $\begingroup$ Thank you very much for the thorough answer, this is more then I could expect! The issue is much more clear to me now. $\endgroup$ – Michael Dorman Jun 5 '20 at 8:43

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.