Correlation

R can perform correlation with the cor() function. Built-in to the base distribution of the program are three routines; for Pearson, Kendal and Spearman Rank correlations.

The pseudo code to get the correlation coefficient would be:

cor(var1, var2, method = "method")

The default method is “pearson” so you may omit this if that is what you want. If you type “kendall” or “spearman” then you will get the appropriate correlation coefficient.

Correlation coefficients - The default correlation returns the pearson correlation coefficient: cor(var1, var2) - If you specify “spearman” you will get the spearman correlation coefficient: cor(var1, var2, method = "spearman") - If you use a datset instead of separate variables you will return a matrix of all the pairwize correlation coefficients: cor(dataset, method = "pearson")

Getting a correlation coefficient is generally only half the story; you will want to know if the relationship is significant. The cor() function in R can be extended to provide the significance testing required. The function is cor.test()

To run a correlation test, the pseudo code would look like this:

cor.test(var1, var2, method = "method")

The default method is “pearson” so you may omit this if that is what you want. If you type “kendall” or “spearman” then you will get the appropriate significance test.

As usual with R it is a good idea to assign a variable name to your result in case you want to perfom additional operations.

Correlation Significance tests

• The default method is “pearson” cor.p = cor.test(var1, var2)
• If you specify “spearman” you will get the spearman correlation coefficient: cor.s = cor.test(var1, var2, method = "spearman")
• To see a summary of your correlation test type the name of the variable e.g. cor.s

So if we want to test the correlation between life expectancy and GDP per capita, using a pearson correlation, then we would do the following:

gapminder_corr <- cor.test(gapminder_2007$lifeExp, gapminder_2007$gdpPercap)
gapminder_corr

    Pearson's product-moment correlation

data:  gapminder_2007$lifeExp and gapminder_2007$gdpPercap
t = 10.933, df = 140, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5786217 0.7585843
sample estimates:
      cor 
0.6786624 

You can also plot the correlation coefficient on top of a scatterplot of the variables:

ggplot(gapminder_2007, aes(x = lifeExp, y = gdpPercap)) + geom_point() + geom_smooth(colour = "red", 
    fill = "lightgreen", method = "lm")

Finally, you can get a correlation matrix, to assess the relationship between all your variables in your data set. For this, we need to select only the quantitative variables in our data set. We should be experts at such subsetting by now, so let’s do this by selecting the quantitative variables:

quant_gap <- gapminder_2007 %>% select(lifeExp, pop, gdpPercap)

Now print a correlation table of all the quantitative variables:

res <- cor(quant_gap)
round(res, 2)

          lifeExp   pop gdpPercap
lifeExp      1.00  0.05      0.68
pop          0.05  1.00     -0.06
gdpPercap    0.68 -0.06      1.00

You can also use the corrplot package to visualise this:

install.packages("corrplot")

library(corrplot)

corrplot 0.84 loaded

corrplot(cor(quant_gap), type = "upper", order = "hclust", tl.col = "black", 
    tl.srt = 45)

The function corrplot() takes the correlation matrix as the first argument. The second argument (type="upper") is used to display only the upper triangular of the correlation matrix.

Chi-square

So you want to look at a cross-tab between two variables. This is as easy as using the table() function, and passing it the two variables as arguments. For example, if we want to look at the relationship between our two dichotomous variables we created about high and low gdp and also below average life expectancy or not, we can just type:

table(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp)

          
           high life expectancy low life expectancy
  high gdp                   44                   3
  low gdp                    41                  54

Excellent, now we can run a chi square test with the… chisq.test() function. Inside the brackets you have to pass a frequency table. Like the one we made below. This can be something you saved as an object, or you can just recreate the table inside the function.

Very straightforward. Now what we want to do is save this to an object:

life_exp_v_gdp <- chisq.test(table(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp))

Now, we can call all sort of information from this chi square test object:

Your test statistic:

life_exp_v_gdp$statistic

X-squared 
 31.25234 

Your degrees of freedom:

life_exp_v_gdp$parameter

df 
 1 

Your p value:

life_exp_v_gdp$p.value

[1] 2.265738e-08

Your expected counts for the cells:

life_exp_v_gdp$expected

          
           high life expectancy low life expectancy
  high gdp              28.1338             18.8662
  low gdp               56.8662             38.1338

Your residuals:

life_exp_v_gdp$residuals

          
           high life expectancy low life expectancy
  high gdp             2.991291           -3.652840
  low gdp             -2.104000            2.569318

Your standardized residuals:

life_exp_v_gdp$stdres

          
           high life expectancy low life expectancy
  high gdp             5.772285           -5.772285
  low gdp             -5.772285            5.772285

However if you want to see it all together, you can use the CrossTable() function from the gmodels package:

install.packages("gmodels")

You can even specify for your table to follow the format you might be familiar with from SPSS or SAS. Use ?CrossTable() to see all the parameters you can specify with this function:

library(gmodels)
CrossTable(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp, prop.t = FALSE, 
    prop.r = FALSE, prop.c = FALSE, expected = TRUE, resid = TRUE, sresid = TRUE, 
    format = "SPSS")

   Cell Contents
|-------------------------|
|                   Count |
|         Expected Values |
| Chi-square contribution |
|                Residual |
|            Std Residual |
|-------------------------|

Total Observations in Table:  142 

                          | gapminder_2007$below_avg_life_exp 
gapminder_2007$gdp_factor | high life expectancy  |  low life expectancy  |            Row Total | 
--------------------------|----------------------|----------------------|----------------------|
                 high gdp |                  44  |                   3  |                  47  | 
                          |              28.134  |              18.866  |                      | 
                          |               8.948  |              13.343  |                      | 
                          |              15.866  |             -15.866  |                      | 
                          |               2.991  |              -3.653  |                      | 
--------------------------|----------------------|----------------------|----------------------|
                  low gdp |                  41  |                  54  |                  95  | 
                          |              56.866  |              38.134  |                      | 
                          |               4.427  |               6.601  |                      | 
                          |             -15.866  |              15.866  |                      | 
                          |              -2.104  |               2.569  |                      | 
--------------------------|----------------------|----------------------|----------------------|
             Column Total |                  85  |                  57  |                 142  | 
--------------------------|----------------------|----------------------|----------------------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  33.31927     d.f. =  1     p =  7.820385e-09 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  31.25234     d.f. =  1     p =  2.265738e-08 

 
       Minimum expected frequency: 18.8662 

Fisher test

Now since we are all stats-minded people, you may have noticed that we have a very low cell count in one of our cells in our contingency table. Since this can call into doubt our findings, it’s always helpful to use tests that account for these issues as well. In this case, the test would be a Fischer exact test, which you carry out in R with the fisher-test() function. This function also just needs a contingency table passed to it as an argument:

fisher.test(table(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp))

    Fisher's Exact Test for Count Data

data:  table(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp)
p-value = 1.457e-09
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   5.479616 102.155242
sample estimates:
odds ratio 
  18.95314 

Odds ratio

To get the odds ratio from a 2 x 2 cross tab, you can use the oddsratio() function from the vcd package.

install.packages("vcd")

And again, all you need to pass it, is a contingency table, and an arument to specify that you want odds ratio, and not log odds

library(vcd)

Loading required package: grid

oddsratio(table(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp), 
    log = FALSE)

 odds ratios for  and  

[1] 19.31707

You can see how it’s blank where it says “and”. If you passes it a table with column names, it would tell you what your odds ratio is being calculated for. But luckily with the code right there, we still see, and if we were to forge the values, we can always just re-print the table. In this case, we know now that the odds of having below-average life expectancy for low gdp countries are 19 times the odds for high gdp countries.

But we know that we had issues with the small cell count, and in general in stats we like to be explicit about our uncertainties, so let’s also include some confidence intervals here. So a final step is to do this with the confint() function:

confint(oddsratio(table(gapminder_2007$gdp_factor, gapminder_2007$below_avg_life_exp), 
    log = FALSE))

                                                             2.5 %
high gdp:low gdp/high life expectancy:low life expectancy 5.601292
                                                            97.5 %
high gdp:low gdp/high life expectancy:low life expectancy 66.61844

T-test

For categorical variables with two possible values, we would use the t-test, if the numeric variable meets the assumption of normal distribution. By now you should be able to guess what functions are called in R! Any ideas for the t-test?

Well… the function is called t.test(). The R function t.test() can be used to perform both one and two sample t-tests on vectors of data. The function contains a variety of options and can be called as follows:

t.test(x ~ y, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, 
    var.equal = FALSE, conf.level = 0.95)

Here x is a numeric vector of data values and y is a binary factor. The option mu provides a number indicating the true value of the mean (or difference in means if you are performing a two sample test) under the null hypothesis. The option alternative is a character string specifying the alternative hypothesis, and must be one of the following: “two.sided” (which is the default), “greater” or “less” depending on whether the alternative hypothesis is that the mean is different than, greater than or less than mu, respectively.

For example the following call:

t.test(x, alternative = "less", mu = 10)

performs a one sample t-test on the data contained in x where the null hypothesis is that =10 and the alternative is that <10. The option paired indicates whether or not you want a paired t-test (TRUE = yes and FALSE = no). If you leave this option out it defaults to FALSE. The option var.equal is a logical variable indicating whether or not to assume the two variances as being equal when performing a two-sample t-test. If TRUE then the pooled variance is used to estimate the variance otherwise the Welch (or Satterthwaite) approximation to the degrees of freedom is used. If you leave this option out it defaults to FALSE. Finally, the option conf.level determines the confidence level of the reported confidence interval for in the one-sample case and 1- 2 in the two-sample case.

So if we want to look at life expectancy between high and low gdp countries, using a t-test, we would:

t.test(gapminder_2007$lifeExp ~ gapminder_2007$gdp_factor, alternative = "two.sided", 
    paired = FALSE, conf.level = 0.95)

    Welch Two Sample t-test

data:  gapminder_2007$lifeExp by gapminder_2007$gdp_factor
t = 9.4568, df = 133.89, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 11.40431 17.43610
sample estimates:
mean in group high gdp  mean in group low gdp 
              76.65474               62.23454 

So now you can say that the difference in mean life expectancy between high and low gdp countries is statistically significant, and is between 11 - 17 years.

ANOVA

And what about the difference in life expectancy between continents? The function here is called aov() and the setup is very similar to the t.test():

fit <- aov(x ~ y, data = mydataframe)

Where x refers to the numeric variable, and y to the categorical variable with multiple possible values. You also specify the data souce, so you only have to use the column names as the x and the y values. For other variations of ANOVA, eg two factor design, see this quick start guide

fit <- aov(lifeExp ~ continent, data = gapminder_2007)
summary(fit)

             Df Sum Sq Mean Sq F value Pr(>F)    
continent     4  13061    3265   59.71 <2e-16 ***
Residuals   137   7491      55                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

You also want to carry out post-hoc testing:

TukeyHSD(fit)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = lifeExp ~ continent, data = gapminder_2007)

$continent
                      diff        lwr       upr     p adj
Americas-Africa  18.802082  13.827042 23.777121 0.0000000
Asia-Africa      15.922446  11.372842 20.472051 0.0000000
Europe-Africa    22.842562  18.155858 27.529265 0.0000000
Oceania-Africa   25.913462  11.183451 40.643472 0.0000305
Asia-Americas    -2.879635  -8.299767  2.540497 0.5844901
Europe-Americas   4.040480  -1.495233  9.576193 0.2630707
Oceania-Americas  7.111380  -7.910342 22.133102 0.6862848
Europe-Asia       6.920115   1.763372 12.076858 0.0027416
Oceania-Asia      9.991015  -4.895221 24.877251 0.3464970
Oceania-Europe    3.070900 -11.857807 17.999607 0.9793776

Density estimate of the residuals

We could plot the density estimate of the residuals:

plot(density(rstudent(fit_2)))

We can see that the right hand side tail is heavier than it should be.

outlierTest(fit_2)

   rstudent unadjusted p-value Bonferonni p
1 -3.721986         0.00019766     0.028068