E is for Effect Sizes

Today is the first in a three-part series this month on how to conduct meta-analysis using R, plus a fourth Statistics Sunday post tying it all together. As I mentioned in yesterday's post, we'll be using the metafor package. If you didn't install the package then, you'll want to install it now.

install.packages("metafor")

## Installing package into '\\marge/users$/slocatelli/My Documents/R/win-library/3.4'
## (as 'lib' is unspecified)

library(metafor)

## Loading required package: Matrix

## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

For this post, we'll focus on computing effect sizes from the summary data. If this sounds like Greek to you, or if you need a refresher, you'll want to first review this post on meta-analysis and this post introducing effect sizes. But to briefly summarize, to conduct a meta-analysis, you're taking information from individual studies to convert into an effect size. You collect these from multiple studies on the same topic to try to estimate the true effect size - the value that all of these individual studies, using different methods, samples, and (yes) flaws are trying to estimate. Later on, we'll go into how you can aggregate those individual study effect sizes.

Let's assume, since all studies are looking at the same topic, that they report the same kind of summary data that we convert into an effect size. In practice, this isn't always true, which is when you have to get into converting between effect sizes. But that's an advanced topic.

Meta-analyses frequently use one of three types of effect sizes: correlation, standardized mean difference (which requires mean and standard deviations for the 2 groups or time points being compared), or odds ratio (or similar metric used for binary data). In some medical meta-analysis, you might also look at event counts over time, such as number of strokes occurring in a set period of time. If you're meta-analyzing correlations, studies often give you your effect size directly, though there are corrections you may want to apply. But if you're meta-analyzing mean differences, ratios, or counts, you have to pull in the summary statistics and convert those to your effect size.

The metafor package will do these calculations easily, using the escalc (effect size calculation) function. You'll state in the function what kind of measure you want, and this will determine what data is needed.

Standardized Mean Difference

If you're computing a standardized mean difference, abbreviated as SMD in metafor, you'll need your study data file to include: mean for each group, standard deviation for each group, and sample size per group. In grad school, I did a meta-analysis on pretrial publicity. Though most studies used guilty/not guilty verdicts as their study outcome, a handful used guilt ratings. Just for fun, I pulled some of those studies out of my study dataset. Here's the data you need to create a data frame we can analyze in metafor:

smd_meta<-data.frame(
  id = c("005","005","029","031","038","041","041","058","058","067","067"),
  study = c(1,2,3,1,1,1,2,1,2,1,2),
  author_year = c("Ruva 2007","Ruva 2007","Chrzanowski 2006","Studebaker 2000",
                  "Ruva 2008","Bradshaw 2007","Bradshaw 2007","Wilson 1998",
                  "Wilson 1998","Locatelli 2011","Locatelli 2011"),
  n1 = c(138,140,144,21,54,78,92,31,29,90,181),
  n2 = c(138,142,234,21,52,20,18,15,13,29,53),
  m1 = c(5.29,5.05,1.97,5.95,5.07,6.22,5.47,6.13,5.69,4.81,4.83),
  m2 = c(4.08,3.89,2.45,3.67,3.96,5.75,4.89,3.80,3.61,4.61,4.51),
  sd1 = c(1.65,1.50,1.08,1.02,1.65,2.53,2.31,2.51,2.51,1.20,1.19),
  sd2 = c(1.67,1.61,1.22,1.20,1.76,2.17,2.59,2.68,2.78,1.39,1.34)
)

ID is a number I assigned to keep track of every source I examined for the meta-analysis, including ones I didn't end up using. Study refers to the study number within the source, since some sources had multiple studies. I used the study number the authors gave, so I could easily refer back to the source if necessary. Not all studies qualified for the meta-analysis, so you'll notice some numbers are skipped.

You may also notice that study ID 067 is mine. Since the meta-analysis, originally done in 2009, was on the same topic as my dissertation, which I completed in 2011, I reran the meta-analysis for one of my dissertation chapters, adding in the data from my own dissertation. In the data above, the group 1s are the treatment group, who saw pretrial publicity, and group 2s are the control group, who did not. This is important to remember when it comes time to interpret the direction of the effect sizes. A positive value means the treatment group gave a higher guilt rating. While each study may use different scales for the guilt rating, since we're standardizing, these differences don't matter. The mean differences are in standard deviation units.

We now apply the escalc function for standardized mean difference, referencing the smd_meta data frame so it appends effect sizes (as a column called "yi") to the data:

smd_meta <- escalc(measure="SMD", m1i=m1, m2i=m2, sd1i=sd1, sd2i=sd2, n1i=n1, n2i=n2,
                   data=smd_meta)

Odds Ratio

Next, we can conduct the same kind of effect size calculation for odds ratio, which is what I used for most of the pretrial publicity studies. As I said, most used verdicts, so each study had a 2x2 table of results:

Once again, here's a data frame pulled from my original meta-analysis dataset:

or_meta<-data.frame(
  id = c("001","003","005","005","011","016","025","025","035","039","045","064","064"),
  study = c(1,5,1,2,1,1,1,2,1,1,1,1,2),
  author_year = c("Bruschke 1999","Finkelstein 1995","Ruva 2007","Ruva 2007",
                  "Freedman 1996","Keelen 1979","Davis 1986","Davis 1986",
                  "Padawer-Singer 1974","Eimermann 1971","Jacquin 2001",
                  "Ruva 2006","Ruva 2006"),
  tg = c(58,26,67,90,36,37,17,17,47,15,133,68,53),
  cg = c(49,39,22,50,12,33,19,17,33,11,207,29,44),
  tn = c(72,60,138,140,99,120,60,55,60,40,136,87,74),
  cn = c(62,90,138,142,54,120,52,57,60,44,228,83,73)
)

I provided guilty counts by group (tg = treatment guilty verdicts, cg = control guilty verdicts), as well as sample sizes per group, which I can use to get not guilty verdicts. We can now request odds ratios (note: metafor gives log odds ratios) with the escalc function:

or_meta <- escalc(measure="OR", ai=tg, bi=(tn-tg), ci=cg, di=(cn-cg), data=or_meta)

But I could request other types of binary effect sizes, such as risk ratios (once again, log-transformed automatically), or risk difference. The metafor package gives more information on the different measures you can request, how they're interpreted, and some sample datasets using the different measures.

In part 2, we'll talk about variances, and in part 3, weights, with a Statistics Sunday post talking about using the results from escalc to generate an aggregate effect size, one of the main goals of meta-analysis. Check back later this month for those posts! And tomorrow, check back for a post on conducting confirmatory factor analysis in R, once again using the Facebook dataset.

D is for Data Frame

Working in R involves dealing with various objects, such as functions and a variety of data structures. One of the objects I work with the most in R is the data frame, a data table in which rows are cases and columns are variables (in the research, rather than programming, sense of the word - while I tend to fall into research speak, column is probably the better word to use in this context). Unlike some other data structures in R, the columns can be of different types - for instance, one column may contain string (text) data, another numeric, another Boolean indicator (TRUE, FALSE), and so on. R and its various libraries comes with many built-in data frames. And creating one with your own data is very easy and can be accomplished in multiple ways. A data frame is then used in statistical analysis. Today, I'll show you some of the main methods for creating data frames from scratch or reading them in from other sources.

Creating a Data Frame From Scratch

Creating a data frame from scratch involves binding together other R objects, such as matrices, arrays, or vectors. This makes sense if you have raw data not yet entered into some other type of data object or, as I've done in some past R posts, when you're generating data. Each row name must be unique, and the same goes for each column, though the name could be as simple as a number. For instance, in my alpha and beta posts, I referred to columns in my Facebook data frame by their column number: e.g., Facebook[,3], which references column 3. (Note: You can use that same notation to refer to row numbers; just put that information within the brackets and before the comma. For instance, Facebook[3,] references row 3. Leaving the space before the comma blank includes all rows, and leaving the space after the comma blank includes all columns, in the data frame.) I could have instead referred by column name, with the dataset$variable format, but number is easier if you're referencing a range of columns.

You can assign whatever names you want - for instance, in my Facebook data frame, column names were pulled from the header in the tab-delimited file - but any duplicates will give you an error. By default, R assigns row names as its number in the data frame, but I could also set that myself; for instance, I might want to assign my unique ID variable as my row names:

Facebook<-read.delim(file="small_facebook_set.txt", header=TRUE)
row.names(Facebook)<-Facebook$ID

Another important thing to remember is that each column in an R data frame must have the same number of values (rows). Those values can be missing, but there still has to be something there, or else R will give you an error when you try to create a data frame.

test<-data.frame(
  meas_id = c(1:6),
  score = c(1,2,5,4,5)
)

## Error in data.frame(meas_id = c(1:6), score = c(1, 2, 5, 4, 5)): arguments imply differing number of rows: 6, 5

So let's start with the first way to create a data frame from scratch - manually type in values. This approach makes sense when the number of rows and columns is small. Let's create a data frame listing the measures I used in my Facebook study, which I'm using as an ongoing example this month:

measures<-data.frame(
  meas_id = c(1:6),
  name = c("Ruminative Response Scale","Savoring Beliefs Inventory",
           "Satisfaction with Life Scale","Ten-Item Personality Measure",
           "Cohen-Hoberman Inventory of Physical Symptoms",
           "Center for Epidemiologic Studies Depression Scale"),
  num_items = c(22,24,5,10,32,16),
  rev_items = c(FALSE, TRUE, FALSE, TRUE, FALSE, TRUE)
)

This code creates a data frame of 6 rows and 4 columns. I've assigned a measurement id, and given the name of the scale, the number of items, and an indicator of whether the scale has any reversed items. Also, I'd like to point out that I use equal signs with the data.frame code. This not only creates that column, it gives the column that name. Using arrows <- instead gives weird column names. If you're curious - hey, the best way I learn is through a combination of curiosity and making (sometimes purposeful) mistakes - change = to <- and see what happens.

Each column is a different data type. I can ask R what the data type is for a specific variable with str followed by the dataset$variable in parentheses, or for all variables with str(dataset):

str(measures)

## 'data.frame': 6 obs. of  4 variables:
##  $ meas_id  : int  1 2 3 4 5 6
##  $ name     : Factor w/ 6 levels "Center for Epidemiologic Studies Depression Scale",..: 3 5 4 6 2 1
##  $ num_items: num  22 24 5 10 32 16
##  $ rev_items: logi  FALSE TRUE FALSE TRUE FALSE TRUE

Int refers to integer, num to number, and logi to Boolean or logical. Strangely, R thinks the measurement name column is a factor, which isn't what I wanted. I can fix this by forcing R to make the column a character vector:

measures$name<-as.character(measures$name)
str(measures$name)

##  chr [1:6] "Ruminative Response Scale" "Savoring Beliefs Inventory" ...

On the other hand, if I'm generating random data, I would write the code to create my data then bind them together into a data frame, again enclosing that code with data.frame(). For examples on how I've done this, see here, here, and here.

Reading Data into R

It's more likely that you already have data in some other form and want to read it into R. I'm a big fan of tab-delimited and CSV files, both of which are easy to create, small and compact even with a lot of data, and really easy to read in. I recommend including variable names in the first row of the file, though you can add them later if they're not there. Reading in a tab-delimited file is easy with the read.delim function, which automatically generates an R data frame - just remember to name it so it will store it in the R environment:

mytabdata <- read.delim("file.txt", header=TRUE)

If you're working with a CSV file, it has it's own function, read.csv:

mycsvdata <- read.csv("file.csv", header=TRUE)

If your files are not saved in the working directory, put the entire path within the file quotes - or better yet, change the working directory with setwd("path"). The defaults for these two functions are tab-separator and comma-separator, respectively, so you don't have to specify, though you could with sep="/t" and sep=",", respectively. These two also default to character data enclosed in quote-marks "". 

If you're accessing a tab-delimited or CSV file saved at a url, just put url(webaddress) where I currently have "file.txt" or "file.csv". There are also ways to read HTML tables into data frames, but I think that's another post for another day.

And if your file does not contain variables in the first line, leave out header=TRUE (or set to header=FALSE). If it does, but you'd rather not use those, you can tell R to skip those by adding skip=# (number of rows to skip). Then add in column names as a vector, with each name enclosed in quote marks and separated by commas - make sure you provide as many names as there are columns or you'll get an error. For instance, if I had instead read in the measure information from above and needed to provide column names separately:

colnames(measures)<-c("meas_id","name","num_items","rev_items")

You can also load a pre-existing dataset in R, and this is frequently done to show examples of what R can do. If a dataset is part of base R, you can load with data(name). A popular example dataset in R is mtcars - a dataset taken from the 1974 Motor Trend US magazine. Sometimes, the data(name) code loads the set but doesn't make it a dataframe. You can quickly coerce it into a dataframe by requesting the head (column names plus first several rows):

data(mtcars)
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Applying a function to the data will also usually coerce it into a data frame. But the head(data) function is a quick way to do that, and is nice if you want to view the data. 

On the other hand, if a dataset is part of an R library, you need to load that library first. For instance, the metafor package I'll be using tomorrow, which is used to conduct meta-analysis, comes with many built-in datasets. The McDaniel 1994 dataset contains studies examining the validity of employment interviews. You can pull it up in the same way as mtcars, after loading (install first if you haven't before) the metafor package:

install.packages("metafor")

## Installing package into '\\marge/users$/slocatelli/My Documents/R/win-library/3.4'
## (as 'lib' is unspecified)

library(metafor)

## Loading required package: Matrix

## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

data(dat.mcdaniel1994)
head(dat.mcdaniel1994)

##   study   ni   ri type struct
## 1     1  123 0.00    j      s
## 2     2   95 0.06    p      u
## 3     3   69 0.36    j      s
## 4     4 1832 0.15    j      s
## 5     5   78 0.14    j      s
## 6     6  329 0.06    j      s

To keep this post from getting obscenely long, I'll wrap up here, but there are a few other ways to read files into data frames, which I hope to blog about later. If you're working with SPSS, SAS, Stata, or other data from statistical packages, you can read them in with the Hmisc (SPSS and SAS) or foreign (all of the above, though not as powerful as Hmisc) package. If you're working directly with Excel files, and don't want to convert to tab-delimited or CSV first, I recommend the XLConnect package, which also lets you move R objects, including graphics, over to Excel.

Finally, you can read in other file types of fixed-width files, XML, or JSON - which have a variety of applications in large-scale testing, particularly computer-adaptive testing. I worked with fixed-width files quite a bit at HMH, since data for some of our cognitive ability tests was sent to us in fixed-width files, and a colleague there frequently worked with JSON, which was used for some of our large-scale computer-adaptive and computer-administered tests. Look for future blog posts on those (and what they are, if these are new terms for you), since there are some nuances to the code that require a bit more depth.

C is for Cross-Tabs Analysis

Cross-tabs is short for cross-tabulation (or cross tables), and are used to display frequencies for combinations of categorical and sometimes ordinal variables. In a previous blog post, I described chi-square, which is frequently used to analyze cross-tabs. As I said in that post, chi-square is a little like an ANOVA in how it examines the relationship between the two variables, but it is used to look at an association between two variables with two (and sometimes more) levels. Basically, you would use a non-parametric test like chi-square when working with variables that aren't continuous or that violate key assumptions of parametric tests.

Today, I'll demonstrate how to generate cross-tabs in R (which I also did in the previous post) and I'll show two ways to analyze cross-tabs: chi-square (again) and a similar test, Fisher's exact test.

Once again, I'll use my Facebook dataset to demonstrate. None of the variables in the set I used previously would really qualify for cross-tabs, since the variables in that dataset are meant to be combined into continuous scales. I included gender, but because of the student body makeup at the school where I collected my data, there are a lot of women and very few men - which isn't optimal for this type of analysis. But I can pull in some additional data collected in the demographics portion of the survey. The goal of the study was to examine Facebook usage patterns, so participants reported whether they used Facebook. The vast majority of the sample did. But they also reported use of other social media, including Twitter and LinkedIn.

A bit more than half the sample were freshmen, since the participant pool is drawn from Introductory Psychology, mostly taken by freshmen. The remaining participants were upper-class or non-traditional students. It makes sense that these two groups might have different usage patterns. That is, older students might be more likely to use LinkedIn, as they prepare for job searches and networking; non-traditional students are likely to already have a job or career. It's unclear, however, whether we might see a similar difference for Twitter users - though keep in mind, these data were collected in 2010, when Twitter may have been a different landscape. So let's generate two cross-tabs, both using the freshmen versus upper-class/non-traditional students (or younger versus older, for simplicity), one to look at Twitter use and one to look at LinkedIn use.

First, I'll read in that data, then redefine the variables I need as factors, which includes age2 (recoded from the continuous age variable) and indicators for using Twitter and LinkedIn. This gives labels to my cross-tabs. In order to make changes to these variables, I have to refer to these variables in my code, first to reflect that I want to change that variable (the information before the <-) and again when I reference what variable to make a factor. I use the dataset$variablename syntax to refer to a specific variable:

age_socialmedia<-read.delim(file="age_usage.txt", header=TRUE)
age_socialmedia$age2<-factor(age_socialmedia$age2, labels=c("Younger","Older"))
age_socialmedia$Twitter<-factor(age_socialmedia$Twitter, labels=c("Non-User","User"))
age_socialmedia$LinkedIn<-factor(age_socialmedia$LinkedIn, labels=c("Non-User","User"))

If I had wanted, I could have created a new variable in the first part of the code. If I wrote the name of a variable that doesn't exist in the dataset, it would be added. But since I'm not recoding, just adding labels, I have no issue with overwriting the existing variable.

I can generate my tables with the following:

Twitter<-table(age_socialmedia$age2, age_socialmedia$Twitter)
Twitter

##          
##           Non-User User
##   Younger      111   29
##   Older         91   25

LinkedIn<-table(age_socialmedia$age2, age_socialmedia$LinkedIn)
LinkedIn

##          
##           Non-User User
##   Younger      139    1
##   Older        109    7

Use of either social media site is not very high in this sample, but much too low for LinkedIn to use chi-square - one of the assumptions of that test is that no cells have counts less than 5. Fisher's exact test will work in that situation, though, so we can use chi-square for Twitter and Fisher's exact test for LinkedIn.

The code for either is very easy, especially if you named your tables:

chisq.test(Twitter)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Twitter
## X-squared = 9.2489e-05, df = 1, p-value = 0.9923

The Yates' continuity correction was developed because chi-square is biased to be significant when samples are large. We can easily turn this feature off, and most people do:

chisq.test(Twitter, correct=FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  Twitter
## X-squared = 0.026729, df = 1, p-value = 0.8701

In this case, the correction made little difference - yes, the p-value is smaller but neither is even close to being significant. So we'd conclude that, in this sample, there is no age difference in Twitter use.

Now let's conduct our Fisher's exact test, which as I mention above, can be used when there are cells with counts less than 5:

fisher.test(LinkedIn)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  LinkedIn
## p-value = 0.02477
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##    1.112105 404.350276
## sample estimates:
## odds ratio 
##   8.863863

This test is significant; LinkedIn users tend to be older in this sample. In fact, older students are over 8 times more likely to be LinkedIn users than younger students in the present sample.

B is for Betas (Standardized Regression Coefficients)

Welcome to Day 2 of Blogging A to Z! As with yesterday, the title of today's post is very similar to Day 2 of Blogging A to Z last year, but once again, a different concept.

Regression is used to predict scores on one variable with one or more predictor variables. As I mentioned previously, regression is similar to correlation, which describes (numerically) the strength of the relationship between two variables. But the goal of regression is not just to describe a relationship but to create an equation that allows one to predict scores, or specifically to see if information from the x variable(s) can be used to generate close approximations of the y variable. It might not be used for prediction in the sense of forecasting future events, though that is one way regression results can be used. The conversion of x variable(s) into y variable is accomplished with regression coefficients - values that are multiplied by the value of x to generate a predicted y, which is hopefully very similar to the observed y. A constant is included in a regression equation to help shift the scale - this constant is equal to the mean value for y when x is equal to 0. When people talk about regression, they usually mean linear regression, which is what I'll focus on today. But there are other types of regression for nonlinear relationships between y and x(s).

In a previous post I predicted my rating of books I read last year with a linear regression that included book length, genre, author gender, and how long it took to read the book in days. In that analysis, I found that book length and fantasy genre predicted higher ratings and YA fiction predicted lower ratings. The other variables were not significant.

Going back to the Facebook file I used yesterday, I have many scales I could use in a linear regression. In that study, among other scales, participants reported rumination on the Ruminative Response Scale (RRS) and depression on the Center for Epidemiologic Studies Depression Scale (CES-D). The relationship between rumination and depression is well-established, and though this is a non-clinical sample, we would expect to find a relationship, such that heightened scores on the RRS should predict higher scores on the CES-D. In fact, here's the scatterplot showing the relationship between RRS and CES-D in this sample; as you can see, it's a positive linear relationship:

We could run a simple linear regression with these two variables. First, we need to generate our scale scores, since at the moment, the file only contains responses to individual items.

In a previous post, I noted that some items are reverse-scored. There are multiple ways I could go about reverse-scoring items and generating scores. Since one of those ways involves a package I plan to discuss more later, for the time being, I'll just write some of my own code to reverse-score. Since I'll be doing that with multiple variables, I'll write a custom function I can reuse.

reverse<-function(max,min,x) {
  y<-(max+min)-x
  return(y)
  }

I can then apply this function to the 3 CES-D items that need to be reverse-scored, providing the max and min rating values, as well as the value (x) I wanted reverse-scored, and create a new variable with "R" added to indicate it is reversed:

Facebook<-read.delim(file="small_facebook_set.txt", header=TRUE)
Facebook$Dep4R<-reverse(3,0,Facebook$Dep4)
Facebook$Dep8R<-reverse(3,0,Facebook$Dep8)
Facebook$Dep12R<-reverse(3,0,Facebook$Dep12)

Now I'll generate my scores. RRS doesn't have any reversed items, so I can just add those columns together; it does, however, have three subscales: Depression-Related Rumination (fixating on one's negative traits or feelings), Brooding (negative thoughts more generally), and Reflection (attempting to understand oneself and one's mood, which could be a positive experience). The CES-D has no subscales.

Facebook$RRS<-rowSums(Facebook[,3:24])
Facebook$RRS_D<-rowSums(Facebook[,c(3,4,5,6,8,10,11,16,19,20,21,24)])
Facebook$RRS_R<-rowSums(Facebook[,c(9,13,14,22,23)])
Facebook$RRS_B<-rowSums(Facebook[,c(7,12,15,17,18)])

Facebook$CESD<-rowSums(Facebook[,c(96,97,98,100,101,102,104,105,106,108,109,110,111,112,
                                   113,114)])

If you don't like scientific notation on your p-values (I don't), be sure to change those options before displaying results; I usually add this code at the beginning of every R session:

options(scipen=999)

I can use the RRS variables in my regression, though I wouldn't want to use the total RRS score in the same regression as the three subscales, since the subscales are derived from the total score. First, let's run a very simple regression with RRS total score and CES-D, which we can do with the lm (for "linear model") function:

RumDep<-lm(CESD~RRS, data=Facebook)
summary(RumDep)

## 
## Call:
## lm(formula = CESD ~ RRS, data = Facebook)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.3020  -3.3885  -0.7835   2.4140  17.2783 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  8.45024    0.86992   9.714 <0.0000000000000002 ***
## RRS          0.19753    0.02132   9.264 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.107 on 255 degrees of freedom
## Multiple R-squared:  0.2518, Adjusted R-squared:  0.2489 
## F-statistic: 85.82 on 1 and 255 DF,  p-value: < 0.00000000000000022

Our regression results are significant, and this is a strong relationship; our R-squared (explained variance) is 0.25. But let's see what happens if we differentiate between the 3 types of rumination.

RumDep2<-lm(CESD~RRS_D+RRS_R+RRS_B, data=Facebook)
summary(RumDep2)

## 
## Call:
## lm(formula = CESD ~ RRS_D + RRS_R + RRS_B, data = Facebook)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.944  -3.308  -0.677   2.572  18.271 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  8.12981    0.86644   9.383 < 0.0000000000000002 ***
## RRS_D        0.36845    0.06312   5.838         0.0000000162 ***
## RRS_R        0.04613    0.09928   0.465                0.643    
## RRS_B       -0.05401    0.12766  -0.423                0.673    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.045 on 253 degrees of freedom
## Multiple R-squared:  0.2757, Adjusted R-squared:  0.2671 
## F-statistic:  32.1 on 3 and 253 DF,  p-value: < 0.00000000000000022

So the relationship between rumination and depression is driven specifically by self-focused rumination, rather than morose thoughts or simply reflecting on one's feelings. Differentiating between types of rumination also gives us a slightly higher R-squared. You may look at the regression coefficients (in the column called "Estimate") and notice that RRS_D is much larger than RRS_R or RRS_B. But remember these three subscales have different numbers of items and are therefore on different scales; I can't really compare them directly. In fact, you may have regression equations containing many predictors on different scales. What if more than one had been significant? How could I compare them?

I can standardize them. Standardized regression coefficients, or betas, are in standard deviation units or Z-scores. This allows you to directly compare your coefficients, because they describe the number of standard deviation units y will change for a 1 standard deviation change in x. To access standardized regression coefficients, we'll need to use another R package called QuantPsyc. Then we use a function in that package called lm.beta to access standardized coefficients for the linear model we ran. (This is why it's always good to name those results with name<- followed by the function; that way we can refer to them again later.)

library(QuantPsyc)

## Warning: package 'QuantPsyc' was built under R version 3.4.4

## Loading required package: boot

## Loading required package: MASS

## 
## Attaching package: 'QuantPsyc'

## The following object is masked from 'package:base':
## 
##     norm

lm.beta(RumDep2)

##       RRS_D       RRS_R       RRS_B 
##  0.53245571  0.03262223 -0.03693183

Because units are standardized (Z-scores), the mean is 0, so we don't have a constant in this equation. These results tell us if our score on the Depression-Related Rumination is equal to 1 standard deviation, our CES-D score will equal 0.53 standard deviation units. In fact, by standardizing the coefficients, we see how much stronger the relationship between Depression-Related Rumination and Depression is than Reflecting and Depression or Brooding and Depression - 16.3 times and 14.4 times stronger, respectively. 

Just for fun, let's see what happens if we add in some additional variables, specifically personality, which was measured by a brief Big-Five measure. Unlike RRS and CES-D, which are scored by summing, the authors of the Ten Item Personality Measure have you average together the two items used to measure each of the 5 traits, after reverse-scoring half of the items. 

Facebook$CritR<-reverse(7,1,Facebook$Critical)
Facebook$AnxR<-reverse(7,1,Facebook$Anxious)
Facebook$ResR<-reverse(7,1,Facebook$Reserved)
Facebook$DisR<-reverse(7,1,Facebook$Disorganized)
Facebook$ConvR<-reverse(7,1,Facebook$Conventional)

Facebook$Extraversion<-(Facebook$Extraverted+Facebook$ResR)/2
Facebook$Agree<-(Facebook$CritR+Facebook$Sympathetic)/2
Facebook$Consc<-(Facebook$Dependable+Facebook$DisR)/2
Facebook$EmoSt<-(Facebook$AnxR+Facebook$Calm)/2
Facebook$Openness<-(Facebook$NewExperiences+Facebook$ConvR)/2

For simplicity, I'll just use total RRS score in this regression.

PersonalityDep<-lm(CESD~RRS+Extraversion+Agree+Consc+EmoSt+Openness, data=Facebook)
summary(PersonalityDep)

## 
## Call:
## lm(formula = CESD ~ RRS + Extraversion + Agree + Consc + EmoSt + 
##     Openness, data = Facebook)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.4740  -3.4643  -0.7196   2.5938  17.7445 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  12.09086    3.20496   3.773             0.000202 ***
## RRS           0.19269    0.02148   8.972 < 0.0000000000000002 ***
## Extraversion -0.65750    0.43047  -1.527             0.127929    
## Agree         0.25896    0.38931   0.665             0.506557    
## Consc         0.49505    0.38714   1.279             0.202182    
## EmoSt         0.12558    0.39897   0.315             0.753197    
## Openness     -1.03245    0.38435  -2.686             0.007710 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.062 on 250 degrees of freedom
## Multiple R-squared:  0.2793, Adjusted R-squared:  0.262 
## F-statistic: 16.15 on 6 and 250 DF,  p-value: 0.000000000000001076

Other than rumination, only openness to experience significantly predicts depression scores, in this case negatively - people with higher scores on openness to experience have lower depression scores. Fortunately, the 5 personality variables are on the same scale, but they're on a different scale than the RRS. Let's request standardized coefficients so we can compare all of our predictor variables.

lm.beta(PersonalityDep)

##          RRS Extraversion        Agree        Consc        EmoSt 
##   0.48951737  -0.08595408   0.03780856   0.07145380   0.01756582 
##     Openness 
##  -0.14954644

The impact of personality variables were quite small, particularly compared to rumination. The effect of rumination on depression is 3.3 times stronger than the strongest personality trait, openness to experience. When I dropped rumination, only two predictors were significant: extraversion and openness to experience, and the relationships still weren't all that strong. Also, the R-squared was small - less than 0.05.

PersonalityDep2<-lm(CESD~Extraversion+Agree+Consc+EmoSt+Openness, data=Facebook)
summary(PersonalityDep2)

## 
## Call:
## lm(formula = CESD ~ Extraversion + Agree + Consc + EmoSt + Openness, 
##     data = Facebook)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.967  -3.746  -1.208   2.777  19.995 
## 
## Coefficients:
##              Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)   17.2529     3.6179   4.769 0.00000314 ***
## Extraversion  -1.0544     0.4913  -2.146     0.0328 *  
## Agree          0.6601     0.4438   1.487     0.1382    
## Consc          0.7030     0.4434   1.585     0.1141    
## EmoSt          0.4017     0.4564   0.880     0.3796    
## Openness      -1.0318     0.4410  -2.340     0.0201 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.809 on 251 degrees of freedom
## Multiple R-squared:  0.04726, Adjusted R-squared:  0.02828 
## F-statistic:  2.49 on 5 and 251 DF,  p-value: 0.03185

Take that, Google dude.

A is for (Cronbach's) Alpha

You may recall that last April, I blogged the A to Z of Statistics, and the very first post was alpha, or the type I error rate. This year, I'm blogging the A to Z of R, and this post is about a different kind of alpha - Cronbach's alpha, a measure of reliability, specifically internal consistency. You can find more detail about Cronbach's alpha conceptually by checking out today's Statistics Sunday post. This post focuses on how to compute Cronbach's alpha in R.

You can find the alpha function in the psych library - a package I'll be blogging about more later this month. Install if you haven't yet and always load the library before you get started (even if you just installed it):

install.packages("psych")

library(psych)

To compute alpha, you'll need a set of items (at least 3) that measure the same concept. Determining if items measure the same concept is a validity issue, but it doesn't really make sense to examine reliability when you have no evidence that the items measure the same thing. For my example throughout much of this month, I'll use the dataset from my Facebook study, since the dataset includes 6 published measures that have previous evidence of reliability and validity. I created a "small" version of the data - one that only includes ID, gender, and responses to the 7 measures.

Because this is only for demonstration purposes, I don't really want to share my original dataset - I'd rather, if people are interested in a copy, that they reach out to me and request it officially - but I did create a simulated dataset, which has similar means and correlations as the original data, and is formatted in the same way. You can download that dataset here to analyze along with me, as well as a mini codebook here, but note that your values will differ slightly. Make sure that when you save it to your computer that you either put it in your R working directory or that you change your working directory to the folder you saved it in.

First, I want to read in the data:

Facebook<-read.delim("small_facebook_set.txt", header=TRUE)

I can start examining my individual measures in the data. To help with analysis, I noted which columns corresponded to each measure, and which items are reversed:

• Ruminative Response Scale, columns 3-24
• Savoring Beliefs Inventory, columns 25-48, even-number items are reversed (e.g., 2nd item, 4th item)
• Satisfaction with Life Scale, columns 49-53
• Ten-Item Personality Measure, columns 54-63
• Cohen-Hoberman Inventory of Physical Symptoms, columns 64-95
• Center for Epidemiologic Studies Depression Scale, columns 96-111, items 4, 8, and 12 are reversed

If you downloaded the simulated dataset referenced above, that file has the same columns and reversed items.

The code for Cronbach's alpha is pretty simple; you include the range of columns that are part of the measure and which items are reverse-scored with the keys statement. For instance:

Savoring<-alpha(Facebook[,25:48], keys=c(seq(from = 2, to = 24, by = 2)))

That command runs Cronbach's alpha on that range of columns - the 24 items of the Savoring Beliefs Inventory - and reverse scores the even-numbered items. The seq command produces the same result as if I had typed keys=c(2,4,6,8,10,12,14,16,18,20,24). I've now created an R object called Savoring that contains the results of Cronbach's alpha. Requesting a summary of that object will give me the raw and standardized alpha, another measure of Cronbach's alpha, the average correlation between items, and a few other descriptive statistics:

summary(Savoring)

Reliability analysis   
 raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd
      0.93      0.93    0.94      0.35  13 0.0066  5.5 0.82

I can get more detailed output (check out a sample here) by requesting the full object (Savoring, without the summary part), which tells me what would happen to alpha if I dropped a particular item and provides statistics for each item (e.g., mean and standard deviation, proportion selecting each response option, and so on).

If you're analyzing a scale with no reversed items, you'd just drop the keys statement:

Rumination<-alpha(Facebook[,3:24])

Not long ago, I blogged that I had misplaced my study codebook, which, among other things, identified which items are reversed. Fortunately, the alpha function has an option that could help if I had been unable to recover that information. The check.keys statement will conduct a principal components analysis of the scale and identify items with negative loadings - that is, items that correlate negatively with the other items and therefore should be reversed. The function will then reverse those items prior to running Cronbach's alpha and provide a warning in red text that lets you know items were reversed:

Savoring_check<-alpha(Facebook[,25:48], check.keys=TRUE)

Another useful feature is delete, which, if set to true, will drop items with no variance (i.e., everyone responded to an item in the exact same way). This is less likely to occur in 1) established measures, where meaningless items have already been weeded out, and 2) measures using rating scales. For instance, when I assisted with psychometric analysis of a new measure of adverse outcomes to opioid prescriptions, we weeded out a few items that had no variance - but then, this was a brand new measure that was essentially a checklist (used a 0, not present and 1, present scale). If one or both of these apply to the measure you're working on, you might want to consider adding delete:

Savoring_drop<-alpha(Facebook[,25:48], delete=TRUE)

As with check.keys, you'll receive a warning.

You can read the document of the alpha function in psych here. There are many other options you can specify, such as number of iterations if you want to bootstrap your confidence intervals, or number of observations if you have a correlation matrix instead of raw data. What I've provided you above are the ones that you're most likely to need when examining alpha in a dataset.

Statistics Sunday: Cronbach's alpha

When developing a measure, there are two constructs that are very important: reliability and validity.

Validity, in this context, means the measure is measuring the right thing and not something else.

Reliability means the measure is consistent in measuring whatever it's measuring.

Obviously, having one does not automatically guarantee you'll have the other. Some measures that assess transient states can be highly valid, but because the thing they measure isn't consistent, the measures will appear to have lower reliability. And using shoe size to measure cognitive ability will have high reliability, even over time, but very low validity. So reliability refers to the consistency of the values while validity refers to the relationship between the values and the construct of interest.

There are different ways you can measure these two constructs. Today, I'll be focusing on reliability, and a specific measure of reliability: Cronbach's alpha. Look for my Blogging A to Z post today for how to compute Cronbach's alpha in R.

The type of reliability you want to measure affects the data you would collect. One type is test-retest reliability; if a person takes a measure on two separate occasions, how well do their scores match up? You would measure this by correlating scores from the first administration with scores from the second time. High test-retest reliability means scores are consistent across time. But, as I mention above, there might be situations where test-retest reliability doesn't make sense - that is, it may not be the type of reliability you're aiming for.

Other types of reliability can be measured with just one administration of a measure; these types of reliability are referred to as internal consistency. The simplest measure of internal consistency is split-half reliability. I literally divide my items in half and correlate scores on one half with scores on the other. I can split my items in a few ways - I could cut off at the halfway point, I could assign even-numbered items to one half and odd-numbered items to the other, or I could split items into halves at random.

Here's what those different kind of splits - halfway, even-odd, and random - might look like for a 6-item measure:

The way I split up my items will affect the overall correlation, though they should all be similar. But there still could be a fluke, where items that correlate poorly with the rest of the items just happen to be placed on the same half. One way you could correct for that is by running every possible split-half combination, then averaging those correlation results together.

That resulting average correlation of all possible split-halves is Cronbach's alpha.

When computing Cronbach's alpha, it's important to make sure items have the same meaning and there are no negative correlations. So if an item measures the reverse of a concept, you'd want to reverse score that item. For instance, say I created a measure of extraversion, where people rate their agreement with statements from 5, Strongly Agree to 1, Strongly Disagree. Since this is a measure of extraversion, we probably want higher scores to mean more extraversion.¹ But say I have one item worded as:

• I prefer to spend my free time alone. 

That item measures the opposite of extraversion. A person who strong disagrees with that item prefers not to spend free time alone, so I'd want to reverse-score that item, with Strongly Disagree being worth 5 points and Strongly Agree worth 1 point. Most software programs can do this kind of recode easily for you, so you shouldn't do it by hand. In fact, the R package I use for this measure will reverse items automatically when computing alpha, so you don't even need to create a new variable.

Because Cronbach's alpha is simply a correlation, and because you don't want any negative correlations (any items that negatively correlate with others should be reversed), Cronbach's alpha ranges from 0 to +1. Closer to 1 is better. There's some disagreement in the literature on how how high Cronbach's alpha needs to be. I usually use 0.8 as a cutoff - Cronbach's alpha below 0.8 suggests poor reliability - with 0.9 being optimal. Essentially, I consider:

• ≥ 0.9 - excellent reliability
• ≥ 0.8 but < 0.9 - acceptable reliability
• < 0.8 - poor reliability

I've seen some people consider values as low as 0.7 acceptable, and for certain measures, that could possibly be the case; as I said, there is some disagreement in the literature on this issue. I wouldn't go any lower than 0.7, though. And if you must sacrifice on reliability, make certain you have really strong evidence for the validity of your measure.

When should you use Cronbach's alpha? 

Cronbach's alpha is a classical test theory approach to reliability, so it makes sense to use it when creating a measure using classical test theory. Item response theory and Rasch measures use a different kind of reliability measure. I have seen Cronbach's alpha used with measures developed using item response theory or Rasch - mainly when the measure uses a fixed form (all examinees receive the exact same items). When creating a measure, I want reliability to be as high as possible, and I'm hesitant to accept anything below 0.9. If my reliability is less than 0.9, I go back to my items and see which ones are poor performers and should potentially be dropped from the final measure.

For instance, I was recently handed a large dataset collected for my company by a marketing research firm. There were groups of items that were believed to assess the same thing, but they weren't developed as a measure (i.e., with psychometric methods and analysis) and the research group did their initial analysis using individual items. (Can you say "p-hacking"?) I did some principal components analyses to make sure all items measure the same thing or to see if there were subscales - this is a data-driven technique, but I had fewer cases than I would have liked for a confirmatory factor analysis. Instead, I adopted a hybrid theory-driven and data-driven approach: I examined items ahead of time and grouped them together as subscales, then confirmed that the PCA found similar results (they did). Then I examined Cronbach's alpha for the subscales that emerged, as further support that they could be grouped together.

You should also compute Cronbach's alpha when using a measure created with classical test theory, and you'll want to make sure the reliability of the measure in your sample is high and comparable to the established reliability from measurement development research. (This is why it's important to track down psychometric articles/reports on the measures you're using. But being a psychometrician, of course I would say that.) In research reports, I've frequently been asked by reviewers to include Cronbach's alpha for all measures I used in a study. When I'm simply using someone else's measure, I'm less worried about super-high reliability; as long as it's close to 0.8, I'm fine with it. Depending on the measure, I've been fine with an alpha of 0.78, for instance. If it's much lower than that, I am usually hesitant to include that measure in my analysis.

Check back later today for the post on conducting Cronbach's alpha in R! And be sure to stop by Deeply Trivial again - there will be statistics-related posts every day this month!

¹One of the assumptions of Rasch and item response theory, the models I use in my psychometrics work, is that higher scores correspond to more of the construct being measured. This is how I approach measurement. Some researchers instead create measures where lower scores indicate more of the construct. There's nothing wrong with that approach if you're using classical test theory, but it feels backwards to me. Even when I use classical test theory, I use the higher scores equals more approach. You may prefer the opposite. The point is that all items need to be expressed in the same way - you don't want some items where higher scores equals more and others where lower scores equals more. So whichever method you adopt, make sure it's consistent and items that don't follow that method are reverse-scored.