ANOVA

Let’s again refer to your knowledge on statistics and sadomania, but what regressionis similar analysis for categorical variables? In fact it is we are trying to compare the differences between a great unifying force and its individual subgroups, which is very similar to analysis of variance (ANOVA - Analysis of Variance). It really is so we can use the data dispersionnogo analysis, using the command anova() on our model.

anova(rmod2)

## Analysis of Variance Table
## 
## Response: earn
##             Df     Sum Sq    Mean Sq F value  Pr(>F)  
## race         3 6.7924e+09 2264121503  2.3241 0.07328 .
## Residuals 1375 1.3395e+12  974196170                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In addition to the ANOVA, there are many similar tests you can use, too:

Let’s look at the dependence of the earnings from the floor. And whether there is a statistically significant difference between the salaries of men and women

smod <- lm(earn ~ sex, data = wages)
coef(smod)

## (Intercept)     sexmale 
##    24245.65    21747.48

wages$sex <- factor(wages$sex, 
levels = c("male", "female"))
smod <- lm(earn ~ sex, data = wages)
coef(smod)

## (Intercept)   sexfemale 
##    45993.13   -21747.48

So, if we take male wages as the base, it will be 45993, and the female will be less on her 21747. The difference was enormous, so let’s check how these data statisticheski reliable:

anova(smod)

## Analysis of Variance Table
## 
## Response: earn
##             Df     Sum Sq    Mean Sq F value    Pr(>F)    
## sex          1 1.5320e+11 1.5320e+11  176.81 < 2.2e-16 ***
## Residuals 1377 1.1931e+12 8.6646e+08                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pr(>F) < 2.2 e-16 *** tells us that the probability of the null hypothesis equal to zero, and thus designed our model a huge difference is statistically significant. Histogram of salaries reassures us that:

qplot(earn, data = wages, geom = "density", color = sex) + theme_bw()

Having constructed a similar chart for the various races we can notice that the difference between races are not so fundamental:

qplot(earn, data = wages, geom = "density", color = race) + theme_bw()

A similar effect of the introduction of the additional variable is widely known and even has its own name - the paradox Simpson

Its essence lies in the fact that the relationship between two variables can change if to consider a third, related variable.

Since this problem is well known, there is a standard solution:

Model each source of variation in the system at the same time

Create a new model that predicts earnings depending on the growth, sex, race, education and age.And try to answer the question is there a relationship between height and income? And between the floor and the earnings?

m3 <- lm(earn ~ height + sex + race + ed + age, data = wages)
coef(m3)

## (Intercept)      height   sexfemale   racewhite   raceblack   raceother 
## -74354.1693    632.7391 -17552.5441   4592.7865   2086.0663   1708.3196 
##          ed         age 
##   4382.0853    287.4744

Symbol . is used in the formula, as a symbol of all the other variables. I.e. we can describe the same model using the code:

lm(earn ~ height + sex + race + ed + age, data = wages)

## 
## Call:
## lm(formula = earn ~ height + sex + race + ed + age, data = wages)
## 
## Coefficients:
## (Intercept)       height    sexfemale    racewhite    raceblack  
##    -74354.2        632.7     -17552.5       4592.8       2086.1  
##   raceother           ed          age  
##      1708.3       4382.1        287.5

lm(earn ~ ., data = wages)

## 
## Call:
## lm(formula = earn ~ ., data = wages)
## 
## Coefficients:
## (Intercept)       height    sexfemale    racewhite    raceblack  
##    -74354.2        632.7     -17552.5       4592.8       2086.1  
##   raceother           ed          age  
##      1708.3       4382.1        287.5

Such a recording, with the use of the mark . you can modify it taking away from the individual variables:

lm(earn ~ height + sex + race + ed, data = wages)

## 
## Call:
## lm(formula = earn ~ height + sex + race + ed, data = wages)
## 
## Coefficients:
## (Intercept)       height    sexfemale    racewhite    raceblack  
##    -46648.9        440.6     -18036.5       6128.5       2915.6  
##   raceother           ed  
##      3379.3       4159.8

lm(earn ~ . - age, data = wages)

## 
## Call:
## lm(formula = earn ~ . - age, data = wages)
## 
## Coefficients:
## (Intercept)       height    sexfemale    racewhite    raceblack  
##    -46648.9        440.6     -18036.5       6128.5       2915.6  
##   raceother           ed  
##      3379.3       4159.8

Important: the order of variables in the formula is of the utmost importance y~x+z and y~z+x are two different models

lm(earn ~ height * sex, data = wages)

## 
## Call:
## lm(formula = earn ~ height * sex, data = wages)
## 
## Coefficients:
##      (Intercept)            height         sexfemale  height:sexfemale  
##         -42677.4            1265.9           30510.4            -701.4

The * character is shorthand notation what we want to consider the impact of these two variables on the dependent variable and want to consider their interaction with each other.

The entry ^2 means that the model needs to take into account all interactions of first order(one-variable - +) and all interactions of second order (interaction of variables with each other - :)

lm(earn ~ height + sex + height:sex, data = w1)

## 
## Call:
## lm(formula = earn ~ height + sex + height:sex, data = w1)
## 
## Coefficients:
##    (Intercept)          height         sexmale  height:sexmale  
##       -31845.8           868.3        139602.2         -1683.4

lm(earn ~ height * sex, data = w1)

## 
## Call:
## lm(formula = earn ~ height * sex, data = w1)
## 
## Coefficients:
##    (Intercept)          height         sexmale  height:sexmale  
##       -31845.8           868.3        139602.2         -1683.4

lm(earn ~ (height + sex)^2, data = w1)

## 
## Call:
## lm(formula = earn ~ (height + sex)^2, data = w1)
## 
## Coefficients:
##    (Intercept)          height         sexmale  height:sexmale  
##       -31845.8           868.3        139602.2         -1683.4

If we want to consider a model with three independent variables,let’s add in our model Russ, we have poyavyatsya interactions and the third order - the height:sex:race

Then we can use the recording ^3, which would indicate that in the model it is necessary to use all interactions of the first, second and third order.

lm(earn ~ (height + sex + race)^3,data = w1)

## 
## Call:
## lm(formula = earn ~ (height + sex + race)^3, data = w1)
## 
## Coefficients:
##                 (Intercept)                       height  
##                     -166395                         2944  
##                     sexmale                 racehispanic  
##                      491439                       183266  
##                   raceother                    racewhite  
##                      385211                       128984  
##              height:sexmale          height:racehispanic  
##                       -7062                        -2437  
##            height:raceother             height:racewhite  
##                       -6064                        -1993  
##        sexmale:racehispanic            sexmale:raceother  
##                          NA                           NA  
##           sexmale:racewhite  height:sexmale:racehispanic  
##                     -396323                           NA  
##    height:sexmale:raceother     height:sexmale:racewhite  
##                          NA                         6052

qplot(height, earn, data = w1, alpha = I(1/10), color = sex) + theme_bw() + geom_line(aes(y = predict(lm(earn ~ (height + sex + race)^3,data = w1))))

Cheat sheet formulas modeling in R

• ‘~’ Separates dependent variables from independent variables from right to left. For example, the formula of dependence y from x, z, and w would look like y ~ x + z + w.
• ‘+’ Separates the independent variables.
• The ‘:’ denotes the interaction between the independent variables. Description y is a value dependent on x, z, and interactions between x and z will look like y ~ x + z + x:z.
• ’’ Short account of all possible interactions between variables. The code y ~ x z * w is equivalent to y ~ x + z + w + x:z + x:w + z:w + x:z:w.
• ‘^’ Denotes the interaction of a certain level. The code y ~ (x + z + w)^2 is equivalent to y ~ x + z + w + x:z + x:w + z:w.
• ‘.’ Brief notation indicating all other variables in the table except for the dependent variable. For example, if the table contains variables x, y, z, and w, then the code y ~ . is equivalent to y ~ x + z + w.
• The ‘-’ sign miinus removes a variable from the formula. For example, y ~ (x + z + w)^2 – x w equivalent to y ~ x + z + w + x:z + z:w.
• ‘-1’ Removes the first coefficient, causing the graph of the model function to pass through the origin .
• ‘I()’ Elements within the braces are interpreted numerically. For example, y ~ x + (z + w)^2 is equivalent to y ~ x + z + w + z:w. But, the code y ~ x + I((z + w)^2) is equivalent to y ~ x + h, where h is a new variable obtained by squaring the sum of z and w.
• ‘function’ is Any mathematical function can be used in the formula. For example, log(y) ~ x + z + w will attempt to describe the dependence of log(y) from x, z, and w.

Analysis of multiple regression models

All of the tasks that stood in front of us to analyze models of linear regressi relevant here. But for models in multiple regression there is another problem - the definition of otdelnoy contribution of each variable in the final model and conclusion about her need to stay in it.

Another problem, besides the fact that the variables that you added may slightly affect your model, is that if present in the model variables correlate, it can contribute to a strong error in their R-values.

In order to avoid such a situation known as multicollinearity, you will have to check the values of correlation for all of your variables

cor(wages$height, wages$ed)

## [1] 0.1140473

cor(wages$height, wages$age)

## [1] -0.1337271

cor(wages$height, as.numeric(wages$sex))

## [1] -0.7036717

If some of your variables are significantly correlated, then you need to choose one of them, abandoning the other.

In addition to these problems, models soerjadi many variables the problem then becomes that an increase in the number of variables increases the chance for some of the variables we polim inadequate R-value is simply the result of a failed sample.

For example, p-value < 0.05 will be every twentieth person among an infinite set of p-values, simply by the law of large numbers.

Then, the probability that at least one p-value < 0.05 = 0.05. The probability that at least one of the two p-values < 0.05 = 0.098. … … … The probability that at least one of the twenty p-values < 0.05 = 0.64

In order not to reject real impact on your model variables as a result of this effect, it is recommended to use one of the following approaches

1. To use a lower α value to determine significance (for example equal to α / p)
2. To check out first what the whole model is statistically significant and therefore make the decision to exclude one of the variables.

For example, for the model m4

summary(m4)

## 
## Call:
## lm(formula = earn ~ height + sex + height:sex, data = wages)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -49699 -20090  -5034  11553 271709 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -42677.4    30488.0  -1.400  0.16180   
## height             1265.9      434.9   2.911  0.00366 **
## sexfemale         30510.4    39644.1   0.770  0.44166   
## height:sexfemale   -701.4      585.8  -1.197  0.23141   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 29340 on 1375 degrees of freedom
## Multiple R-squared:  0.1205, Adjusted R-squared:  0.1186 
## F-statistic: 62.82 on 3 and 1375 DF,  p-value: < 2.2e-16

R-value of the model will be less than 2.2 e-16, which means that your model significantly better predicts the value of a variable than completely arbitrary overkill.

To compare models you can use the functions anova()

anova(m1, m4)

## Analysis of Variance Table
## 
## Model 1: earn ~ height
## Model 2: earn ~ height + sex + height:sex
##   Res.Df        RSS Df  Sum of Sq      F    Pr(>F)    
## 1   1377 1.2318e+12                                   
## 2   1375 1.1840e+12  2 4.7797e+10 27.753 1.528e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The findings indicate that the model m4 has a significant rasnz in the accuracy of predictions than model m1. In other words, because the difference between the models was the presence of the m4 model variable gender(sex), the inclusion of a gender variable in our model improves its accuracy, which means that the variable gender has a right to attend.

If we consider the data obtained from the function anova() for only one model m4, we see that

anova(m4)

## Analysis of Variance Table
## 
## Response: earn
##              Df     Sum Sq    Mean Sq  F value    Pr(>F)    
## height        1 1.1448e+11 1.1448e+11 132.9410 < 2.2e-16 ***
## sex           1 4.6562e+10 4.6562e+10  54.0720 3.307e-13 ***
## height:sex    1 1.2343e+09 1.2343e+09   1.4334    0.2314    
## Residuals  1375 1.1840e+12 8.6112e+08                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value for the variable height:sex more marginal p-znaeniya, and therefore the model contains the variable heights + sex + heights:sex is not significantly better than the model containing the variables heights + sex, meaning of a variable height sex, we can refuse.

The final algorithm of the research model will be to obtain and interpret the following values

• The coefficients coef(mod)
• Balances - resid(mod)
• The confidence interval confint(mod)
• P-values - summary(mod)
• Comparison of variants model - anova(mod1, mod2)