2
$\begingroup$

I have a data set that includes 6 quantitative variables. In those variables, I have determined that each of LYMPHO and NEUTRO have statistically significant correlations with WBC. So the question I'm not understand is, "is it necessary to include both LYMPHO and NEUTRO as predictor variables in the regression model if we treat WBC as the response variable?"

I know there is a certain concept where if you put a certain variable in first (i.e. Length) then there isn't much left for another variable (i.e. Circumference) to explain for a given response variable (i.e. Weight). Would that be applicable here as well?

cor.test(WBC, LYMPHO)
# this gives an R = 0.8269 and an R^2 = 0.6838
cor.test(WBC, NEUTRO) 
# this gives a R = 0.5443 and an R^2 = 0.2962

Therefore can explain 68.3% of the variation in WBC by LYMPHO and 29.6% of the variation in WBC by NEUTRO. Does NEUTRO even need to be included then?

Sample of the data set:

    HAEMO PCV  WBC LYMPHO NEUTRO LEAD
1    13.4  39 4100     14     25   17
2    14.6  46 5000     15     30   20
3    13.5  42 4500     19     21   18
4    15.0  46 4600     23     16   18
5    14.6  44 5100     17     31   19
6    14.0  44 4900     20     24   19
7    16.4  49 4300     21     17   18
8    14.8  44 4400     16     26   29
9    15.2  46 4100     27     13   27
10   15.5  48 8400     34     42   36   
$\endgroup$

1 Answer 1

1
$\begingroup$

You can use R stats::anova two ways:

Make up an lm model where you include the following (inclusion order important): covariate mostly correlated with output; next covariate that also strongly correlated with output. You can play with this manually how you want. Fit the model lm_model, and use anova(lm_model).

dat <- data.frame(y = rnorm(100, 0, 1), x1 = rnorm(100, 0, 1), x2 = rnorm(100, 0, 1))

lm_model <- lm(y ~ x1 + x2, data = dat)

summary(lm_model)

anova(lm_model)

> anova(lm_model)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value  Pr(>F)  
x1         1  0.026  0.0260  0.0288 0.86550  
x2         1  3.828  3.8283  4.2485 0.04196 *
Residuals 97 87.407  0.9011                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

If your model consists of more than 1 coefficient, the anova will make estimation of:

  • whether x1 adds significant residual elimination given intercept term

    whether x2 adds significant residual elimination given intercept and x1.

You can just look at Pr(>F) respective to the variable added.

The second way is to feed two models to anova:

lm_model_1 <- lm(y ~ x1, data = dat)

lm_model_2 <- lm(y ~ x1 + x2, data = dat)


anova(lm_model_1, lm_model_2)


   > anova(lm_model_1, lm_model_2)
Analysis of Variance Table

Model 1: y ~ x1
Model 2: y ~ x1 + x2
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     98 91.235                              
2     97 87.407  1    3.8283 4.2485 0.04196 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

You get the same p-value for the second term x2 after looking at both approaches to the same thing. If the p-value is lower than precpecified alhpa (significance level) you are fine with inclusion of both predictors.

Also consider this:

lm_model_0 <- lm(y ~ 1, data = dat)

lm_model_1 <- lm(y ~ x1, data = dat)

lm_model_2 <- lm(y ~ x1 + x2, data = dat)


anova(lm_model_0, lm_model_1, lm_model_2)


> anova(lm_model_0, lm_model_1, lm_model_2)
Analysis of Variance Table

Model 1: y ~ 1
Model 2: y ~ x1
Model 3: y ~ x1 + x2
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     99 91.261                              
2     98 91.235  1    0.0260 0.0288 0.86550  
3     97 87.407  1    3.8283 4.2485 0.04196 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  • model_0 is based on a constant term only;
  • model_1 includes one factor besides the constant;
  • model_2 includes one more factor besides what was inside, so it is nested.

If the same result apperas. If you just use model_1 and model_2, then contrasting to a (unseen) model with a constant term is included in anova by default unless you specify the model without intercept.

And this:

> lm_model <- lm(y ~ x1 + x2 - 1, data = dat)
> 
> #summary(lm_model)
> 
> anova(lm_model)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value  Pr(>F)  
x1         1  0.017  0.0166  0.0186 0.89183  
x2         1  3.893  3.8926  4.3480 0.03965 *
Residuals 98 87.735  0.8953                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I removed intercept, and results changed.

$\endgroup$
3
  • $\begingroup$ Alexey Burnakov+ Thank you for the reply! That is really helpful. I'm just curious how come you used: dat <- data.frame(y = rnorm(100, 0, 1), x1 = rnorm(100, 0, 1), x2 = rnorm(100, 0, 1)) $\endgroup$ Oct 6, 2017 at 1:37
  • $\begingroup$ This is a dummy dataset created purely for code exampling. It will not give you any significance in y ~ x, so if you look at the F-statistics value, and the probability that given true null-hypothesis the observed F will be at least so big, Pr(>F), you see its not very close to zero. You can use Your data and play with permutating first, second, and maybe third covariates in the formula. If you found my reply really helpful please mark as the reply or ask for more clarification. ) $\endgroup$ Oct 6, 2017 at 9:05
  • $\begingroup$ I modified my answer slightly so that you could see that both approaches produce the same. $\endgroup$ Oct 6, 2017 at 9:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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