Data Subset Prep

I’m going to focus on the World Values Survey because it’s awesome. You can download it in many formats, including an rdata file. I prepared a subset for you all.

It includes

• some basic demographics
• a question about how much time one spends thinking about meaning and purpose of life
• ratings on importance of various things in life (e.g., family, friends, work)

setwd('C:/Users/jackflibb/Documents/UO/rClub/611_dplyr_wvs_etc')
load("WorldValuesSurvey.rdata")
ls()

##  [1] "degree"                                       
##  [2] "df"                                           
##  [3] "DV"                                           
##  [4] "gender"                                       
##  [5] "height"                                       
##  [6] "IQ"                                           
##  [7] "n"                                            
##  [8] "RT1"                                          
##  [9] "RT2"                                          
## [10] "WorldValuesSurvey-Wave6-2010-2014_v2014-11-07"

wvs<-`WorldValuesSurvey-Wave6-2010-2014_v2014-11-07`
dim(wvs)

## [1] 85070   387

countryNames<-cbind(names(attr(wvs,"label.table")[[2]]),attr(wvs,"label.table")[[2]])

library(dplyr)

smallerData<-wvs %>%
transmute(
    countryCode=V2,
    sex=V240, #1 male 2 female
    age=V242,
    nationalIncomeLadder=V239,
    meaningInLife=V143,
    imptcOfFamily=V4,
    imptcOfFriends=V5,
    imptcOfLeisure=V6,
    imptcOfPolitics=V7,
    imptcOfWork=V8,
    imptcOfReligion=V9)

write.csv(smallerData,'wvs_small.csv',row.names=F)
write.csv(countryNames,'country_codes.csv',row.names=F)

Don’t worry if you don’t understand some of the syntax yet – I just wanted to show you how I prepare the dataset.

Now we’re going to read in the datasets that you downloaded. Change the working directory to wherever you put the downloaded data files.

setwd('C:/Users/jackflibb/Documents/UO/rClub/611_dplyr_wvs_etc') # Change this
wvs<-read.csv('wvs_small.csv')
countryNames<-read.csv('country_codes.csv',stringsAsFactors=F)

head(wvs)

##   countryCode sex age nationalIncomeLadder meaningInLife imptcOfFamily
## 1          12   1  21                    5             2             1
## 2          12   2  24                    6             1             1
## 3          12   2  26                    6             1             1
## 4          12   2  28                    5             2             1
## 5          12   2  35                    7             2             1
## 6          12   1  36                    5             2             1
##   imptcOfFriends imptcOfLeisure imptcOfPolitics imptcOfWork
## 1              1              1              -2           1
## 2              2              3               4           2
## 3              3              2               4           2
## 4              1              3               4           3
## 5              1              1               2           1
## 6              2              2               2           4
##   imptcOfReligion
## 1               1
## 2               2
## 3               1
## 4               1
## 5               1
## 6               1

Using the chaining operator

As you may have noticed, we have a lot of missing data coded unhelpfully as negative values. The first thing I’ll show you is how to easily recode them to NA using the dplyr function mutate_each. Before I do that, let me introduce the new operator introduced by dplyr: %>%. It’s weird at first, but once you get used to it, it’s great.

Basically, %>% sends the output data frame of one function to the first argument of the next function. For example:

library(dplyr)

wvs %>% names()

##  [1] "countryCode"          "sex"                  "age"                 
##  [4] "nationalIncomeLadder" "meaningInLife"        "imptcOfFamily"       
##  [7] "imptcOfFriends"       "imptcOfLeisure"       "imptcOfPolitics"     
## [10] "imptcOfWork"          "imptcOfReligion"

We can easily add another step to the chain:

wvs %>% names() %>% paste('is a variable in the WVS data')

##  [1] "countryCode is a variable in the WVS data"         
##  [2] "sex is a variable in the WVS data"                 
##  [3] "age is a variable in the WVS data"                 
##  [4] "nationalIncomeLadder is a variable in the WVS data"
##  [5] "meaningInLife is a variable in the WVS data"       
##  [6] "imptcOfFamily is a variable in the WVS data"       
##  [7] "imptcOfFriends is a variable in the WVS data"      
##  [8] "imptcOfLeisure is a variable in the WVS data"      
##  [9] "imptcOfPolitics is a variable in the WVS data"     
## [10] "imptcOfWork is a variable in the WVS data"         
## [11] "imptcOfReligion is a variable in the WVS data"

This isn’t very useful but, it’s a little more interpretable than:

paste(names(wvs),'is a variable in the WVS data')

This will make things a lot more interpretable as we produce complex chains of transformations.

Code missing values using mutate_each

Okay, let’s set the negative values to NA. We’ll use mutate_each, which takes a data frame and applies the functions given to it by the function funs. So you could divide everything in half by doing yourData %>% mutate_each(funs(half = . / 2)).

wvs_NAs<-wvs %>%
mutate_each(funs(ifelse(.<=-1,NA,.)))

Basic summary stats with summarise, summarise_each, select, and filter

Remember describe from last week? The function summarise allows you to easily make your own sort of describe, tailored exactly to your data set.

Let’s look at the structure of our WVS data to decide what we want to do.

str(wvs_NAs)

## 'data.frame':    85070 obs. of  11 variables:
##  $ countryCode         : int  12 12 12 12 12 12 12 12 12 12 ...
##  $ sex                 : int  1 2 2 2 2 1 1 2 1 2 ...
##  $ age                 : int  21 24 26 28 35 36 41 44 59 72 ...
##  $ nationalIncomeLadder: int  5 6 6 5 7 5 7 2 5 4 ...
##  $ meaningInLife       : int  2 1 1 2 2 2 2 2 2 2 ...
##  $ imptcOfFamily       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ imptcOfFriends      : int  1 2 3 1 1 2 1 1 1 1 ...
##  $ imptcOfLeisure      : int  1 3 2 3 1 2 1 1 1 1 ...
##  $ imptcOfPolitics     : int  NA 4 4 4 2 2 1 1 2 2 ...
##  $ imptcOfWork         : int  1 2 2 3 1 4 1 2 2 1 ...
##  $ imptcOfReligion     : int  1 2 1 1 1 1 1 2 2 1 ...

Looks like the demographics should be summarized differently than the ratings.

categorical_summary <- wvs_NAs %>%
summarise(
    N=n(), #a dplyr helper function
    N_countries=length(unique(countryCode)),
    Prop_Male=sum(sex==1,na.rm=T)/N # WOAH, you can use variables created within this function call to build other variables in the same funciton call!!
    )
categorical_summary

##       N N_countries Prop_Male
## 1 85070          59 0.4766428

Let’s look at the numbers for women only…

categorical_summary_F <- wvs_NAs %>%
filter(sex==2) %>% 
summarise(
    N=n(),
    N_countries=length(unique(countryCode)),
    Prop_Male=sum(sex==1,na.rm=T)/N 
    )
categorical_summary_F

##       N N_countries Prop_Male
## 1 44434          59         0

We can use summarise_each to take care of the rest. Like mutate_each, summarise_each takes a comma separated list of functions in funs(...), and optionally a specification of which variables/columns should be summarized. This is really easy in dplyr – we’re going to tell it to summarize all variables from age to imptcOfReligion just by passing in age:imptcOfReligion after our list of functions.

continuous_summary <- wvs_NAs %>% 
select(age:imptcOfReligion) %>% # you can use the : operator directly on variable names!
summarise_each( # glorious ..._each functions!
    funs(
        min=min(.,na.rm=T),
        max=max(.,na.rm=T),
        mean=mean(.,na.rm=T),
        median=median(.,na.rm=T),
        sd=sd(.,na.rm=T),
        Prop_missing=sum(is.na(.))/length(.)
        )
    )
continuous_summary

##   age_min nationalIncomeLadder_min meaningInLife_min imptcOfFamily_min
## 1      16                        1                 1                 1
##   imptcOfFriends_min imptcOfLeisure_min imptcOfPolitics_min
## 1                  1                  1                   1
##   imptcOfWork_min imptcOfReligion_min age_max nationalIncomeLadder_max
## 1               1                   1      99                       10
##   meaningInLife_max imptcOfFamily_max imptcOfFriends_max
## 1                 4                 4                  4
##   imptcOfLeisure_max imptcOfPolitics_max imptcOfWork_max
## 1                  4                   4               4
##   imptcOfReligion_max age_mean nationalIncomeLadder_mean
## 1                   4 42.05654                  4.908784
##   meaningInLife_mean imptcOfFamily_mean imptcOfFriends_mean
## 1           1.844841           1.107936            1.684175
##   imptcOfLeisure_mean imptcOfPolitics_mean imptcOfWork_mean
## 1             1.88321             2.612757         1.545598
##   imptcOfReligion_mean age_median nationalIncomeLadder_median
## 1             1.914088         40                           5
##   meaningInLife_median imptcOfFamily_median imptcOfFriends_median
## 1                    2                    1                     2
##   imptcOfLeisure_median imptcOfPolitics_median imptcOfWork_median
## 1                     2                      3                  1
##   imptcOfReligion_median   age_sd nationalIncomeLadder_sd meaningInLife_sd
## 1                      2 16.54851                2.104927        0.8582036
##   imptcOfFamily_sd imptcOfFriends_sd imptcOfLeisure_sd imptcOfPolitics_sd
## 1        0.3758205         0.7397572         0.8341943          0.9822105
##   imptcOfWork_sd imptcOfReligion_sd age_Prop_missing
## 1      0.8116806           1.064065      0.001798519
##   nationalIncomeLadder_Prop_missing meaningInLife_Prop_missing
## 1                        0.03605266                   0.015787
##   imptcOfFamily_Prop_missing imptcOfFriends_Prop_missing
## 1                0.003714588                 0.005442577
##   imptcOfLeisure_Prop_missing imptcOfPolitics_Prop_missing
## 1                  0.01120254                   0.01702128
##   imptcOfWork_Prop_missing imptcOfReligion_Prop_missing
## 1               0.01702128                   0.01498766

Dang, that’s kind of ugly. We can reshape these using the reshape2 package and output it to a table using kable from knitr.

library(knitr)
library(reshape2)

continuous_summary_melted <- melt(continuous_summary)

## No id variables; using all as measure variables

kable(categorical_summary,
    caption='Summary of Categorical Variables')

kable(continuous_summary_melted,
    col.names=c('Statistic','Value'),
    digits=2,
    caption='Summary of Continuous Variables')

continuous_summary_melted %>% 
transmute(
    item=sub('(\\w+?)_(.*)','\\1',variable),
    statistic=sub('(\\w+?)_(.*)','\\2',variable),
    value=value) %>%
arrange(item,statistic) %>%
kable(
    col.names=c('Item','Statistic','Value'),
    digits=2,
    caption='Summary of Continuous Variables')

Get country names using left_join

We also have a file with country names to go with our country codes. We’re going to start doing some analysis within-country, so let’s get those names into the data set. This is what some call a one-to-many merge. That is, we have a data frame with one row per country name, and we want to merge each of these rows with every row in the main data set that has a matching country code.

There are 4 join functions (?left_join for more info), of which left_join is our best choice. According to the docs it “return[s] all rows from x, and all columns from x and y. If there are multiple matches between x and y, all combination of the matches are returned.”

str(countryNames)

## 'data.frame':    195 obs. of  2 variables:
##  $ V1: chr  "Missing; Unknown" "Not asked in survey" "Not applicable" "No answer" ...
##  $ V2: int  -5 -4 -3 -2 -1 4 8 12 16 20 ...

names(wvs_NAs)

##  [1] "countryCode"          "sex"                  "age"                 
##  [4] "nationalIncomeLadder" "meaningInLife"        "imptcOfFamily"       
##  [7] "imptcOfFriends"       "imptcOfLeisure"       "imptcOfPolitics"     
## [10] "imptcOfWork"          "imptcOfReligion"

#We use the same country code variable name, countryCode, as the main data set
names(countryNames)<-c('countryName','countryCode') 

#left_join will automatically find matching column names
wvs_NAs_cnames <- left_join(wvs_NAs,countryNames)

## Joining by: "countryCode"

head(wvs_NAs_cnames) #success!

##   countryCode sex age nationalIncomeLadder meaningInLife imptcOfFamily
## 1          12   1  21                    5             2             1
## 2          12   2  24                    6             1             1
## 3          12   2  26                    6             1             1
## 4          12   2  28                    5             2             1
## 5          12   2  35                    7             2             1
## 6          12   1  36                    5             2             1
##   imptcOfFriends imptcOfLeisure imptcOfPolitics imptcOfWork
## 1              1              1              NA           1
## 2              2              3               4           2
## 3              3              2               4           2
## 4              1              3               4           3
## 5              1              1               2           1
## 6              2              2               2           4
##   imptcOfReligion countryName
## 1               1     Algeria
## 2               2     Algeria
## 3               1     Algeria
## 4               1     Algeria
## 5               1     Algeria
## 6               1     Algeria

Summaries by country

We’re going to use the summarise_each function again, but this time instead of using select before hand, we can give it the columns we want as the second argument. I’ll use the helper function starts_with to get every column that starts with ‘imptcOf’ (our “Importance in life of…” variables).

wvs_NAs_cnames %>%
group_by(countryName) %>%
summarise_each(
    funs(
        mean=mean(.,na.rm=T)),
    starts_with('imptcOf')) %>%
kable(caption='Importance of ___ Item Means by Country')

Histograms of each variable by country

You can apply this same type of logic for running models for each country. We melt the data frame so that we get a row for every variable/value pair indexed by country. This makes it easy for ggplot to give us a a histogram for each variable in the same plot, using facet_wrap. If we didn’t use facet_wrap below, we’d get one histogram with the values of all the variables. Check out the first six, last, and a random sample of six lines of the melted data frame:

m_temp <- melt(wvs_NAs_cnames,id.vars='countryName')
kable(head(m_temp))

kable(tail(m_temp))

kable(sample_n(m_temp,6),row.names=F)

library(ggplot2)
#plotprintlist is a variable that is basically trash, but we could use it to print the plots again if we wanted
plotprintlist <- wvs_NAs_cnames %>%
melt(id.vars='countryName') %>%
group_by(countryName) %>%
do(
     plot=print(
        ggplot(.,aes(x=value))+ # set up the plot with the data frame '.' passed to `do`
            geom_histogram()+ # We want a histogram of the values from column 'value'
            facet_wrap(~variable,scales='free')+ # We want a facet (i.e., subplot) for every unique value of the column 'variable'
            ggtitle(unique(.$countryName)) # give it a title with the country name from the countrName column of '.'
            )
)

All the analysis

We can use the same technique to get a coefficients from a linear model for every country

wvs_NAs_cnames_modulus <- wvs_NAs_cnames %>%
mutate(age_ending = age %% 10, age_ending_9 = age_ending %in% 9)

library(biglm)

## Warning: package 'biglm' was built under R version 3.1.2

## Loading required package: DBI

## Warning: package 'DBI' was built under R version 3.1.2

models <- wvs_NAs_cnames_modulus %>%
group_by(countryName) %>%
do(
     mod=biglm(
            meaningInLife ~ 1 + 
            sex +
            age +
            nationalIncomeLadder +
            imptcOfFamily +
            imptcOfFriends +
            imptcOfLeisure +
            imptcOfPolitics +
            imptcOfWork +
            imptcOfReligion +
            age_ending_9,
            data=.)
     )

modelCoefs <- models %>%
do(
    data.frame(
        country=.$countryName,
        n=.$mod$n,
        ci=summary(.$mod)$mat[,'Coef']-summary(.$mod)$mat[,'(95%'],
        var=names(coef(.$mod)),
        coef=coef(.$mod))
    )

## Warning in rbind_all(out[[1]]): Unequal factor levels: coercing to
## character

kable(modelCoefs[1:20,])

modelCoefs %>%
ggplot(aes(x=reorder(country,n),y=coef))+
geom_point(aes(size=n),alpha=.3,color='blue')+
geom_point(size=1,alpha=1,color='black')+
geom_errorbar(aes(ymin=coef-ci,ymax=coef+ci))+
facet_wrap(~var,scales='free_y',ncol=1)+
geom_hline(xintercept = 0,color='red',size=.5)+
theme(
    axis.text=element_text(size=10),
    axis.text.x=element_text(angle=315,hjust=0,vjust=1),
    panel.background=element_rect(fill='white'),
    strip.text=element_text(size=10))+
labs(x='Country, ordered by sample size',y='Coefficient Estimate',title='How often do you contemplate the meaning of life?')