Monday, April 30, 2018

Z is for Z-Scores and Standardizing

Z is for Z-Scores and Standardizing Last April, I wrapped up the A to Z of Statistics with a post about Z-scores. It seems only fitting that I'm wrapping this April A to Z with the same topic. Z-scores are frequently used, sometimes when you don't even realize it. When you take your child to the doctor and they say he's at the x percentile on height, or when you take a standardized test and are told you scored in the y percentile, those values are being derived from Z-scores. You create a Z-score when you subtract the population mean from a value and divide that result by the population standard deviation.

Of course, we often will standardize variables in statistics, and the results are similar to Z-scores (though technically not the same if the mean and standard deviation aren't population values). In fact, when I demonstrated the GLM function earlier this month, I skipped a very important step when conducting an analysis with interactions. I should have standardized my continuous predictors first, which means subtracting the variable mean and dividing by the variable standard deviation, creating a new variable with a mean of 0 and a standard deviation of 1 (just like the normal distribution).

Let's revisit that GLM analysis. I was predicting verdict (guilty, not guilty) with binomial regression. I did one analysis where I used a handful of attitude items and the participant's guilt rating, and a second analysis where I created interactions between each attitude item and the guilt rating. The purpose was to see if an attitude impacts the threshold - how high the guilt rating needed to be before a participant selected "guilty". When working with interactions, the individual variables are highly correlated with the interaction variables based on them, so we solve that problem, and make our analysis and output a bit cleaner, by centering our variables and using those centered values to create interactions.

I'll go ahead and load my data. Also, since I know I have some missing values, which caused an error when I tried to work with predicted values and residuals Saturday, I'll also go ahead and identify that case/those cases.

dissertation<-read.delim("dissertation_data.txt",header=TRUE)
predictors<-c("obguilt","reasdoubt","bettertolet","libertyvorder",
              "jurevidence","guilt")
library(psych)
## Warning: package 'psych' was built under R version 3.4.4
describe(dissertation[predictors])
##               vars   n mean   sd median trimmed  mad min max range  skew
## obguilt          1 356 3.50 0.89      4    3.52 0.00   1   5     4 -0.50
## reasdoubt        2 356 2.59 1.51      2    2.68 1.48  -9   5    14 -3.63
## bettertolet      3 356 2.36 1.50      2    2.38 1.48  -9   5    14 -3.41
## libertyvorder    4 355 2.74 1.31      3    2.77 1.48  -9   5    14 -3.89
## jurevidence      5 356 2.54 1.63      2    2.66 1.48  -9   5    14 -3.76
## guilt            6 356 4.80 1.21      5    4.90 1.48   2   7     5 -0.59
##               kurtosis   se
## obguilt          -0.55 0.05
## reasdoubt        26.92 0.08
## bettertolet      25.47 0.08
## libertyvorder    34.58 0.07
## jurevidence      25.39 0.09
## guilt            -0.54 0.06
dissertation<-subset(dissertation, !is.na(libertyvorder))

R has a built-in function that will do this for you: scale. The scale function has 3 main arguments - the variable or variables to be scaled, and whether you want those variables centered (subtract mean) and/or scaled (divided by standard deviation). For regression with interactions, we want to both center and scale. For instance, to center and scale the guilt rating:

dissertation$guilt_c<-scale(dissertation$guilt, center=TRUE, scale=TRUE)

I have a set of predictors I want to do this to, so I want to apply a function across those specific columns:

dissertation[46:51]<-lapply(dissertation[predictors], function(x) {
  y<-scale(x, center=TRUE, scale=TRUE)
  }
)

Now, let's rerun that binomial regression, this time using the centered variables in the model.

pred_int<-'verdict ~ obguilt.1 + reasdoubt.1 + bettertolet.1 + libertyvorder.1 + 
                  jurevidence.1 + guilt.1 + obguilt.1*guilt.1 + reasdoubt.1*guilt.1 +
                  bettertolet.1*guilt.1 + libertyvorder.1*guilt.1 + jurevidence.1*guilt.1'
model2<-glm(pred_int, family="binomial", data=dissertation)
summary(model2)
## 
## Call:
## glm(formula = pred_int, family = "binomial", data = dissertation)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6101  -0.5432  -0.1289   0.6422   2.2805  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.47994    0.16264  -2.951  0.00317 ** 
## obguilt.1                0.25161    0.16158   1.557  0.11942    
## reasdoubt.1             -0.09230    0.20037  -0.461  0.64507    
## bettertolet.1           -0.22484    0.20340  -1.105  0.26899    
## libertyvorder.1          0.05825    0.21517   0.271  0.78660    
## jurevidence.1            0.07252    0.19376   0.374  0.70819    
## guilt.1                  2.31003    0.26867   8.598  < 2e-16 ***
## obguilt.1:guilt.1        0.14058    0.23411   0.600  0.54818    
## reasdoubt.1:guilt.1     -0.61724    0.29693  -2.079  0.03764 *  
## bettertolet.1:guilt.1    0.02579    0.30123   0.086  0.93178    
## libertyvorder.1:guilt.1 -0.27492    0.29355  -0.937  0.34899    
## jurevidence.1:guilt.1    0.27601    0.36181   0.763  0.44555    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 490.08  on 354  degrees of freedom
## Residual deviance: 300.66  on 343  degrees of freedom
## AIC: 324.66
## 
## Number of Fisher Scoring iterations: 6

The results are essentially the same; the constant and slopes of the predictors variables are different but the variables that were significant before still are. So standardizing doesn't change the results, but it is generally recommended because it makes results easier to interpret. The variables are centered around the mean, so negative numbers are below the mean, and positive numbers are above the mean.

Hard to believe A to Z is over! Don't worry, I'm going to keep blogging about statistics, R, and whatever strikes my fancy. I almost kept this post going by applying the prediction work from yesterday to the binomial model, but decided that would make for a fun future post. And I'll probably sprinkle in posts in the near future on topics I didn't have room for this month or that I promised to write a future post on. Thanks for reading and hope you keep stopping by, even though April A to Z is officially over!

Sunday, April 29, 2018

Statistics Sunday: Conducting Meta-Analysis in R

Here it is, everyone! The promised 4th post on meta-analysis, and my second video for Deeply Trivial! In this video, I walk through conducting a basic meta-analysis, both fixed and random effects, in the metafor package:



See these previous posts and links for more information:
You can access the code I used in the video here as well as code to do similar analysis with the or_meta dataset here.

Saturday, April 28, 2018

Y is for Ys, Y-hats, and Residuals

Y is for Ys, Y-hats, and Residuals When working with a prediction model, like a linear regression, there are a few Ys you need to concern yourself with: the ys (observed outcome variable), the y-hats (predicted outcome variables based on the equation), and the residuals (y minus y-hat). Today, I'll dig into the different flavors of y and how you might work with them when conducting linear regression.

As my data example, I'll use my dissertation dataset and the linear model I introduced in the GLM post.

dissertation<-read.delim("dissertation_data.txt",header=TRUE)
guilt_lm_full<-lm(guilt ~ illev + circumst + deftest + policepower +
                    suspicious + overzealous + upstanding,
                  data=dissertation)
options(scipen = 999)
summary(guilt_lm_full)
## 
## Call:
## lm(formula = guilt ~ illev + circumst + deftest + policepower + 
##     suspicious + overzealous + upstanding, data = dissertation)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0357 -0.7452  0.1828  0.9706  2.5013 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  4.16081    0.38966  10.678 < 0.0000000000000002 ***
## illev        0.11111    0.05816   1.911              0.05689 .  
## circumst    -0.08779    0.06708  -1.309              0.19147    
## deftest     -0.02020    0.05834  -0.346              0.72942    
## policepower  0.02828    0.06058   0.467              0.64090    
## suspicious   0.17286    0.06072   2.847              0.00468 ** 
## overzealous -0.03298    0.04792  -0.688              0.49176    
## upstanding   0.08941    0.05374   1.664              0.09706 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.169 on 347 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.07647, Adjusted R-squared:  0.05784 
## F-statistic: 4.105 on 7 and 347 DF,  p-value: 0.0002387

In this model, the outcome variable is a guilt rating, ranging from 1 to 7. This is our y, which our regression model is trying to recreate through the linear relationship between our x's and our y. Using the coefficients in the output, we could compute the predicted y (y-hat) - what a person's score would be if the linear model fit the data perfectly. Fortunately, R has a built-in function that will compute y-hat for a dataset: predict. This function requires two arguments: a regression model and the dataset to use to predict values. Let's have R predict values for the dissertation dataset, and add it on as a new variable.

dissertation$predicted<-predict(guilt_lm_full, dissertation)

In this application, we don't care as much about the predicted values - we will later on in this post - but we probably do care about the residuals: the difference between the observed value and the predicted value. This gives us an idea of how our model is doing and whether it fits reasonably well. It can also tell us if the model falls apart at certain values or ranges of values.

In the residuals post, I showed that you can easily request residuals from the model. As we did with predicted, let's create a new variable in the dataset that contains our residuals.

dissertation$residual<-resid(guilt_lm_full)
## Error in `$<-.data.frame`(`*tmp*`, residual, value = structure(c(0.0326393185592984, : replacement has 355 rows, data has 356

Ruh-roh, we got an error. Our dataset contains 356 observations, but we only have 355 residuals. This is because someone has a missing value on one of the variables in the regression model and was dropped from the analysis. There are a variety of ways we could find out which case is missing a value, but since I'm only working with a handful of variables, I'll just run descriptives and look for the variable with only 355 values.

library(psych)
## Warning: package 'psych' was built under R version 3.4.4
describe(dissertation[c(13,15,18,21,27,29,31,44)])
##             vars   n mean   sd median trimmed  mad min max range  skew
## illev          1 356 2.98 1.13      3    3.02 1.48   1   5     4 -0.17
## circumst       2 356 2.99 0.95      3    2.97 1.48   1   5     4  0.13
## deftest        3 356 3.39 1.46      4    3.57 0.00  -9   5    14 -5.25
## policepower    4 355 3.86 1.41      4    4.02 0.00  -9   5    14 -6.40
## suspicious     5 356 2.09 1.14      2    2.01 1.48  -9   5    14 -1.97
## overzealous    6 356 3.34 1.34      4    3.41 1.48  -9   5    14 -4.49
## upstanding     7 356 3.09 1.29      3    3.11 1.48  -9   5    14 -2.31
## guilt          8 356 4.80 1.21      5    4.90 1.48   2   7     5 -0.59
##             kurtosis   se
## illev          -1.04 0.06
## circumst       -0.51 0.05
## deftest        40.74 0.08
## policepower    55.05 0.08
## suspicious     23.52 0.06
## overzealous    38.44 0.07
## upstanding     19.66 0.07
## guilt          -0.54 0.06

The variable policepower is the culprit. I can drop that missing value then rerun the residual code.

dissertation<-subset(dissertation, !is.na(policepower))
dissertation$residual<-resid(guilt_lm_full)

Now I can plot my observed values and residuals.

library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
qplot(guilt,residual,data=dissertation)

We want our residuals to fall around 0, which is only happening for guilt ratings near the midpoint of the scale. This suggests that ordinary least squares regression may not be the best fit for the data, as the model shows a tendency to overpredict (negative residual) guilt rating for people with lower observed ratings and underpredict (positive residual) for people with higher observed ratings.

But, as I often do on this blog for the sake of demonstration, let's pretend the model is doing well. One way we could use a regression model is to predict scores in a new sample. For instance, there are multiple rumors that different graduate schools have prediction equations they use to predict a candidate's anticipated graduate school GPA, based on a combination of factors asked about in the application packet, to determine if a person is grad school-ready (and ultimately, to decide if they should be admitted). Schools generally won't confirm they do this, nor would they ever share the prediction equation, should such an equation exist. But this is one way regression results from a training sample could be used to make decisions on a testing sample. So let's do that.

Unfortunately, I don't have a second dissertation dataset laying around that I could apply this equation to, but I could take a note from the data science playbook, and randomly divide my sample into training and testing datasets. I use the training dataset to generate my equation, and I use the testing dataset to apply my equation and predict values. Since I have outcome variable data in the testing dataset too, I can see how well my model did. Once I have a well-performing model, I could then apply it to new data, maybe to predict how highly you'll rate a book or movie, or to generate recommendations, or even to determine if I should let you in to the super-elite Monstersori school I want to create.

First, I'll split my dataset in half.

smp_size <- floor(0.50 * nrow(dissertation))

set.seed(42)
train_ind <- sample(seq_len(nrow(dissertation)), size = smp_size)

train <- dissertation[train_ind, ]
test <- dissertation[-train_ind, ]

Now I have a train dataset, with 177 observations, and a test dataset with 178. I can rerun the same linear model as before, this time only using the training data.

guilt_lm_train<-lm(guilt ~ illev + circumst + deftest + policepower +
               suspicious + overzealous + upstanding,
             data=train)
summary(guilt_lm_train)
## 
## Call:
## lm(formula = guilt ~ illev + circumst + deftest + policepower + 
##     suspicious + overzealous + upstanding, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9420 -0.8359  0.1641  0.9371  2.3151 
## 
## Coefficients:
##             Estimate Std. Error t value       Pr(>|t|)    
## (Intercept)  5.28874    0.77150   6.855 0.000000000128 ***
## illev        0.08866    0.08485   1.045        0.29759    
## circumst    -0.13018    0.09917  -1.313        0.19109    
## deftest     -0.25726    0.10699  -2.405        0.01727 *  
## policepower  0.01758    0.12316   0.143        0.88665    
## suspicious   0.25716    0.08857   2.903        0.00419 ** 
## overzealous -0.11683    0.08240  -1.418        0.15807    
## upstanding   0.10371    0.07574   1.369        0.17273    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.194 on 169 degrees of freedom
## Multiple R-squared:  0.1265, Adjusted R-squared:  0.09027 
## F-statistic: 3.495 on 7 and 169 DF,  p-value: 0.001586

I can use my predict function to predict scores for the testing dataset. Remember, all this function needs is the linear model name and a dataset to use for the prediction function - and it can be any dataset, as long as it contains the same variables from the model.

test$predicted2<-predict(guilt_lm_train, test)

The original predicted value (from when I was working with the full dataset) is still in this set. I could have replaced values by using the same variable name, but just for fun, decided to keep those values and create a second prediction variable.

Because we have observed and predicted2 for our training dataset, let's see how well our model did, by creating a new residual variable, residual2. We can't use the resid function, because we didn't have R perform a linear regression on the testing dataset, but we can easily generate this variable by subtracting the predicted score from the observed score. Then we can once again plot our observed and residual values.

test$residual2<-test$guilt-test$predicted2
qplot(guilt,residual2,data=test)

We're still seeing similar issues with the residuals as we did for the full dataset. If we wanted to actually apply our linear model, we'd want to do more research and pilot work to get the best equation we can. As with many things in statistics, the process is heavily determined by what you plan to do with the results. If you want to report variables that have a strong linear relationship with an outcome, we might be fine with these regression results. If we want to build an equation that predicts an outcome with a minimum of error, we'd want to do more exploratory work, pulling in new variables and dropping low-performing ones. Many of the predictors in the model were not significant, so we could start model-building from those results. We may need to build multiple training datasets, to ensure we aren't only picking up chance relationships. And for much larger applications, such as recommendation systems on services like Amazon and Netflix, machine learning may be a better, more powerful method.

One more A to Z post left!

Friday, April 27, 2018

WRiTE CLUB 2018

This month is an awesome writing contest called WRiTE CLUB, which faces off 500-word excerpts each day over 15 days. (This is round 1 of the event. There will be new bouts and authors who advance will also have the chance to share additional excerpts.) Ten bouts have been posted so far with 5 more next week. Readers are invited to vote each day, with the chance to win prizes just for voting. And the winning author also wins prizes.

Voting on each bout is open for a week, so voting is still open for many of the bouts. You can access them here.

X is for By

X is for By Today's post will be rather short, demonstrating a set of functions from the psych package, which allows you to conduct analysis by group. These commands add "By" to the end of existing functions. But first, a word of caution: With great power comes great responsibility. This function could very easily turn into a fishing expedition (also known as p-hacking). Conducting planned group comparisons is fine. Conducting all possible group comparisons and cherry-picking any differences is problematic. So use these group by functions with care.

Let's pull up the Facebook dataset for this.

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

This is the full dataset, which includes all the variables I collected. I don't want to run analyses on all variables, so I'll pull out the ones most important for this blog post demonstration.

smallFB<-Facebook[,c(1:2,77:80,105:116,122,133:137,170,187)]

First, I'll run descriptives on this smaller data frame by gender.

library(psych)
## Warning: package 'psych' was built under R version 3.4.4
describeBy(smallFB,smallFB$gender)
## 
##  Descriptive statistics by group 
## group: 0
##              vars  n      mean      sd   median   trimmed     mad      min
## RespondentId    1 73 164647.77 1711.78 164943.0 164587.37 2644.96 162373.0
## gender          2 73      0.00    0.00      0.0      0.00    0.00      0.0
## Rumination      3 73     37.66   14.27     37.0     37.41   13.34      8.0
## DepRelat        4 73     21.00    7.86     21.0     20.95    5.93      4.0
## Brood           5 73      8.49    3.76      9.0      8.42    2.97      1.0
## Reflect         6 73      8.16    4.44      8.0      8.24    4.45      0.0
## SavorPos        7 73     64.30   10.93     65.0     64.92    8.90     27.0
## SavorNeg        8 73     33.30   11.48     33.0     33.08   13.34     12.0
## SavorTot        9 73     31.00   20.15     34.0     31.15   19.27    -10.0
## AntPos         10 73     20.85    3.95     21.0     20.93    4.45     10.0
## AntNeg         11 73     11.30    4.23     11.0     11.22    4.45      4.0
## AntTot         12 73      9.55    6.90     10.0      9.31    7.41     -3.0
## MomPos         13 73     21.68    3.95     22.0     21.90    2.97      9.0
## MomNeg         14 73     11.45    4.63     11.0     11.41    5.93      4.0
## MomTot         15 73     10.23    7.63     11.0     10.36    8.90    -11.0
## RemPos         16 73     21.77    4.53     23.0     22.20    4.45      8.0
## RemNeg         17 73     10.55    4.39      9.0     10.27    4.45      4.0
## RemTot         18 73     11.22    8.05     14.0     11.68    7.41     -8.0
## LifeSat        19 73     24.63    6.80     25.0     24.93    7.41     10.0
## Extravert      20 73      4.32    1.58      4.5      4.33    1.48      1.5
## Agreeable      21 73      4.79    1.08      5.0      4.85    1.48      1.0
## Conscient      22 73      5.14    1.34      5.0      5.19    1.48      2.0
## EmotStab       23 73      5.10    1.22      5.0      5.15    1.48      1.0
## OpenExp        24 73      5.11    1.29      5.5      5.20    1.48      2.0
## Health         25 73     28.77   19.56     25.0     26.42   17.79      0.0
## Depression     26 73     10.26    7.27      9.0      9.56    5.93      0.0
##                 max  range  skew kurtosis     se
## RespondentId 168279 5906.0  0.21    -1.36 200.35
## gender            0    0.0   NaN      NaN   0.00
## Rumination       71   63.0  0.12    -0.53   1.67
## DepRelat         42   38.0  0.10    -0.04   0.92
## Brood            17   16.0  0.15    -0.38   0.44
## Reflect          19   19.0 -0.12    -0.69   0.52
## SavorPos         84   57.0 -0.69     0.76   1.28
## SavorNeg         57   45.0  0.14    -0.95   1.34
## SavorTot         72   82.0 -0.17    -0.75   2.36
## AntPos           28   18.0 -0.24    -0.46   0.46
## AntNeg           22   18.0  0.27    -0.55   0.49
## AntTot           24   27.0  0.11    -0.76   0.81
## MomPos           28   19.0 -0.69     0.55   0.46
## MomNeg           22   18.0  0.08    -0.98   0.54
## MomTot           24   35.0 -0.25    -0.55   0.89
## RemPos           28   20.0 -0.88     0.35   0.53
## RemNeg           22   18.0  0.56    -0.66   0.51
## RemTot           24   32.0 -0.53    -0.77   0.94
## LifeSat          35   25.0 -0.37    -0.84   0.80
## Extravert         7    5.5 -0.09    -0.93   0.19
## Agreeable         7    6.0 -0.60     1.04   0.13
## Conscient         7    5.0 -0.24    -0.98   0.16
## EmotStab          7    6.0 -0.60     0.28   0.14
## OpenExp           7    5.0 -0.49    -0.55   0.15
## Health           91   91.0  1.13     1.14   2.29
## Depression       36   36.0  1.02     0.95   0.85
## -------------------------------------------------------- 
## group: 1
##              vars   n      mean      sd    median   trimmed     mad
## RespondentId    1 184 164373.49 1515.34 164388.00 164253.72 1891.80
## gender          2 184      1.00    0.00      1.00      1.00    0.00
## Rumination      3 184     38.09   15.28     40.00     38.16   17.05
## DepRelat        4 184     21.67    8.78     21.00     21.66    8.90
## Brood           5 184      8.57    4.14      8.50      8.47    3.71
## Reflect         6 184      7.84    4.06      8.00      7.73    4.45
## SavorPos        7 184     67.22    9.63     68.00     67.71    8.90
## SavorNeg        8 184     29.75   11.62     27.50     28.72    9.64
## SavorTot        9 184     37.47   19.30     40.00     38.66   20.02
## AntPos         10 184     22.18    3.37     23.00     22.28    2.97
## AntNeg         11 184     10.10    4.44      9.00      9.78    4.45
## AntTot         12 184     12.08    6.85     14.00     12.36    5.93
## MomPos         13 184     22.28    3.88     23.00     22.59    2.97
## MomNeg         14 184     10.60    4.88      9.50     10.13    5.19
## MomTot         15 184     11.68    7.75     13.00     12.29    7.41
## RemPos         16 184     22.76    3.85     23.00     23.10    2.97
## RemNeg         17 184      9.05    3.79      8.00      8.68    2.97
## RemTot         18 184     13.71    6.97     15.00     14.34    5.93
## LifeSat        19 184     23.76    6.25     24.00     24.18    7.41
## Extravert      20 184      4.66    1.57      5.00      4.74    1.48
## Agreeable      21 184      5.22    1.06      5.50      5.26    1.48
## Conscient      22 184      5.32    1.24      5.50      5.42    1.48
## EmotStab       23 184      4.70    1.31      4.75      4.75    1.11
## OpenExp        24 184      5.47    1.08      5.50      5.56    0.74
## Health         25 184     32.54   16.17     30.00     31.43   16.31
## Depression     26 184     12.19    8.48      9.00     11.09    5.93
##                   min    max  range  skew kurtosis     se
## RespondentId 162350.0 167714 5364.0  0.46    -0.90 111.71
## gender            1.0      1    0.0   NaN      NaN   0.00
## Rumination        3.0     74   71.0 -0.05    -0.60   1.13
## DepRelat          0.0     42   42.0  0.00    -0.46   0.65
## Brood             0.0     19   19.0  0.19    -0.62   0.31
## Reflect           0.0     19   19.0  0.25    -0.48   0.30
## SavorPos         33.0     84   51.0 -0.59     0.36   0.71
## SavorNeg         12.0     64   52.0  0.79     0.25   0.86
## SavorTot        -18.0     72   90.0 -0.57    -0.10   1.42
## AntPos            9.0     28   19.0 -0.49     0.41   0.25
## AntNeg            4.0     22   18.0  0.63    -0.39   0.33
## AntTot           -8.0     24   32.0 -0.43    -0.48   0.50
## MomPos           10.0     28   18.0 -0.81     0.54   0.29
## MomNeg            4.0     24   20.0  0.81    -0.03   0.36
## MomTot          -13.0     24   37.0 -0.69    -0.03   0.57
## RemPos            9.0     28   19.0 -0.87     0.81   0.28
## RemNeg            4.0     21   17.0  0.83     0.33   0.28
## RemTot           -9.0     24   33.0 -0.82     0.50   0.51
## LifeSat           8.0     35   27.0 -0.53    -0.32   0.46
## Extravert         1.0      7    6.0 -0.36    -0.72   0.12
## Agreeable         2.5      7    4.5 -0.27    -0.63   0.08
## Conscient         1.0      7    6.0 -0.70     0.13   0.09
## EmotStab          1.5      7    5.5 -0.35    -0.73   0.10
## OpenExp           1.5      7    5.5 -0.91     0.62   0.08
## Health            2.0     85   83.0  0.60    -0.05   1.19
## Depression        0.0     39   39.0  1.14     0.66   0.62

In this dataset, I coded men as 0 and women as 1. The descriptive statistics table generated includes all scale and subscale scores, and gives me mean, standard deviation, median, a trimmed mean (dropping very low and very high values), median absolute deviation, minimum and maximum values, range, skewness, and kurtosis. I'd need to run t-tests to find out if differences were significant, but this still gives me some idea of how men and women might differ on these measures.

There are certain measures I included that we might hypothesize would show gender differences. For instance, some research suggests gender differences for rumination and depression. In addition to running descriptives by group, I might also want to display these differences in a violin plot. The psych package can quickly generate such a plot by group.

violinBy(smallFB,"Rumination","gender",grp.name=c("M","F"))
violinBy(smallFB,"Depression","gender",grp.name=c("M","F"))

ggplot2 will generate a violin plot by group, so this feature might not be as useful for final displays, but could help in quickly visualizing the data during analysis. And you may find that you prefer the appearance of this plots. To each his own.

Another function is error.bars.by, which plots means and confidence intervals by group for multiple variables. Again, this is a way to get some quick visuals, though differences in scale among measures should be taken into consideration when generating this plot. One set of variables for which this display might be useful is the 5 subscales of the Five-Factor Personality Inventory. This 10-item measure assesses where participants fall on the so-called Big Five personality traits: Openness to Experience, Conscientiousness, Extraversion, Agreeableness, and Neuroticism (Emotional Stability). These subscales are all on the same metric.

error.bars.by(smallFB[,c(20:24)],group=smallFB$gender,xlab="Big Five Personality Traits",ylab="Score on Subscale")

Finally, we have the statsBy function, which gives descriptive statistics by group as well as between group statistics. This functions generates a lot of output, and you can read more about everything it gives you here.
FBstats<-statsBy(smallFB[,2:26],"gender",cors=TRUE,method="pearson",use="pairwise")
print(FBstats,short=FALSE)
## Statistics within and between groups  
## Call: statsBy(data = smallFB[, 2:26], group = "gender", cors = TRUE, 
##     method = "pearson", use = "pairwise")
## Intraclass Correlation 1 (Percentage of variance due to groups) 
##     gender Rumination   DepRelat      Brood    Reflect   SavorPos 
##       1.00      -0.01      -0.01      -0.01      -0.01       0.03 
##   SavorNeg   SavorTot     AntPos     AntNeg     AntTot     MomPos 
##       0.03       0.04       0.05       0.02       0.05       0.00 
##     MomNeg     MomTot     RemPos     RemNeg     RemTot    LifeSat 
##       0.00       0.01       0.02       0.05       0.04       0.00 
##  Extravert  Agreeable  Conscient   EmotStab    OpenExp     Health 
##       0.01       0.05       0.00       0.03       0.03       0.01 
## Depression 
##       0.01 
## Intraclass Correlation 2 (Reliability of group differences) 
##     gender Rumination   DepRelat      Brood    Reflect   SavorPos 
##       1.00     -22.34      -2.06     -50.93      -2.21       0.77 
##   SavorNeg   SavorTot     AntPos     AntNeg     AntTot     MomPos 
##       0.80       0.83       0.86       0.75       0.86       0.19 
##     MomNeg     MomTot     RemPos     RemNeg     RemTot    LifeSat 
##       0.39       0.46       0.68       0.87       0.84      -0.04 
##  Extravert  Agreeable  Conscient   EmotStab    OpenExp     Health 
##       0.60       0.88       0.05       0.80       0.81       0.60 
## Depression 
##       0.66 
## eta^2 between groups  
## Rumination.bg   DepRelat.bg      Brood.bg    Reflect.bg   SavorPos.bg 
##          0.00          0.00          0.00          0.00          0.02 
##   SavorNeg.bg   SavorTot.bg     AntPos.bg     AntNeg.bg     AntTot.bg 
##          0.02          0.02          0.03          0.02          0.03 
##     MomPos.bg     MomNeg.bg     MomTot.bg     RemPos.bg     RemNeg.bg 
##          0.00          0.01          0.01          0.01          0.03 
##     RemTot.bg    LifeSat.bg  Extravert.bg  Agreeable.bg  Conscient.bg 
##          0.02          0.00          0.01          0.03          0.00 
##   EmotStab.bg    OpenExp.bg     Health.bg Depression.bg 
##          0.02          0.02          0.01          0.01 
## Correlation between groups 
##               Rmnt. DpRl. Brd.b Rflc. SvrP. SvrN. SvrT. AntP. AntN. AntT.
## Rumination.bg  1                                                         
## DepRelat.bg    1     1                                                   
## Brood.bg       1     1     1                                             
## Reflect.bg    -1    -1    -1     1                                       
## SavorPos.bg    1     1     1    -1     1                                 
## SavorNeg.bg   -1    -1    -1     1    -1     1                           
## SavorTot.bg    1     1     1    -1     1    -1     1                     
## AntPos.bg      1     1     1    -1     1    -1     1     1               
## AntNeg.bg     -1    -1    -1     1    -1     1    -1    -1     1         
## AntTot.bg      1     1     1    -1     1    -1     1     1    -1     1   
## MomPos.bg      1     1     1    -1     1    -1     1     1    -1     1   
## MomNeg.bg     -1    -1    -1     1    -1     1    -1    -1     1    -1   
## MomTot.bg      1     1     1    -1     1    -1     1     1    -1     1   
## RemPos.bg      1     1     1    -1     1    -1     1     1    -1     1   
## RemNeg.bg     -1    -1    -1     1    -1     1    -1    -1     1    -1   
## RemTot.bg      1     1     1    -1     1    -1     1     1    -1     1   
## LifeSat.bg    -1    -1    -1     1    -1     1    -1    -1     1    -1   
## Extravert.bg   1     1     1    -1     1    -1     1     1    -1     1   
## Agreeable.bg   1     1     1    -1     1    -1     1     1    -1     1   
## Conscient.bg   1     1     1    -1     1    -1     1     1    -1     1   
## EmotStab.bg   -1    -1    -1     1    -1     1    -1    -1     1    -1   
## OpenExp.bg     1     1     1    -1     1    -1     1     1    -1     1   
## Health.bg      1     1     1    -1     1    -1     1     1    -1     1   
## Depression.bg  1     1     1    -1     1    -1     1     1    -1     1   
##               MmPs. MmNg. MmTt. RmPs. RmNg. RmTt. LfSt. Extr. Agrb. Cnsc.
## MomPos.bg      1                                                         
## MomNeg.bg     -1     1                                                   
## MomTot.bg      1    -1     1                                             
## RemPos.bg      1    -1     1     1                                       
## RemNeg.bg     -1     1    -1    -1     1                                 
## RemTot.bg      1    -1     1     1    -1     1                           
## LifeSat.bg    -1     1    -1    -1     1    -1     1                     
## Extravert.bg   1    -1     1     1    -1     1    -1     1               
## Agreeable.bg   1    -1     1     1    -1     1    -1     1     1         
## Conscient.bg   1    -1     1     1    -1     1    -1     1     1     1   
## EmotStab.bg   -1     1    -1    -1     1    -1     1    -1    -1    -1   
## OpenExp.bg     1    -1     1     1    -1     1    -1     1     1     1   
## Health.bg      1    -1     1     1    -1     1    -1     1     1     1   
## Depression.bg  1    -1     1     1    -1     1    -1     1     1     1   
##               EmtS. OpnE. Hlth. Dprs.
## EmotStab.bg    1                     
## OpenExp.bg    -1     1               
## Health.bg     -1     1     1         
## Depression.bg -1     1     1     1   
## Correlation within groups 
##               Rmnt. DpRl. Brd.w Rflc. SvrP. SvrN. SvrT. AntP. AntN. AntT.
## Rumination.wg  1.00                                                      
## DepRelat.wg    0.95  1.00                                                
## Brood.wg       0.88  0.78  1.00                                          
## Reflect.wg     0.80  0.63  0.59  1.00                                    
## SavorPos.wg   -0.20 -0.20 -0.18 -0.15  1.00                              
## SavorNeg.wg    0.43  0.43  0.36  0.30 -0.64  1.00                        
## SavorTot.wg   -0.36 -0.36 -0.31 -0.25  0.89 -0.92  1.00                  
## AntPos.wg     -0.06 -0.05 -0.08 -0.03  0.86 -0.49  0.73  1.00            
## AntNeg.wg      0.32  0.32  0.28  0.21 -0.54  0.89 -0.80 -0.50  1.00      
## AntTot.wg     -0.23 -0.23 -0.21 -0.15  0.78 -0.82  0.89  0.83 -0.89  1.00
## MomPos.wg     -0.26 -0.26 -0.22 -0.19  0.86 -0.60  0.80  0.60 -0.47  0.61
## MomNeg.wg      0.46  0.46  0.39  0.35 -0.51  0.88 -0.78 -0.33  0.66 -0.59
## MomTot.wg     -0.42 -0.42 -0.36 -0.32  0.75 -0.85  0.89  0.51 -0.65  0.68
## RemPos.wg     -0.20 -0.19 -0.17 -0.15  0.89 -0.56  0.79  0.66 -0.44  0.62
## RemNeg.wg      0.34  0.35  0.28  0.23 -0.65  0.87 -0.85 -0.49  0.69 -0.69
## RemTot.wg     -0.29 -0.30 -0.25 -0.21  0.85 -0.79  0.90  0.63 -0.62  0.72
## LifeSat.wg    -0.47 -0.47 -0.43 -0.31  0.54 -0.50  0.57  0.39 -0.33  0.41
## Extravert.wg  -0.20 -0.19 -0.11 -0.20  0.34 -0.35  0.38  0.21 -0.29  0.29
## Agreeable.wg  -0.18 -0.18 -0.20 -0.10  0.35 -0.45  0.45  0.28 -0.39  0.39
## Conscient.wg  -0.25 -0.30 -0.20 -0.10  0.24 -0.21  0.25  0.16 -0.14  0.17
## EmotStab.wg   -0.48 -0.44 -0.49 -0.34  0.34 -0.44  0.43  0.20 -0.33  0.32
## OpenExp.wg    -0.16 -0.14 -0.21 -0.10  0.37 -0.31  0.37  0.27 -0.27  0.31
## Health.wg      0.44  0.47  0.36  0.29 -0.30  0.34 -0.35 -0.21  0.26 -0.27
## Depression.wg  0.57  0.58  0.49  0.38 -0.44  0.55 -0.55 -0.27  0.39 -0.39
##               MmPs. MmNg. MmTt. RmPs. RmNg. RmTt. LfSt. Extr. Agrb. Cnsc.
## MomPos.wg      1.00                                                      
## MomNeg.wg     -0.56  1.00                                                
## MomTot.wg      0.86 -0.91  1.00                                          
## RemPos.wg      0.65 -0.42  0.59  1.00                                    
## RemNeg.wg     -0.55  0.63 -0.67 -0.65  1.00                              
## RemTot.wg      0.66 -0.58  0.69  0.91 -0.91  1.00                        
## LifeSat.wg     0.55 -0.55  0.62  0.48 -0.42  0.49  1.00                  
## Extravert.wg   0.39 -0.37  0.43  0.28 -0.25  0.29  0.27  1.00            
## Agreeable.wg   0.33 -0.43  0.43  0.31 -0.36  0.37  0.25  0.12  1.00      
## Conscient.wg   0.25 -0.16  0.22  0.23 -0.26  0.26  0.33  0.03  0.29  1.00
## EmotStab.wg    0.40 -0.50  0.51  0.27 -0.32  0.32  0.44  0.12  0.41  0.27
## OpenExp.wg     0.39 -0.26  0.36  0.30 -0.28  0.32  0.34  0.29  0.36  0.14
## Health.wg     -0.30  0.33 -0.36 -0.27  0.29 -0.31 -0.42 -0.10 -0.25 -0.24
## Depression.wg -0.45  0.56 -0.58 -0.41  0.49 -0.50 -0.65 -0.24 -0.29 -0.26
##               EmtS. OpnE. Hlth. Dprs.
## EmotStab.wg    1.00                  
## OpenExp.wg     0.24  1.00            
## Health.wg     -0.31 -0.18  1.00      
## Depression.wg -0.54 -0.28  0.56  1.00
## 
## Many results are not shown directly. To see specific objects select from the following list:
##  mean sd n F ICC1 ICC2 ci1 ci2 r within pooled sd.r raw rbg pbg rwg nw pwg etabg etawg nwg nG Call

The variance explained by gender is quite small for all of the variables. Instead, the relationships between the variables seem to be more meaningful.

A to Z is almost done! Just Y and Z, plus look for an A-to-Z-influenced Statistics Sunday post!

Thursday, April 26, 2018

Predictive Analytics and Veteran Suicide Prevention

One of the best things about working for the Department of Veterans Affairs was the vast amount of data available on Veterans receiving care through VA. While my research center often used surveys, focus groups, and interviews to collect data on Veterans, we frequently pulled in data from Veterans' medical records (with their permission, of course). And other researchers were accessing Veteran data directly to understand and improve care.

An issue we frequently heard about, and sometimes dealt with firsthand, was the high rate of suicide among our Veterans. Veterans are at high risk for many physical and mental conditions, and are at heightened risk for suicide. The National Suicide Prevention Lifeline was created to help anyone, including Veterans, who is feeling helpless. Through partnerships with VA and other federal agencies, we heard about many success stories of the Lifeline.

But with the vast amount of data available on our Veterans, it would be great if we could intervene and help before someone gets to the crisis point. Today, I read about how the REACH Vet program is using predictive analytics to identify Veterans at risk for suicide:
The REACH Vet program draws on the agency’s vast trove of electronic health records and uses predictive analytics to identify patients who might be at risk of suicide. It alerts VA clinicians of veterans who could benefit from more attention, and the program prompts clinicians to call and check in with their patients.

“What we found … not surprisingly is that veterans at highest risk of suicide are also at very high risk of some other things,” Aaron Eagan, VA’s deputy director for innovation, said Thursday at ACT-IAC’s Health Innovation Day in Washington. “They’re at significantly increased rates of all-cause mortality, accident morality, overdoses, violence … [and] opioids.”

Veterans who engaged with REACH Vet were admitted to mental health inpatient units less often, showed up to more mental health and primary care appointments and visited the VA more frequently, compared to veterans who weren’t part of the program.

Eagan said he expected veterans would be frustrated by the phone calls, but his team hasn’t gotten any complaints.

“It’s a great reminder that people really feel good about us caring about them, and that’s what the response generally is,” he said.

The REACH Vet team is updating its predictive model for the program now, and it’s starting a new collaboration with the Energy Department’s super computer, Eagan said.
 

W is for (Meta-Analysis) Weights

Weights in Meta-Analysis Yesterday, I talked about the variance calculation for meta-analysis, which is based on sample size - the exact calculation varies by the type of effect size you're using. It's a good idea to inspect those variance calculations. There are many places where your numbers for the meta-analysis can be incorrect.

One place is when extracting that information from the studies themselves, which is why having more than one person extract information and calculating some sort of interrater reliability is a good idea. Another place is at data entry. You should always double-check those numbers, either with a visual examination of the dataset or, my preference especially for large meta-analyses, through data visualization, as I did yesterday with ggplot. Values that don't show an inverse relationship between sample size and variance should be checked. This was the case with study 045 in the or_meta file. A visual check, and knowledge of the formula for variance, confirmed that this variance was correct; it was large because there were a lot more people in the control group than treatment group.

Once we're satisfied that our numbers are correct, we can move forward with the meta-analysis. Part of that calculation involves study weights. Weights and variance have an inverse relationship; in fact, one value can be used to compute the other:

wi = 1/vi
vi = 1/wi

When you take the inverse of a small number, such as variance, you get a large number, weight. And when you take the inverse of a large number, you get a small number. metafor calculates those weights for you once you conduct the actual meta-analysis, using the rma function. So let's do that. First, I'll load my data frames, and calculate effect sizes and variances.

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)
)

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)
)

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).
smd_meta <- escalc(measure="SMD", m1i=m1, m2i=m2, sd1i=sd1, sd2i=sd2, n1i=n1,
                   n2i=n2,data=smd_meta)
or_meta <- escalc(measure="OR", ai=tg, bi=(tn-tg), ci=cg, di=(cn-cg),
                  data=or_meta)

These two datasets use different effect size metrics, so I can't combine them unless I convert the effect sizes from one dataset. For now, we'll meta-analyze them separately. First up will be the dataset using standardized mean difference (Cohen's d).

There are two overall methods of meta-analysis - fixed effects and random effects. Fixed effects means that you suspect there is only one underlying effect size that your studies are trying to estimate. You can access this by setting method="FE" (for fixed effects) in the rma function.

smd.rma<-rma(yi,vi,method="FE",data=smd_meta)
summary(smd.rma)
## 
## Fixed-Effects Model (k = 11)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -41.4601   97.3274   84.9202   85.3181   85.3647  
## 
## Test for Heterogeneity: 
## Q(df = 10) = 97.3274, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub     
##   0.3492  0.0526  6.6388  <.0001  0.2461  0.4523  ***
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I'll go into the various parts of this output on Sunday. For now, we'll look at the Model Results, which includes the estimate (the overall effect size).

smd.rma$beta
##              [,1]
## intrcpt 0.3491923

The overall Cohen's d for this set of studies is 0.349, a moderate effect. In plain English, the average guilt rating among people who saw pretrial publicity is about 0.349 standard deviations higher than people in the control group. There's an associated Z-value and p-value, which tests whether that metric is significantly different from 0.

options(scipen=999)
smd.rma$zval
## [1] 6.638784
smd.rma$pval
## [1] 0.00000000003162821

The Z-test confirms that this value is significantly different from 0. I'll go into conducting a random effects model (and how to tell if you should) on Sunday. Now, let's run the meta-analysis for the log odds ratio dataset.

or.rma<-rma(yi,vi,method="FE",data=or_meta)
summary(or.rma)
## 
## Fixed-Effects Model (k = 13)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -23.0842   47.8107   48.1684   48.7334   48.5321  
## 
## Test for Heterogeneity: 
## Q(df = 12) = 47.8107, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub     
##   0.7456  0.0987  7.5566  <.0001  0.5522  0.9389  ***
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
or.rma$beta
##              [,1]
## intrcpt 0.7455659
or.rma$zval
## [1] 7.556647
or.rma$pval
## [1] 0.00000000000004135935

Our estimated log odds ratio for this dataset is 0.7456, which is significantly different from 0. This number is difficult to interpret in isolation, so we want to convert it back to odds ratio, which is easier to interpret.

exp(or.rma$beta)
##             [,1]
## intrcpt 2.107634

Overall, people who saw pretrial publicity were over 2 times as likely to convict as people in the control group.

metafor will also let you supply user-defined weights, using the weights argument in the rma function. If you're applying a correction to your effect sizes - such as converting your Pearson's r coefficients to Fisher's Z - you'll probably want to supply custom weights, because the formula is different. If you leave that argument out, as I did, rma will create weights as the inverse of variance by default.

I should note that I only picked a handful of studies to include for this example, so these don't match the results I found in my actual meta-analysis, which included many more studies. (If you're conducting research on pretrial publicity and would like the actual results of my meta-analysis, please contact me. You can find contact info on my about page. I'm happy to share my actual results. These numbers are for demonstration purposes only. I still hope to publish this work someday, which is why I'm not using the full dataset or actual results here, but at this point, I need to collect more data on the new pretrial publicity studies that have been conducted since 2011.)

More on meta-analysis Sunday! And tune in tomorrow for the letter X!

Wednesday, April 25, 2018

V is for (Meta-Analysis) Variance

Variance in Meta-Analysis For the letter E, I introduced the metafor package to compute effect sizes. That is, you provide a data frame with the study information and the data needed to compute the effect size(s), and metafor does that for you. But what I didn't mention is that metafor also gives you additional information that will become very important when conducting an actual meta-analysis - and that information will be the focus of today's and tomorrow's posts.

Today, I'll talk about the concept of variance in meta-analysis and discuss how metafor gives you that information. But first, we should set up our data frames with study information and data. I'll be reusing the two data frames I created for the E post - one to compute standardized mean difference (Cohen's d) and one to compute log odds ratio.

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)
)

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)
)

Now we'll rerun the code from that previous post that we used to generate effect sizes. We reference the original data frame when we run this code if we want the computed variables to be appended to the existing data frame.

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).
smd_meta <- escalc(measure="SMD", m1i=m1, m2i=m2, sd1i=sd1, sd2i=sd2, n1i=n1,
                   n2i=n2,data=smd_meta)
or_meta <- escalc(measure="OR", ai=tg, bi=(tn-tg), ci=cg, di=(cn-cg),
                  data=or_meta)

Let's take a quick look at the updated data frame, starting with smd_meta.

smd_meta
##     id study      author_year  n1  n2   m1   m2  sd1  sd2      yi     vi
## 1  005     1        Ruva 2007 138 138 5.29 4.08 1.65 1.67  0.7269 0.0154
## 2  005     2        Ruva 2007 140 142 5.05 3.89 1.50 1.61  0.7433 0.0152
## 3  029     3 Chrzanowski 2006 144 234 1.97 2.45 1.08 1.22 -0.4099 0.0114
## 4  031     1  Studebaker 2000  21  21 5.95 3.67 1.02 1.20  2.0087 0.1433
## 5  038     1        Ruva 2008  54  52 5.07 3.96 1.65 1.76  0.6464 0.0397
## 6  041     1    Bradshaw 2007  78  20 6.22 5.75 2.53 2.17  0.1893 0.0630
## 7  041     2    Bradshaw 2007  92  18 5.47 4.89 2.31 2.59  0.2444 0.0667
## 8  058     1      Wilson 1998  31  15 6.13 3.80 2.51 2.68  0.8927 0.1076
## 9  058     2      Wilson 1998  29  13 5.69 3.61 2.51 2.78  0.7867 0.1188
## 10 067     1   Locatelli 2011  90  29 4.81 4.61 1.20 1.39  0.1592 0.0457
## 11 067     2   Locatelli 2011 181  53 4.83 4.51 1.19 1.34  0.2603 0.0245

You'll notice that, in addition to the information from the data frame creation code above, I have two new variables: yi and vi. yi refers to the effect size for that study. vi is the variance associated with the effect size for that study. This information is important for computing the study weight, which I'll talk about tomorrow, as well as generating confidence intervals around the overall effect size, which I'll talk about Sunday. But you might be wondering where this variance comes from. For that, stand back - I'm about to use math.

Meta-analysis uses weights based on sample size. The reason for that is because meta-analysis is attempting to take information from multiple studies on a topic and estimate the underlying effect size that these studies are trying to estimate. We weight on sample size because studies done on more people are expected to be better able to estimate the true population value, so larger studies get more weight. On the flip side, larger studies are expected to have more precision in their estimate, so larger studies have smaller variance. In meta-analysis, the variance is estimated based on sample size. For standardized mean difference, the variance is computed as:

vi = ((n1+n2)/(n1*n2)) + (yi2/(2*(n1+n2)))

So for the first effect size, which is 0.7269, we can recreate the computed variance (within some rounding error):

((138+138)/(138*138)) + (0.7269/(2*(138+138)))
## [1] 0.0158096

That calculation is done for each effect size, and fortunately, metafor will do it for you. And since sample size is part of the calculation, you'll note that variances are larger for smaller studies and smaller for larger studies.

library(ggplot2)
ggplot(smd_meta, aes(x=vi, y=(n1+n2))) + geom_point() +
  geom_text(aes(label=id),hjust=0, vjust=0) +
  scale_x_continuous(breaks=seq(0,0.15,0.01)) +
  scale_y_continuous(breaks=seq(0,400,25)) +
  labs(x = "Variance", y = "Total Sample Size") + theme_bw()


I did the same calculations for the or_meta data frame. Variance for odds ratio is calculated based on the sample size in each of the four cells (treatment-outcome1 (a), treatment-outcome2 (b), control-outcome1 (c), and control-outcome2 (d)):

vi = (1/na) + (1/nb) + (1/nc) + (1/nd)

So again, larger studies will have smaller variances, though imbalance among those four cells (for instance, having a lot more people in one group than another) will inflate variance somewhat. We can see that for study 045, which has a lot more people in the control group than treatment group.

or_meta
##     id study         author_year  tg  cg  tn  cn      yi     vi
## 1  001     1       Bruschke 1999  58  49  72  62  0.0945 0.1860
## 2  003     5    Finkelstein 1995  26  39  60  90  0.0000 0.1131
## 3  005     1           Ruva 2007  67  22 138 138  1.6046 0.0831
## 4  005     2           Ruva 2007  90  50 140 142  1.1976 0.0620
## 5  011     1       Freedman 1996  36  12  99  54  0.6931 0.1508
## 6  016     1         Keelen 1979  37  33 120 120  0.1615 0.0809
## 7  025     1          Davis 1986  17  19  60  52 -0.3759 0.1650
## 8  025     2          Davis 1986  17  17  55  57  0.0513 0.1690
## 9  035     1 Padawer-Singer 1974  47  33  60  60  1.0845 0.1655
## 10 039     1      Eimermann 1971  15  11  40  44  0.5878 0.2279
## 11 045     1        Jacquin 2001 133 207 136 228  1.5035 0.3933
## 12 064     1           Ruva 2006  68  29  87  83  1.8968 0.1203
## 13 064     2           Ruva 2006  53  44  74  73  0.5089 0.1237
ggplot(or_meta, aes(x=vi, y=(tn+cn))) + geom_point() +
  geom_text(aes(label=id),hjust=0, vjust=0) +
  scale_x_continuous(breaks=seq(0,0.4,0.05)) +
  scale_y_continuous(breaks=seq(0,375,25)) +
  labs(x = "Variance", y = "Total Sample Size") + theme_bw()


Tune in tomorrow, when I'll talk about weights for meta-analysis, which is based on this variance calculation!