Tables and statistics

This week, we’re actually going to be social scientists, and run a study from start to finish: We’ll generate hypotheses, gather data, run a regression model, and write up the results.

Understanding linear regression

This is not a statistics class, but for this class we’re going to dive into it a little bit, running a linear regression, and understanding what it means.

Formulating a hypothesis

I study the news media, and want to understand what factors might influence newspaper circulation. Thinking for a few seconds, I believe that a richer countries have more newspapers, and that people in those countries read more newspapers.

Of course, we have to be very careful not to talk about causation here. We can say that there is a correlation between GDP per capita and newspaper circulation, but we can’t say that one causes the other. Maybe richer countries have more newspapers, or reading a newspaper makes you richer. Or maybe both are caused by a third factor, like education or literacy.

We first write out a hypothesis, which is a statement that we can test with data. In this case, I might say:

Hypothesis 1: Countries with a higher GDP per capita have a higher newspaper circulation.

Gathering our data

Ok, so we have a hypothesis. Now we need to gather some data to test it. We need two things: Newspaper circulation data, and GDP per capita data.

Newspaper circulation

We can get newspaper circulation data from the UN data website. The data is in a zip file, so we need to download it and unzip it.

GDP per capita

We also need GDP per capita data. We can get this from the World Bank. The World Bank has a nice API discussed in chapter 4 that we can use to get the data.

Finally, we use the iso3c code to join the two data frames together. We drop the NA values here, because we don’t want to include countries that don’t have data for both variables. This results in fewer countries in our study than before, but we can just add that to the limitations section, don’t worry.

Running a linear regression with `lm()`

Now, we want to test our hypothesis. We can do this with a linear regression model. Remember way back when, we used the geom_smooth() function to fit a linear regression line to our scatter plot? This is the same thing, but we can get more information out of it.

The lm() function in R is used to fit linear models, and it takes a formula as its first argument. The formula is of the form y ~ x1 + x2 + ..., where y is the dependent variable, and x1, x2, etc. are the independent variables. In this case, we have only one independent variable, gdp_per_capita, so the formula is circulation ~ gdp_per_capita.

When we run the lm() function, it returns an object that contains the results of the regression. We can use the summary() function to get a summary of the results.

Reading a regression table

What are we looking at here? When an outsider is looking at your regression table, they will focus on a few things in particular:

1. The Coefficients table shows the estimated coefficients for each independent variable. In this case, we have one independent

2. Pr(>|t|) columns shows the p-value for each coefficient. The p-value tells us whether the coefficient is statistically

3. The p-value is emphasized by some little stars. The stars are a shorthand for the p-value, and they tell you how statistically

4. The Multiple R-squared value is a measure of how well the model fits the data. It is the proportion of variance in the dependent

The rest of this stuff you should learn from a statistics class. Just know that when reviewers are looking at your scientific paper, these are the things that they’ll zoom in on.

Adding a third variable: Literacy

Ok, so we’ve tested our hypothesis, and we found that there is a statistically significant relationship between GDP per capita and newspaper circulation. But we also know that there are other factors that might influence newspaper circulation, like literacy.

Countries with people that know how to read and write are more likely to have newspapers, and people that can read and write are more likely to read newspapers. So we want to add literacy as a third variable to our model.

First, we need to formulate a new hypothesis:

Hypothesis 2: Countries with a higher literacy rate have a higher newspaper circulation.

Then we need to get the data. We can get literacy data from the World Bank, just like we did with GDP per capita. The indicator we want is SE.ADT.LITR.ZS, which is the literacy rate for people aged 15 and older.

We get the data through our API, and clean it up a little bit.

We then join it back into our main data frame.We should be pros at this by now.

Now, let’s take a look at a univariate regression model with literacy rate as the only independent variable.

Hell yeah! We have a statistically significant relationship between literacy rate and newspaper circulation.

But what happens when we make a multivariate regression model, with both GDP per capita and literacy rate as independent variables?

Noooooo! Literacy rate is no longer statistically significant. Because GDP per capita and literacy rate are correlated, the model is telling us that GDP per capita is a better predictor of newspaper circulation than literacy rate. This is called multicollinearity, and it can be a problem in regression models.

Sometimes, the opposite can happen, as well, where a variable that is not statistically significant in a univariate model becomes significant in a multivariate model.

Because papers that include statistically significant results are more likely to be published, it’s a common temptation to just keep adding variables to your model until you get a statistically significant result. This is called p-hacking.

To avoid this, you should always start with a theoretical model, and then test it with data. Even better is to pre-register your hypotheses and methods before you collect your data.

Adding a categorical predictor

We’re not limited to just continuous variables in our regression model. We can also add categorical variables, like the UN region of a country.

If we look at the regression_df data frame, we can see that we have a variable called un_region. We can use this as a categorical predictor in our regression model.

From this, we can see that there are some differences in newspaper circulation by UN region. However, we can also (annoyingly) notice that one of the UN’s five regions is missing. What’s going on?

The first region is always the reference category, and the other regions are compared to it. This means that the coefficients for the other regions tell us how much higher or lower newspaper circulation is in those regions compared to it.

Because lm() always chooses the first category as the reference category, we can change the reference category by reordering the levels of the factor. For example, but putting the Americas first, we can use that as a baseline.

Note that comparing from a different baseline can affect the significance levels, and this is another way that unethical researchers can hack their results.

Creating a multivariate model

Last, we can create a multivariate model with all of our independent variables. This model will take into account the relationship of GDP per capita, literacy rate, and UN region with newspaper circulation.

Great! Let’s practice this a little.

Classwork: Make and test some hypotheses

1. Think up up with 2 more indicators that you think might be related to a country’s newspaper circulation, and write down a hypothesis

2. Download the data from any source you like.
3. Join the data into your main regression_df data frame.
4. Test them with a bivariate regression model (just the dependent variable and one independent variable).
5. Test them with a multivariate regression model (the dependent variable and all the independent variables).
6. Save this model to your environment.

Here’s an example:

1. I think that old people read more newspapers, so the median age of a country might be related to newspaper circulation. I formulate this

Hypothesis 3: “Countries with a higher median age have a higher newspaper circulation.”

2. I found some data on the UN website

3. I can then read it in with read_csv() and clean it up, then join it to my main data frame.

4. I can then test it first with a bivariate regression model.

5. And then with a multivariate regression model.

That’s science, baby! Send it off to Nature and buy some Prosecco for the apero.

Making good-looking tables

But wait, you might say, do I have to copy and paste all this stuff into a Word document? No, The reason I’ve taught you this was to give you some shortcuts to make your life easier.

A good scientific paper has good-looking tables, and they often follow the same basic structure: A table of summary statistics of all the variables going into the model, a table of regression results, and some stuff in the appendix for more detailed explanations.

Table 1: Summary statistics

First, we need to create a table of summary statistics for all the variables in our model. This is a good way to show the reader what the data looks like, and it can help them understand the results of the regression. This is good to have, because it keeps you from hiding some bad, weird data in your model.

the library gtsummary has a function called tbl_summary() that makes this easy. It takes a data frame as its first argument, and it will automatically calculate the summary statistics for all the variables in the data frame.

We need to tell it which variables we want to include in the table. Because each line of the table is a different country, we want to exclude those from the results. Do this for all the variables in your model.

We also have some options to customize the table. We can change the labels of the variables, and we can change the number of decimal places for the continuous variables.

I don’t expect you to memorize all the options, but you can look them up in the documentation when you need it.

By default, tbl_summary() will show the median, first and third quartiles for continuous variables, and the number of observations and percentage for categorical variables. We can change this with the statistic argument.

Finally, we can use the as_gt() function to convert the table to a gt object, which is the underlying package for making our tables. We can then use the tab_header() function to add a title to the table.

You can do a lot more to customize the table, but I don’t want to overwhelm you with options. You can look at the gt documentation for more information.

We also save it to a variable here.

Finally, we can save the table to a Word document. If you’re not a fan of Word (I’m sure not), you can also save it to a PDF, HTML, LaTeX, or png file. The world is your oyster.

Table 2: Regression results

Now we want to create a table of the regression results. We can use the tbl_regression() function from the gtsummary package to do this. It takes a model object as its first argument, and it will automatically create a table of the regression results.

I’ve added the following options to the table:

1. Just like before, we can customize the table labels with the label argument.
2. We can make significant p-values bold with the bold_p() function. This is a good way to draw attention to the most important results

3. We can make the labels bold with the bold_labels() function.
4. We can italicize the levels of the categorical variables with the italicize_levels() function.
5. We can add a source note to the table with the add_glance_source_note() function. This contains all the other

6. We can use the as_gt() function to convert the table to a gt object, and add a title with the tab_header() function.

When you’re doing this for real, you should look at the documentation and see what you can do.

Table 3: Miscellaneous Data

Great, you’ve added these to your paper, submitted it, and moved on with your life. A few months later, reviewer no. 2 comes back and says that they want to see the full data set added to the appendix.

You can do this with the gt() function, which takes a data frame as its first argument. If you have too much data for one page, you can use the gt_split() function to split the table into multiple pages. Handy!

Automating your results section with `report()`

That handled my tables, but what about the text? I don’t want to have to write all that stuff out by hand. Luckily, I have a solution to that too: the report package.

The report() function takes a model object as its first argument, and it will automatically create a report of the results. It will include the coefficients, p-values, and confidence intervals for each independent variable, as well as the R-squared value for the model. Just take this, and re-word it a little bit so it sounds less robotic, and you have half of your results section written up.

Classwork:

1. Follow the basic steps to create the tables and summaries for your own model.
2. Answer the following questions:
   1. Which independent variables were statistically significant?