What is R?

Other Software

PSPP

Pros and Cons

Zoltan Fazekas’s analogy:

Many of you have probably seen Fight Club and this course or working in R might be similar:

In terms of R:

R and R Studio. (available for mac, linux and windows)

Download links R R Studio

IDE (Integrated Development Environments). There are others: emacs and vi and sublime text (even steeper learning curve)

Scripting:

So:

#writewhateverhere

Does nothing

Very simple commands:

2 + 2
## [1] 4
2*4 + 2
## [1] 10
2*(4+2)
## [1] 12
print('Hello World')
## [1] "Hello World"

You can use both double or simple quotes, just keep it consistent

Check logical arguments:

26 == (3 * 13) - 13
## [1] TRUE
26 == (3 * 13) - 12
## [1] FALSE

Working directory

Get using getw:

getwd()
## [1] "/Users/levi/data"

Set using setwd()

setwd('~/data/') # on Windows, it is 'c:/data/'
getwd()
## [1] "/Users/levi/data"

Objects

R -Object-oriented programming language

Operands <- and -> Operation

y <- -10
y
## [1] -10
print(y) #more elegant
## [1] -10

Also works with logical statements

x <- T
x
## [1] TRUE

R is not going to stop you from running invalid calculations (but it still won’t perform for example calculations with strings)

z <- x + y
z
## [1] -9

Beware of rewriting objects, you can’t undo it

Mathematical operands + - * / ^

6^2
## [1] 36

also < <= > >= == (equals) != (does not equal), & (intersection), | (union)

Vectors

Defined with the c command, like this:

numericVector <- c(1, 2, 3, 4, 5) 
print(numericVector)
## [1] 1 2 3 4 5
numericVector <- c(1:5) #or 
print(numericVector)
## [1] 1 2 3 4 5
charvector <- c('Bob', 'Jane', 'Jack')
print(charvector)
## [1] "Bob"  "Jane" "Jack"
logicalvector <- c(T, F, T) ## T = TRUE, F = FALSE)
print(logicalvector)
## [1]  TRUE FALSE  TRUE

ls command - lists all objects

ls()
## [1] "charvector"    "logicalvector" "numericVector" "x"            
## [5] "y"             "z"

rm command - removes objects

rm()

You can, of course, do math with vectors (only your algebra skills are the limit)

Protip: You can also print your object right away by putting it in brackets like this:

(numericVector <- c(1, 2, 3, 4, 5))   #look! No separate print command
## [1] 1 2 3 4 5

Matrices

Order is rows by columns - ALWAYS!

trialMat <- matrix(1:20, nrow = 2, ncol = 10, byrow = F)
trialMat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    3    5    7    9   11   13   15   17    19
## [2,]    2    4    6    8   10   12   14   16   18    20
trialMat <- matrix(1:20, nrow = 2, ncol = 10, byrow = T)
trialMat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    2    3    4    5    6    7    8    9    10
## [2,]   11   12   13   14   15   16   17   18   19    20
trialMat[1, ]  # get the row with all columns if coordinate 2 blank
##  [1]  1  2  3  4  5  6  7  8  9 10
trialMat[ , 5]  # get the column with all rows if coordinate 1 blank
## [1]  5 15
trialMat[2, 2]  # get exact elements
## [1] 12

You can also transpose a matrix like this:

trialMat <- matrix(1:20, nrow = 2, ncol = 10, byrow = T)
tranMat <- t(trialMat)
print(tranMat)
##       [,1] [,2]
##  [1,]    1   11
##  [2,]    2   12
##  [3,]    3   13
##  [4,]    4   14
##  [5,]    5   15
##  [6,]    6   16
##  [7,]    7   17
##  [8,]    8   18
##  [9,]    9   19
## [10,]   10   20
propMat <- t(tranMat)
print(propMat)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    2    3    4    5    6    7    8    9    10
## [2,]   11   12   13   14   15   16   17   18   19    20

Fun with as.vector function

trialMat <- matrix(1:20, nrow = 2, ncol = 10, byrow = T)
as.vector (trialMat)
##  [1]  1 11  2 12  3 13  4 14  5 15  6 16  7 17  8 18  9 19 10 20
as.vector (t(trialMat)) # transposed
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

Lost?

Getting help - just use question mark - like

?sd

Need help with syntax? Check arguments using “args” like this

args(sd)
## function (x, na.rm = FALSE) 
## NULL

Continuting with matrices

a <- c(1:25) # Lets make some data
sd(a) # Lets check the standard deviation
## [1] 7.359801

Missing values

b <- c(1, 2, 3, NA, -5) # Make some more data.  Note NA!
sd(b)
## [1] NA
sd(b, na.rm = T) # The na.rm argument - missing not (always) automatically removed by R
## [1] 3.593976

Histogram

Lets make some more data

d <- rnorm(300) #rnorm produces random normal numbers with mean of 0 and sd of 1
hist(d)  # Shows us a histogram

Packages

Package ‘foreign’ = awesome. It reads Data Stored by Minitab, S, SAS, SPSS, Stata, Systat, Weka, dBase…you name it (though to read Stata 13 or newer files, you need a different package: readstata13)

install.packages('foreign') # This is how you download a package
library(foreign) # This is how you indicate to R it's used in your script

To update the packages

update.packages()

What packages do I have installed?

(.packages(all.available=TRUE))
##   [1] "abind"            "acepack"          "AER"             
##   [4] "alr4"             "Amelia"           "antiword"        
##   [7] "arm"              "assertthat"       "backports"       
##  [10] "base"             "base64enc"        "bdsmatrix"       
##  [13] "betareg"          "BH"               "bindr"           
##  [16] "bindrcpp"         "bitops"           "boot"            
##  [19] "broom"            "btergm"           "car"             
##  [22] "carData"          "caTools"          "cellranger"      
##  [25] "checkmate"        "class"            "cli"             
##  [28] "cluster"          "coda"             "codetools"       
##  [31] "colorspace"       "compiler"         "config"          
##  [34] "covr"             "coxme"            "crayon"          
##  [37] "crosstalk"        "curl"             "DAMisc"          
##  [40] "data.table"       "datasets"         "DEoptimR"        
##  [43] "descr"            "digest"           "doParallel"      
##  [46] "dotwhisker"       "dplyr"            "effects"         
##  [49] "ergm"             "ergm.count"       "estimability"    
##  [52] "evaluate"         "fansi"            "filehash"        
##  [55] "flexmix"          "forcats"          "foreach"         
##  [58] "foreign"          "Formula"          "gdata"           
##  [61] "geepack"          "ggplot2"          "ggstance"        
##  [64] "glue"             "GPArotation"      "gplots"          
##  [67] "graphics"         "grDevices"        "grid"            
##  [70] "gridExtra"        "gsl"              "gtable"          
##  [73] "gtools"           "haven"            "heplots"         
##  [76] "highr"            "Hmisc"            "hms"             
##  [79] "htmlTable"        "htmltools"        "htmlwidgets"     
##  [82] "httpuv"           "httr"             "igraph"          
##  [85] "interactionTest"  "interplot"        "iterators"       
##  [88] "itertools"        "jsonlite"         "keras"           
##  [91] "kerasformula"     "KernSmooth"       "knitr"           
##  [94] "labeling"         "later"            "lattice"         
##  [97] "latticeExtra"     "lavaan"           "lazyeval"        
## [100] "leaps"            "lme4"             "lmtest"          
## [103] "lpSolve"          "magrittr"         "manipulateWidget"
## [106] "maptools"         "markdown"         "MASS"            
## [109] "Matching"         "MatchIt"          "Matrix"          
## [112] "matrixcalc"       "MatrixModels"     "maxLik"          
## [115] "mcmc"             "MCMCpack"         "mediation"       
## [118] "methods"          "mfx"              "mgcv"            
## [121] "mi"               "mime"             "miniUI"          
## [124] "minqa"            "miscTools"        "mnormt"          
## [127] "modeltools"       "munsell"          "mvtnorm"         
## [130] "network"          "networkDynamic"   "nFactors"        
## [133] "nlme"             "nloptr"           "NLP"             
## [136] "nnet"             "nose"             "numDeriv"        
## [139] "openssl"          "openxlsx"         "ordinal"         
## [142] "parallel"         "pbivnorm"         "pbkrtest"        
## [145] "pdftools"         "pillar"           "pkgconfig"       
## [148] "plogr"            "plyr"             "poLCA"           
## [151] "praise"           "prettyunits"      "processx"        
## [154] "progress"         "promises"         "ps"              
## [157] "pscl"             "psych"            "purrr"           
## [160] "QRM"              "quantreg"         "R6"              
## [163] "RColorBrewer"     "Rcpp"             "RcppArmadillo"   
## [166] "RcppEigen"        "Rcsdp"            "RCurl"           
## [169] "readr"            "readstata13"      "readxl"          
## [172] "rematch"          "remotes"          "reshape2"        
## [175] "reticulate"       "rex"              "rgl"             
## [178] "rio"              "rlang"            "rmarkdown"       
## [181] "robustbase"       "ROCR"             "rpart"           
## [184] "rprojroot"        "RSiena"           "rstudioapi"      
## [187] "sandwich"         "scales"           "scatterplot3d"   
## [190] "sem"              "shiny"            "slam"            
## [193] "sm"               "sna"              "SnowballC"       
## [196] "sourcetools"      "sp"               "SparseM"         
## [199] "spatial"          "speedglm"         "splines"         
## [202] "statnet"          "statnet.common"   "stats"           
## [205] "stats4"           "stringi"          "stringr"         
## [208] "SuppDists"        "survey"           "survival"        
## [211] "swirl"            "sys"              "tcltk"           
## [214] "tensorflow"       "tergm"            "testthat"        
## [217] "texreg"           "tfestimators"     "tfruns"          
## [220] "tibble"           "tidyr"            "tidyselect"      
## [223] "timeDate"         "timeSeries"       "tinytex"         
## [226] "tm"               "tools"            "trust"           
## [229] "ucminf"           "utf8"             "utils"           
## [232] "VGAM"             "viridis"          "viridisLite"     
## [235] "webshot"          "whisker"          "withr"           
## [238] "xergm.common"     "xfun"             "XML"             
## [241] "xml2"             "xtable"           "yaml"            
## [244] "zeallot"          "Zelig"            "zip"             
## [247] "zoo"

Lists

Are kind of an abstract way of storing information in Objects

simpleList <- list(course = "Intro to R", start = 17, credit = 2)
print (simpleList)
## $course
## [1] "Intro to R"
## 
## $start
## [1] 17
## 
## $credit
## [1] 2
compList <- list(course = c("Intro to R", "Intro to Stata"), start = c(17, 
    17), credit = c(2, NA))
print (compList)
## $course
## [1] "Intro to R"     "Intro to Stata"
## 
## $start
## [1] 17 17
## 
## $credit
## [1]  2 NA

Getting info from the list

compList[1]
## $course
## [1] "Intro to R"     "Intro to Stata"
compList[[1]]
## [1] "Intro to R"     "Intro to Stata"
compList[[1]][1]
## [1] "Intro to R"
max(compList$start) # ?max to see what max does
## [1] 17

Data Frames

(newList <- list(x = c(2, 5), y = c(6, 8.1)))
## $x
## [1] 2 5
## 
## $y
## [1] 6.0 8.1
newList$x 
## [1] 2 5
(newData <- as.data.frame(newList))
##   x   y
## 1 2 6.0
## 2 5 8.1
newData[[1]] == newData[1, ]  ## you can access elements intuitively, why false for one element?
##      x     y
## 1 TRUE FALSE
newData[[1]] == newData[, 1]
## [1] TRUE TRUE
newData$x
## [1] 2 5
mean(newData$x)  ## you can see that we are asking R to calculate the mean of the first variable
## [1] 3.5
mean(newData[, 1])  ## this will do the same
## [1] 3.5
mean(newData[1, ])  ## remember, these are data frames and some functions are applied to columns, not rows.
## Warning in mean.default(newData[1, ]): argument is not numeric or logical:
## returning NA
## [1] NA

Lets make up some data

x1 <- 1:25
x2 <- 51:75
exData <- as.data.frame(cbind(x1, x2)) # cbind - binds columns together.  (see also rbind)

Ways to look at your data

head(exData)  ## display first 6 rows
##   x1 x2
## 1  1 51
## 2  2 52
## 3  3 53
## 4  4 54
## 5  5 55
## 6  6 56
tail(exData)  ## display last 6 rows
##    x1 x2
## 20 20 70
## 21 21 71
## 22 22 72
## 23 23 73
## 24 24 74
## 25 25 75
exData[10:13, ]  ## display rows of your preference
##    x1 x2
## 10 10 60
## 11 11 61
## 12 12 62
## 13 13 63
dim(exData)  ## again, know the dimensions! (That is what the dim commands does)
## [1] 25  2
str(exData)  ## very useful to get an idea what you are up against, but also to check whether you have the data frame and format with the variables you wanted to have (Probably most important R command.)
## 'data.frame':    25 obs. of  2 variables:
##  $ x1: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ x2: int  51 52 53 54 55 56 57 58 59 60 ...
length(exData$x1)  ## Observations 25 we have, as Yoda would articulate it. Suppose your observations are individuals, and you want a simple row name: Individual and the number. Let us generate a character vector of length 25 to do this
## [1] 25
(rNames <- paste("Individual", 1:25, sep = " "))  
##  [1] "Individual 1"  "Individual 2"  "Individual 3"  "Individual 4" 
##  [5] "Individual 5"  "Individual 6"  "Individual 7"  "Individual 8" 
##  [9] "Individual 9"  "Individual 10" "Individual 11" "Individual 12"
## [13] "Individual 13" "Individual 14" "Individual 15" "Individual 16"
## [17] "Individual 17" "Individual 18" "Individual 19" "Individual 20"
## [21] "Individual 21" "Individual 22" "Individual 23" "Individual 24"
## [25] "Individual 25"

or

rNames <- paste("Individual", 1:length(exData$x1), sep = " ")  ## why is this better?
rownames(exData) <- rNames  ## as you can see, assign our new character vector to the row names of our data. You can use return values from a function on both sides of an assignment
exData[1:3, ]
##              x1 x2
## Individual 1  1 51
## Individual 2  2 52
## Individual 3  3 53
dim(exData)  ## Note, length(data frame) is not the way to go. length is good for vectors
## [1] 25  2
length(exData)
## [1] 2
length(exData$x2)
## [1] 25
names(exData)  ## use names to get your variable names. How many should there be? Again, you can store it out if you wish for further reference, or you can change it. It is a character vector!
## [1] "x1" "x2"
names(exData)[1] <- "var1"  ## it gives an error if it is not a character you are supplying
names(exData)
## [1] "var1" "x2"
names(exData) <- c("newVar1", "newVar2")
names(exData)  ## always know how many values you should supply
## [1] "newVar1" "newVar2"
summary(exData)
##     newVar1      newVar2  
##  Min.   : 1   Min.   :51  
##  1st Qu.: 7   1st Qu.:57  
##  Median :13   Median :63  
##  Mean   :13   Mean   :63  
##  3rd Qu.:19   3rd Qu.:69  
##  Max.   :25   Max.   :75
str(exData)
## 'data.frame':    25 obs. of  2 variables:
##  $ newVar1: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ newVar2: int  51 52 53 54 55 56 57 58 59 60 ...

Importing data

  1. Set working directory
setwd("~/data/") #or whatever you want to use

2.Put the data from the elearning in that folder.** (All of them.)

  1. Now make sure you have the foreign package installed
## install.packages('foreign', dependencies = T) ## uncomment this line if
## you are working on your personal computer and do not have the foreign
## package
library(foreign)
  1. Now you can import data

SPSS

SPSS <- read.spss("example.por")  ## very basic, and accordingly not working
str(SPSS)  ## again, this is not what you expect. It is a list.
## List of 10
##  $ VERSION : chr [1:2208] "ANES2008TS_VERSION:20120217" "ANES2008TS_VERSION:20120217" "ANES2008TS_VERSION:20120217" "ANES2008TS_VERSION:20120217" ...
##  $ V080001 : num [1:2208] 1 2 3 4 5 6 7 8 9 10 ...
##  $ V080101 : num [1:2208] 2 2.16 2.46 2.16 2.93 ...
##  $ V080101A: num [1:2208] 175097 189642 215132 189642 256749 ...
##  $ V080102 : Factor w/ 807 levels "0. No Post-election interview",..: 725 1 777 755 803 1 770 743 771 295 ...
##  $ V080102A: Factor w/ 856 levels "0. No Post-election interview",..: 774 1 826 804 852 1 819 792 820 333 ...
##  $ V080103 : num [1:2208] 67986 67986 67986 67986 67986 ...
##  $ V081001 : Factor w/ 2 levels "0. Pre-election only",..: 2 1 2 2 2 1 2 2 2 2 ...
##  $ V081101 : Factor w/ 2 levels "1. Male respondent selected",..: 1 2 2 1 1 2 2 1 2 1 ...
##  $ V081102 : Factor w/ 9 levels "-9. Refused in household listing",..: 6 3 3 3 3 3 3 3 3 4 ...
##  - attr(*, "label.table")=List of 10
##   ..$ VERSION : NULL
##   ..$ V080001 : NULL
##   ..$ V080101 : NULL
##   ..$ V080101A: NULL
##   ..$ V080102 : Named num 0
##   .. ..- attr(*, "names")= chr "0. No Post-election interview"
##   ..$ V080102A: Named num 0
##   .. ..- attr(*, "names")= chr "0. No Post-election interview"
##   ..$ V080103 : NULL
##   ..$ V081001 : Named num [1:2] 1 0
##   .. ..- attr(*, "names")= chr [1:2] "1. Pre-election and Post-election" "0. Pre-election only"
##   ..$ V081101 : Named num [1:2] 2 1
##   .. ..- attr(*, "names")= chr [1:2] "2. Female respondent selected" "1. Male respondent selected"
##   ..$ V081102 : Named num [1:9] 7 6 5 4 3 2 1 -4 -9
##   .. ..- attr(*, "names")= chr [1:9] "7. White, black and another race" "6. Black and another race" "5. White and another race" "4. Other race" ...
##  - attr(*, "variable.labels")= Named chr [1:10] " " "ID.1. CASE ID" "WT.1. PRE CROSS-SECTION SAMPLE WEIGHT - POST-STRAT, centered" "WT.1. PRE CROSS-SECTION SAMPLE WEIGHT - POST-STRATIFIED" ...
##   ..- attr(*, "names")= chr [1:10] "VERSION" "V080001" "V080101" "V080101A" ...
##  - attr(*, "missings")=List of 10
##   ..$ VERSION :List of 1
##   .. ..$ type: chr "none"
##   ..$ V080001 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V080101 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V080101A:List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V080102 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V080102A:List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V080103 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V081001 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V081101 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
##   ..$ V081102 :List of 2
##   .. ..$ type : chr "low"
##   .. ..$ value: num -1
is.data.frame(SPSS) ## but we want data frames
## [1] FALSE
SPSS <- read.spss("example.por", to.data.frame = T)  ## it is closer to what we want
summary(SPSS[8:10])  ## but those labels?! Maybe better without
##                               V081001    
##  0. Pre-election only             : 212  
##  1. Pre-election and Post-election:1995  
##  NA's                             :   1  
##                                          
##                                          
##                                          
##                                          
##                           V081101                          V081102    
##  1. Male respondent selected  : 941   1. White                 :1355  
##  2. Female respondent selected:1266   2. Black/African-American: 565  
##  NA's                         :   1   4. Other race            : 253  
##                                       5. White and another race:  16  
##                                       6. Black and another race:   6  
##                                       (Other)                  :   2  
##                                       NA's                     :  11
#hist(SPSS[, 10]) You can try what this does

SPSS <- read.spss("example.por", to.data.frame = T, use.value.labels = F)  ##
summary(SPSS[8:10])  ## now working with this and the codebook will be better
##     V081001          V081101         V081102     
##  Min.   :0.0000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:1.0000   1st Qu.:1.000   1st Qu.:1.000  
##  Median :1.0000   Median :2.000   Median :1.000  
##  Mean   :0.9039   Mean   :1.574   Mean   :1.651  
##  3rd Qu.:1.0000   3rd Qu.:2.000   3rd Qu.:2.000  
##  Max.   :1.0000   Max.   :2.000   Max.   :7.000  
##  NA's   :1        NA's   :1       NA's   :11
hist(SPSS[, 10])

Stata

STATA <- read.dta("example.dta")  ## try Stata files
## Warning in read.dta("example.dta"): value labels ('V080001_') for 'V080001'
## are missing
str(STATA)  ## here, we did better by default, as it is a data.frame
## 'data.frame':    2322 obs. of  10 variables:
##  $ Version : chr  "ANES2008TS_VERSION:20120217" "ANES2008TS_VERSION:20120217" "ANES2008TS_VERSION:20120217" "ANES2008TS_VERSION:20120217" ...
##  $ V080001 : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ V080101 : num  2 2.16 2.46 2.16 2.93 ...
##  $ V080101a: num  175097 189642 215132 189642 256749 ...
##  $ V080102 : num  1.92 0 2.41 2.13 3 ...
##  $ V080102a: num  186174 0 232918 206027 290691 ...
##  $ V080103 : num  67986 67986 67986 67986 67986 ...
##  $ V081001 : Factor w/ 2 levels "0. Pre-election only",..: 2 1 2 2 2 1 2 2 2 2 ...
##  $ V081101 : Factor w/ 2 levels "1. Male respondent selected",..: 1 2 2 1 1 2 2 1 2 1 ...
##  $ V081102 : Factor w/ 9 levels "-9. Refused in household listing",..: 6 3 3 3 3 3 3 3 3 4 ...
##  - attr(*, "datalabel")= chr ""
##  - attr(*, "time.stamp")= chr "12 Feb 2013 14:56"
##  - attr(*, "formats")= chr  "%27s" "%9.0g" "%10.0g" "%10.0g" ...
##  - attr(*, "types")= int  27 254 255 255 255 255 255 251 251 254
##  - attr(*, "val.labels")= chr  "" "V080001_" "" "" ...
##  - attr(*, "var.labels")= chr  "" "ID.1. CASE ID" "WT.1. PRE CROSS-SECTION SAMPLE WEIGHT - POST-STRAT, centered" "WT.1. PRE CROSS-SECTION SAMPLE WEIGHT - POST-STRATIFIED" ...
##  - attr(*, "version")= int 12
##  - attr(*, "label.table")=List of 5
##   ..$ V080102_: Named int 0
##   .. ..- attr(*, "names")= chr "0. No Post-election interview"
##   ..$ V080102a: Named int 0
##   .. ..- attr(*, "names")= chr "0. No Post-election interview"
##   ..$ V081001_: Named int  0 1
##   .. ..- attr(*, "names")= chr  "0. Pre-election only" "1. Pre-election and Post-election"
##   ..$ V081101_: Named int  1 2
##   .. ..- attr(*, "names")= chr  "1. Male respondent selected" "2. Female respondent selected"
##   ..$ V081102_: Named int  -9 -4 1 2 3 4 5 6 7
##   .. ..- attr(*, "names")= chr  "-9. Refused in household listing" "-4. NA (blank recorded)" "1. White" "2. Black/African-American" ...
#hist(STATA[, 10])  ## we still have labels, and some variables are factors. It is not bad, just confusing

STATA <- read.dta("example.dta", convert.factors = F)  ## although with a different function attribute, we fix that as well
hist(STATA[, 10])  ## but of course, -9 was supposed to be a missing value, but it was not coded as system missing in STATA, whereas if you remember, in the SPSS version of this data this was done (in the original data)

When working with data, we need to make sure missing values are coded correctly. More on this later.

# Aditionaly you might encounter some data in stata 13 and more recent formats. Unfortunately, foreign package cannot load 
# this data so you will need to install additional packages. The package you are looking for is readstata13 and the function is very similar to the foreign function, read.dta13.

#install.packages("readstata13")
library(readstata13)

STATA13 <- read.dta13("example13.dta")  ## try Stata 13 files
str(STATA13)  
## 'data.frame':    6000 obs. of  9 variables:
##  $ x1: num  -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ x2: num  -0.475 -0.425 -0.375 -0.325 -0.275 ...
##  $ x3: num  -0.1252 -0.0597 -0.4811 0.1341 -0.2678 ...
##  $ x4: num  -0.484 0.1251 0.1075 0.1713 -0.0225 ...
##  $ x5: num  -0.4765 -0.5894 -0.0468 -0.1133 -0.1585 ...
##  $ x6: num  -0.03671 -0.0015 0.03864 -0.00021 0.01257 ...
##  $ id: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ y : num  -0.9082 -0.0575 -1.0565 0.2599 0.5356 ...
##  $ z : int  0 0 0 1 1 0 0 1 0 0 ...
##  - attr(*, "datalabel")= chr ""
##  - attr(*, "time.stamp")= chr " 3 Mar 2013 13:07"
##  - attr(*, "formats")= chr  "%9.0g" "%9.0g" "%9.0g" "%9.0g" ...
##  - attr(*, "types")= int  65527 65527 65527 65527 65527 65527 65528 65527 65530
##  - attr(*, "val.labels")= Named chr  "" "" "" "" ...
##   ..- attr(*, "names")= chr  "" "" "" "" ...
##  - attr(*, "var.labels")= chr  "" "" "" "" ...
##  - attr(*, "version")= int 117
##  - attr(*, "label.table")= list()
##  - attr(*, "expansion.fields")=List of 1
##   ..$ : chr  "_dta" "iis" "id"
##  - attr(*, "byteorder")= chr "LSF"
##  - attr(*, "orig.dim")= int  6000 9

When working with data, we need to make sure missing values are coded correctly. More on this later.

How about opening up those text files? R can actually do this without a package.

(TABDATA <- read.table("example.txt"))  ## Something is wrong. It seems that we already have the variable names in the data, but R does not know that
##         V1    V2    V3     V4    V5    V6   V7    V8
## 1  country polar   vol nopart  disp   GDP demo nogov
## 2       AT   1.4 12.39   3.41  2.47 28900   56     8
## 3       BE  2.97 13.93   9.05  4.46 27200   80     9
## 4       BG  6.43 42.09   2.52  7.37  2600   20     9
## 5      CYP 33.47  6.89   3.52  2.82 16500   50     4
## 6       CZ  9.97 29.29   3.71  5.21  8100   20    12
## 7       DK  3.74 10.95   4.71   1.5 35600   80     9
## 8       EE   8.1 36.37    5.5  4.62  7900   18    12
## 9      FIN  7.28  9.22   5.15  3.31 30900   80     6
## 10      FR  7.29 18.96   3.54 18.37 25700   80     8
## 11     GER  4.07   8.1   3.31  3.28 27300   62     5
## 12     GRE 11.66  9.36   2.21  7.07 16400   36     7
## 13     HUN  5.75 23.61   4.08  9.34  6600   20     9
## 14     IRE  5.91 10.28   2.99  6.03 35100   80     8
## 15     ITA  6.59 30.92   6.07  5.62 21700   63    12
## 16     LAT  8.57 53.16   5.49  3.93  6800   17    16
## 17     LIT   2.5 48.72   3.52  8.34  6300   18    16
## 18     LUX  3.87  8.43   4.34  4.15 61500   80     5
## 19      MT    23   2.8   1.99  0.66 11700   35     6
## 20      NL  1.85 21.31   4.81  0.95 29200   80     7
## 21      PL   0.7 43.15   2.95  8.38  6400   20    21
## 22      PT 11.75 11.69   2.62  5.29 12400   34     6
## 23      RO  4.59 36.08   3.53  7.01  2900   20    16
## 24     SVK  6.03 40.62   4.76  5.64  6200   20    13
## 25     SLO  3.66 30.89   4.55   3.8 14400   19    12
## 26     SPA  6.03     8   2.48  6.33 17800   33     6
## 27     SWE     3 15.06   4.29  1.78 35300   80     6
## 28      UK   1.4   7.7   2.12 16.16 31400   80     5
(TABDATA <- read.table("example.txt", header = T))  ## now, the variable names are in their right place
##    country polar   vol nopart  disp   GDP demo nogov
## 1       AT  1.40 12.39   3.41  2.47 28900   56     8
## 2       BE  2.97 13.93   9.05  4.46 27200   80     9
## 3       BG  6.43 42.09   2.52  7.37  2600   20     9
## 4      CYP 33.47  6.89   3.52  2.82 16500   50     4
## 5       CZ  9.97 29.29   3.71  5.21  8100   20    12
## 6       DK  3.74 10.95   4.71  1.50 35600   80     9
## 7       EE  8.10 36.37   5.50  4.62  7900   18    12
## 8      FIN  7.28  9.22   5.15  3.31 30900   80     6
## 9       FR  7.29 18.96   3.54 18.37 25700   80     8
## 10     GER  4.07  8.10   3.31  3.28 27300   62     5
## 11     GRE 11.66  9.36   2.21  7.07 16400   36     7
## 12     HUN  5.75 23.61   4.08  9.34  6600   20     9
## 13     IRE  5.91 10.28   2.99  6.03 35100   80     8
## 14     ITA  6.59 30.92   6.07  5.62 21700   63    12
## 15     LAT  8.57 53.16   5.49  3.93  6800   17    16
## 16     LIT  2.50 48.72   3.52  8.34  6300   18    16
## 17     LUX  3.87  8.43   4.34  4.15 61500   80     5
## 18      MT 23.00  2.80   1.99  0.66 11700   35     6
## 19      NL  1.85 21.31   4.81  0.95 29200   80     7
## 20      PL  0.70 43.15   2.95  8.38  6400   20    21
## 21      PT 11.75 11.69   2.62  5.29 12400   34     6
## 22      RO  4.59 36.08   3.53  7.01  2900   20    16
## 23     SVK  6.03 40.62   4.76  5.64  6200   20    13
## 24     SLO  3.66 30.89   4.55  3.80 14400   19    12
## 25     SPA  6.03  8.00   2.48  6.33 17800   33     6
## 26     SWE  3.00 15.06   4.29  1.78 35300   80     6
## 27      UK  1.40  7.70   2.12 16.16 31400   80     5
summary(TABDATA)  ## we have meaningful values
##     country       polar             vol            nopart     
##  AT     : 1   Min.   : 0.700   Min.   : 2.80   Min.   :1.990  
##  BE     : 1   1st Qu.: 3.330   1st Qu.: 9.29   1st Qu.:2.970  
##  BG     : 1   Median : 5.910   Median :15.06   Median :3.540  
##  CYP    : 1   Mean   : 7.096   Mean   :21.85   Mean   :3.971  
##  CZ     : 1   3rd Qu.: 7.695   3rd Qu.:33.50   3rd Qu.:4.735  
##  DK     : 1   Max.   :33.470   Max.   :53.16   Max.   :9.050  
##  (Other):21                                                   
##       disp             GDP             demo           nogov      
##  Min.   : 0.660   Min.   : 2600   Min.   :17.00   Min.   : 4.00  
##  1st Qu.: 3.295   1st Qu.: 7350   1st Qu.:20.00   1st Qu.: 6.00  
##  Median : 5.210   Median :16500   Median :36.00   Median : 8.00  
##  Mean   : 5.700   Mean   :19733   Mean   :47.44   Mean   : 9.37  
##  3rd Qu.: 7.040   3rd Qu.:29050   3rd Qu.:80.00   3rd Qu.:12.00  
##  Max.   :18.370   Max.   :61500   Max.   :80.00   Max.   :21.00  
## 
is.numeric(TABDATA[ , 2])  ## as we would want it
## [1] TRUE
CSVDATA <- read.csv("example.csv")
CSVDATA <- read.csv("example.csv", sep = ";")  ## now simply everything is one variable in the data frame, because one element will be taken up until we find the separator
CSVDATA <- read.csv("example.csv", sep = ",")  ## If you know your data and system, this is redundant. Probably one of the few cases where I would recommend not to be parsimonious and use the separator (which works for the read.table() as well). We will be working with this data set, so let us call it something meaningful
EU <- CSVDATA  ## Note that this does not mean that the older data frame vanished. How do we check?

ls()  ## if we are a bit to crowded in memory, we drop some data frames (although they are trivially small, so for the sake of the exercize)
##  [1] "a"             "b"             "charvector"    "compList"     
##  [5] "CSVDATA"       "d"             "EU"            "exData"       
##  [9] "logicalvector" "newData"       "newList"       "numericVector"
## [13] "propMat"       "rNames"        "simpleList"    "SPSS"         
## [17] "STATA"         "STATA13"       "TABDATA"       "tranMat"      
## [21] "trialMat"      "x"             "x1"            "x2"           
## [25] "y"             "z"
rm(list = c("CSVDATA", "TABDATA"))
ls()  ## and we checked
##  [1] "a"             "b"             "charvector"    "compList"     
##  [5] "d"             "EU"            "exData"        "logicalvector"
##  [9] "newData"       "newList"       "numericVector" "propMat"      
## [13] "rNames"        "simpleList"    "SPSS"          "STATA"        
## [17] "STATA13"       "tranMat"       "trialMat"      "x"            
## [21] "x1"            "x2"            "y"             "z"

Working with Data

is.na(EU)  ## this will be too hard if larger dataset.
##       country polar   vol nopart  disp   GDP  demo nogov
##  [1,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [2,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [3,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [4,]   FALSE  TRUE  TRUE   TRUE  TRUE  TRUE  TRUE  TRUE
##  [5,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [6,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [7,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [8,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
##  [9,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [10,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [11,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [12,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [13,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [14,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [15,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [16,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [17,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [18,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [19,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [20,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [21,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [22,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [23,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [24,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [25,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [26,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [27,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
## [28,]   FALSE FALSE FALSE  FALSE FALSE FALSE FALSE FALSE
summary(EU)  ## we know that some variables have missing values
##     country       polar             vol            nopart     
##  AT     : 1   Min.   : 0.700   Min.   : 2.80   Min.   :1.990  
##  BE     : 1   1st Qu.: 3.330   1st Qu.: 9.29   1st Qu.:2.970  
##  BG     : 1   Median : 5.910   Median :15.06   Median :3.540  
##  BLR    : 1   Mean   : 7.096   Mean   :21.85   Mean   :3.971  
##  CYP    : 1   3rd Qu.: 7.695   3rd Qu.:33.50   3rd Qu.:4.735  
##  CZ     : 1   Max.   :33.470   Max.   :53.16   Max.   :9.050  
##  (Other):22   NA's   :1        NA's   :1       NA's   :1      
##       disp             GDP             demo           nogov      
##  Min.   : 0.660   Min.   : 2600   Min.   :17.00   Min.   : 4.00  
##  1st Qu.: 3.295   1st Qu.: 7350   1st Qu.:20.00   1st Qu.: 6.00  
##  Median : 5.210   Median :16500   Median :36.00   Median : 8.00  
##  Mean   : 5.700   Mean   :19733   Mean   :47.44   Mean   : 9.37  
##  3rd Qu.: 7.040   3rd Qu.:29050   3rd Qu.:80.00   3rd Qu.:12.00  
##  Max.   :18.370   Max.   :61500   Max.   :80.00   Max.   :21.00  
##  NA's   :1        NA's   :1       NA's   :1       NA's   :1
EU[is.na(EU$polar), ]  ## so, we want to see all the information for the observation where polarization is missing. Note: in R variable/observation == NA does not work. For system missing you always use is.na(). Let us drop the observation where we have a missing value. There some alternatives.
##   country polar vol nopart disp GDP demo nogov
## 4     BLR    NA  NA     NA   NA  NA   NA    NA
EUNM <- EU[is.na(EU$polar) == F, ]  ## keep all information if polarization has valid values
EUNM <- subset(EU, is.na(polar) == F)  ## use the subset function in R. It is a very useful function, where you stipulate first which data do you want to subset, and then a condition (note: it is enough to use the variable name, without the data frame name)
EUNM <- na.omit(EU)  ## Simplest way through the canned R function that will drop EVERYTHING that is missing. Note: everything here means that if you have an observation with a missing value for only one variable, the whole observation will be dropped! If you want to keep observations with some missing values, this is not the way to go.
mean(EU$polar)
## [1] NA
mean(EU$polar, na.rm = T) == mean(EUNM$polar)
## [1] TRUE

A few words on the attach function. DON’T USE IT! (Discuss why!)

Recoding

Now lets do some recoding.

EUNM$govDummy <- NA  ## a new variable in our dataset that only contains missing values
EUNM$govDummy[EUNM$nogov <= 8] <- 0  ## assign the value 0 to the new variable, but only if a condition holds true: nogov is smaller or equal to 8. If you do not pre-define the variable and already start with a condition, you will run into the problem of equal length. Why?
EUNM  ## one step done
##    country polar   vol nopart  disp   GDP demo nogov govDummy
## 1       AT  1.40 12.39   3.41  2.47 28900   56     8        0
## 2       BE  2.97 13.93   9.05  4.46 27200   80     9       NA
## 3       BG  6.43 42.09   2.52  7.37  2600   20     9       NA
## 5      CYP 33.47  6.89   3.52  2.82 16500   50     4        0
## 6       CZ  9.97 29.29   3.71  5.21  8100   20    12       NA
## 7       DK  3.74 10.95   4.71  1.50 35600   80     9       NA
## 8       EE  8.10 36.37   5.50  4.62  7900   18    12       NA
## 9      FIN  7.28  9.22   5.15  3.31 30900   80     6        0
## 10      FR  7.29 18.96   3.54 18.37 25700   80     8        0
## 11     GER  4.07  8.10   3.31  3.28 27300   62     5        0
## 12     GRE 11.66  9.36   2.21  7.07 16400   36     7        0
## 13     HUN  5.75 23.61   4.08  9.34  6600   20     9       NA
## 14     IRE  5.91 10.28   2.99  6.03 35100   80     8        0
## 15     ITA  6.59 30.92   6.07  5.62 21700   63    12       NA
## 16     LAT  8.57 53.16   5.49  3.93  6800   17    16       NA
## 17     LIT  2.50 48.72   3.52  8.34  6300   18    16       NA
## 18     LUX  3.87  8.43   4.34  4.15 61500   80     5        0
## 19      MT 23.00  2.80   1.99  0.66 11700   35     6        0
## 20      NL  1.85 21.31   4.81  0.95 29200   80     7        0
## 21      PL  0.70 43.15   2.95  8.38  6400   20    21       NA
## 22      PT 11.75 11.69   2.62  5.29 12400   34     6        0
## 23      RO  4.59 36.08   3.53  7.01  2900   20    16       NA
## 24     SVK  6.03 40.62   4.76  5.64  6200   20    13       NA
## 25     SLO  3.66 30.89   4.55  3.80 14400   19    12       NA
## 26     SPA  6.03  8.00   2.48  6.33 17800   33     6        0
## 27     SWE  3.00 15.06   4.29  1.78 35300   80     6        0
## 28      UK  1.40  7.70   2.12 16.16 31400   80     5        0
EUNM$govDummy[EUNM$nogov > 8] <- 1  ## assign the value 1 to the new variable, but only if a condition holds true: nogov is larger to 8. If you are absolutely sure that there are no missing values, you can take a shortcut: define the new variable with one of the values, and set up only one additional condition
EUNM$govDummyIF <- ifelse(EUNM$nogov <= 8, 0, 1)  ## If the new variable has only two values, ifelse() does the trick. You tell the condition, and then the first argument is the value for the new variable if condition is true, and the last one if the condition is false.


#An example of mathematical recoding. New variable is created by taking an average of two already existing variables - polarization and volatility.
EUNM$newVar <- (EUNM$polar + EUNM$vol) / 2

head(EUNM)
##   country polar   vol nopart disp   GDP demo nogov govDummy govDummyIF
## 1      AT  1.40 12.39   3.41 2.47 28900   56     8        0          0
## 2      BE  2.97 13.93   9.05 4.46 27200   80     9        1          1
## 3      BG  6.43 42.09   2.52 7.37  2600   20     9        1          1
## 5     CYP 33.47  6.89   3.52 2.82 16500   50     4        0          0
## 6      CZ  9.97 29.29   3.71 5.21  8100   20    12        1          1
## 7      DK  3.74 10.95   4.71 1.50 35600   80     9        1          1
##   newVar
## 1  6.895
## 2  8.450
## 3 24.260
## 5 20.180
## 6 19.630
## 7  7.345

With more categories and more conditions, logic is the same. We want a new variable that will be:

Now, we will take the shortcut here, because it is a good proxy for checking whether you understand the process:

Frequency and tables

EUNM$stability <- "Low"
EUNM$stability[(EUNM$polar > 5 & EUNM$polar <= 8) & (EUNM$nogov > 8 & EUNM$nogov <= 
    11)] <- "Medium"
EUNM$stability[(EUNM$polar <= 5) & (EUNM$nogov <= 8)] <- "High"  ## Let us discuss these lines
table(EUNM$stability, EUNM$polar)  ## simple cross-table
##         
##          0.7 1.4 1.85 2.5 2.97 3 3.66 3.74 3.87 4.07 4.59 5.75 5.91 6.03
##   High     0   2    1   0    0 1    0    0    1    1    0    0    0    0
##   Low      1   0    0   1    1 0    1    1    0    0    1    0    1    2
##   Medium   0   0    0   0    0 0    0    0    0    0    0    1    0    0
##         
##          6.43 6.59 7.28 7.29 8.1 8.57 9.97 11.66 11.75 23 33.47
##   High      0    0    0    0   0    0    0     0     0  0     0
##   Low       0    1    1    1   1    1    1     1     1  1     1
##   Medium    1    0    0    0   0    0    0     0     0  0     0
table(EUNM$stability, EUNM$nogov)  ## simple cross-table. Three-way tables are a bit harder to overview in this case. But there are possibilities to get at them. First, let us get percentages instead of counts
##         
##          4 5 6 7 8 9 12 13 16 21
##   High   0 3 1 1 1 0  0  0  0  0
##   Low    1 0 4 1 2 2  4  1  3  1
##   Medium 0 0 0 0 0 2  0  0  0  0
tTab <- table(EUNM$stability, EUNM$nogov)
prop.table(tTab, 1)  ## row percentages. You get column percentages by telling R you want the calculating on the second dimension of the table matrix (prop.table(tTab, 2))
##         
##                   4          5          6          7          8          9
##   High   0.00000000 0.50000000 0.16666667 0.16666667 0.16666667 0.00000000
##   Low    0.05263158 0.00000000 0.21052632 0.05263158 0.10526316 0.10526316
##   Medium 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 1.00000000
##         
##                  12         13         16         21
##   High   0.00000000 0.00000000 0.00000000 0.00000000
##   Low    0.21052632 0.05263158 0.15789474 0.05263158
##   Medium 0.00000000 0.00000000 0.00000000 0.00000000
## For three way tables let us create another dichotmous variable
EUNM$dispDummy <- 0
EUNM$dispDummy[EUNM$disp > 3.5] <- 1  ## Clear what we did?

tTabThree <- table(EUNM$stability, EUNM$govDummy, EUNM$dispDummy)
ftable(tTabThree)  ## You can see why it can get very complicated if not low number of categories. Also, as you saw: the table returned by the table function can be stored (it is an object). It is useful if you want to work with these frequencies. You can transform them simply in matrices or data frames
##            0  1
##                
## High   0   4  2
##        1   0  0
## Low    0   3  5
##        1   1 10
## Medium 0   0  0
##        1   0  2

Factors

class(EUNM$stability)  ## We have created this variable, and given the values we have assigned it is a character variable. But we know more about this, as there is some orderding, so this might not be good enough for us
## [1] "character"
table(EUNM$stability)
## 
##   High    Low Medium 
##      6     19      2
EUNM$stabFactor <- factor(EUNM$stability)  ## We should treat this as a factor variable, that has levels
levels(EUNM$stabFactor)  ## But our goal was to incorporate some ordering in this variable. What is the ordering right now? Actually, there is no real ordering.
## [1] "High"   "Low"    "Medium"
is.ordered(EUNM$stabFactor)
## [1] FALSE
EUNM$stabFactor <- factor(EUNM$stability, ordered = T)  ## Now, there will be ordering, but is it correct?
is.ordered(EUNM$stabFactor)
## [1] TRUE
levels(EUNM$stabFactor)  ## Well, if your goal is to have alphabetical ordering, sure. But we have qualitative knowledge about the order
## [1] "High"   "Low"    "Medium"
EUNM$stabFactor <- factor(EUNM$stability, levels = c("Low", "Medium", "High"), 
    ordered = T)
levels(EUNM$stabFactor)  ## This is more what we are looking for. It is important for modeling and plotting. Character vectors can be included in regression without a problem, as R does the factor transformation for you. But assumes no order, and you have no control on the baseline category.
## [1] "Low"    "Medium" "High"
## One more useful functions here:
EUNM$stabFactor <- factor(EUNM$stability)  ## Get the simplest version
EUNM$stabFactor <- relevel(EUNM$stabFactor, ref = "Medium")  ## If you do not necessarily want an ordered factor (happens many times), but you care about the reference category, this is the way to go
levels(EUNM$stabFactor)  ## see, not ordered, but you exerted control on the baseline. Factors are easy to understand and manipulate, just a matter of exercise
## [1] "Medium" "High"   "Low"

Saving Data

write.table(EU, file = "eu-data.txt", sep = " ", row.names = F)  ## it will save in your working directory. Object and file name are always mandatory, for other settings, such as default separator you can check the help file. Append is useful if you are iterating and saving out some parts separately. If you do not wish to save out variable names (maybe preparing something for the good old Mplus), set col.names = F
write.csv(EU, file = "eu-data.csv", row.names = F)  ## here, the separator is assumed, but arguments this function takes are the same

Factors

class(EUNM$stability)  ## We have created this variable, and given the values we have assigned it is a character variable. But we know more about this, as there is some orderding, so this might not be good enough for us
## [1] "character"
table(EUNM$stability)
## 
##   High    Low Medium 
##      6     19      2
EUNM$stabFactor <- factor(EUNM$stability)  ## We should treat this as a factor variable, that has levels
levels(EUNM$stabFactor)  ## But our goal was to incorporate some ordering in this variable. What is the ordering right now? Actually, there is no real ordering.
## [1] "High"   "Low"    "Medium"
EUNM$stabFactor <- factor(EUNM$stability, ordered = T)  ## Now, there will be ordering, but is it correct?
is.ordered(EUNM$stabFactor)
## [1] TRUE
levels(EUNM$stabFactor)  ## Well, if your goal is to have alphabetical ordering, sure. But we have qualitative knowledge about the order
## [1] "High"   "Low"    "Medium"
EUNM$stabFactor <- factor(EUNM$stability, levels = c("Low", "Medium", "High"), ordered = T)
                          
levels(EUNM$stabFactor)  ## This is more what we are looking for. It is important for modeling and plotting. Character vectors can be included in regression without a problem, as R does the factor transformation for you. But assumes no order, and you have no control on the baseline category.
## [1] "Low"    "Medium" "High"

One more useful functions here:

EUNM$stabFactor <- factor(EUNM$stability)  ## Get the simplest version
EUNM$stabFactor <- relevel(EUNM$stabFactor, ref = "Medium")  ## If you do not necessarily want an ordered factor (happens many times), but you care about the reference category, this is the way to go
levels(EUNM$stabFactor)  ## see, not ordered, but you exerted control on the baseline. Factors are easy to understand and manipulate, just a matter of exercise
## [1] "Medium" "High"   "Low"
# OK.  So lets hget rid of these text variables here.  (Did someone ask how to drop variables?)
EUNM
##    country polar   vol nopart  disp   GDP demo nogov govDummy govDummyIF
## 1       AT  1.40 12.39   3.41  2.47 28900   56     8        0          0
## 2       BE  2.97 13.93   9.05  4.46 27200   80     9        1          1
## 3       BG  6.43 42.09   2.52  7.37  2600   20     9        1          1
## 5      CYP 33.47  6.89   3.52  2.82 16500   50     4        0          0
## 6       CZ  9.97 29.29   3.71  5.21  8100   20    12        1          1
## 7       DK  3.74 10.95   4.71  1.50 35600   80     9        1          1
## 8       EE  8.10 36.37   5.50  4.62  7900   18    12        1          1
## 9      FIN  7.28  9.22   5.15  3.31 30900   80     6        0          0
## 10      FR  7.29 18.96   3.54 18.37 25700   80     8        0          0
## 11     GER  4.07  8.10   3.31  3.28 27300   62     5        0          0
## 12     GRE 11.66  9.36   2.21  7.07 16400   36     7        0          0
## 13     HUN  5.75 23.61   4.08  9.34  6600   20     9        1          1
## 14     IRE  5.91 10.28   2.99  6.03 35100   80     8        0          0
## 15     ITA  6.59 30.92   6.07  5.62 21700   63    12        1          1
## 16     LAT  8.57 53.16   5.49  3.93  6800   17    16        1          1
## 17     LIT  2.50 48.72   3.52  8.34  6300   18    16        1          1
## 18     LUX  3.87  8.43   4.34  4.15 61500   80     5        0          0
## 19      MT 23.00  2.80   1.99  0.66 11700   35     6        0          0
## 20      NL  1.85 21.31   4.81  0.95 29200   80     7        0          0
## 21      PL  0.70 43.15   2.95  8.38  6400   20    21        1          1
## 22      PT 11.75 11.69   2.62  5.29 12400   34     6        0          0
## 23      RO  4.59 36.08   3.53  7.01  2900   20    16        1          1
## 24     SVK  6.03 40.62   4.76  5.64  6200   20    13        1          1
## 25     SLO  3.66 30.89   4.55  3.80 14400   19    12        1          1
## 26     SPA  6.03  8.00   2.48  6.33 17800   33     6        0          0
## 27     SWE  3.00 15.06   4.29  1.78 35300   80     6        0          0
## 28      UK  1.40  7.70   2.12 16.16 31400   80     5        0          0
##    newVar stability dispDummy stabFactor
## 1   6.895      High         0       High
## 2   8.450       Low         1        Low
## 3  24.260    Medium         1     Medium
## 5  20.180       Low         0        Low
## 6  19.630       Low         1        Low
## 7   7.345       Low         0        Low
## 8  22.235       Low         1        Low
## 9   8.250       Low         0        Low
## 10 13.125       Low         1        Low
## 11  6.085      High         0       High
## 12 10.510       Low         1        Low
## 13 14.680    Medium         1     Medium
## 14  8.095       Low         1        Low
## 15 18.755       Low         1        Low
## 16 30.865       Low         1        Low
## 17 25.610       Low         1        Low
## 18  6.150      High         1       High
## 19 12.900       Low         0        Low
## 20 11.580      High         0       High
## 21 21.925       Low         1        Low
## 22 11.720       Low         1        Low
## 23 20.335       Low         1        Low
## 24 23.325       Low         1        Low
## 25 17.275       Low         1        Low
## 26  7.015       Low         1        Low
## 27  9.030      High         0       High
## 28  4.550      High         1       High
EUNM <- EUNM[, 1:8]
EUNM
##    country polar   vol nopart  disp   GDP demo nogov
## 1       AT  1.40 12.39   3.41  2.47 28900   56     8
## 2       BE  2.97 13.93   9.05  4.46 27200   80     9
## 3       BG  6.43 42.09   2.52  7.37  2600   20     9
## 5      CYP 33.47  6.89   3.52  2.82 16500   50     4
## 6       CZ  9.97 29.29   3.71  5.21  8100   20    12
## 7       DK  3.74 10.95   4.71  1.50 35600   80     9
## 8       EE  8.10 36.37   5.50  4.62  7900   18    12
## 9      FIN  7.28  9.22   5.15  3.31 30900   80     6
## 10      FR  7.29 18.96   3.54 18.37 25700   80     8
## 11     GER  4.07  8.10   3.31  3.28 27300   62     5
## 12     GRE 11.66  9.36   2.21  7.07 16400   36     7
## 13     HUN  5.75 23.61   4.08  9.34  6600   20     9
## 14     IRE  5.91 10.28   2.99  6.03 35100   80     8
## 15     ITA  6.59 30.92   6.07  5.62 21700   63    12
## 16     LAT  8.57 53.16   5.49  3.93  6800   17    16
## 17     LIT  2.50 48.72   3.52  8.34  6300   18    16
## 18     LUX  3.87  8.43   4.34  4.15 61500   80     5
## 19      MT 23.00  2.80   1.99  0.66 11700   35     6
## 20      NL  1.85 21.31   4.81  0.95 29200   80     7
## 21      PL  0.70 43.15   2.95  8.38  6400   20    21
## 22      PT 11.75 11.69   2.62  5.29 12400   34     6
## 23      RO  4.59 36.08   3.53  7.01  2900   20    16
## 24     SVK  6.03 40.62   4.76  5.64  6200   20    13
## 25     SLO  3.66 30.89   4.55  3.80 14400   19    12
## 26     SPA  6.03  8.00   2.48  6.33 17800   33     6
## 27     SWE  3.00 15.06   4.29  1.78 35300   80     6
## 28      UK  1.40  7.70   2.12 16.16 31400   80     5

Functions

Writing your own functions - If you have to do something more than twice, write a function for it!) Here’s a simple example. (Embrace it or suffer…)

standardize <- function(x) {
  ## {} marks code chunk that is executed together. If you are doing multiple operations, you can use return(object of interest) at the end of the function. In the () we specify the arguments the function takes, in this case a vector that we name x
  (x - mean(x, na.rm = T))/sd(x, na.rm = T)
}  ## You select the code chunk of the function, and run it. By this, it is defined and R will know it as an object

hist(EUNM$nopart)

# So lets standardize it
EUNM$nopartStand <- standardize(EUNM$nopart)  ## we have created a new variable in our data that contains the standardized age
hist(EUNM$nopartStand)

x - Yes, this is about environments: when looking in our global environment where we work, x is not defined. x is only defined within the local environment of the function. You can use various names

Now I believe you are armed with Basic data management skills to starts doing stats

Just think:

Are we ready for stats? :D

Let me show you a fun trick now that you know packages. (But we need a package for it)

#install.packages("psych", dependencies = T)
library(psych)

(EUNMDesc <- psych::describe(EUNM))  ## Another great tool to check your data. Asterix marks non-numeric variables, so of course some statistics offered here are misleading.
##             vars  n     mean       sd   median  trimmed      mad     min
## country*       1 27    14.89     8.12    15.00    14.96    10.38    1.00
## polar          2 27     7.10     6.92     5.91     5.78     3.34    0.70
## vol            3 27    21.85    14.91    15.06    20.80    10.91    2.80
## nopart         4 27     3.97     1.51     3.54     3.83     1.50    1.99
## disp           5 27     5.70     4.06     5.21     5.12     2.82    0.66
## GDP            6 27 19733.33 13821.75 16500.00 18704.35 15122.52 2600.00
## demo           7 27    47.44    26.84    36.00    47.22    26.69   17.00
## nogov          8 27     9.37     4.24     8.00     9.00     2.97    4.00
## nopartStand    9 27     0.00     1.00    -0.29    -0.10     0.99   -1.32
##                  max    range  skew kurtosis      se
## country*       28.00    27.00 -0.06    -1.28    1.56
## polar          33.47    32.77  2.39     5.95    1.33
## vol            53.16    50.36  0.58    -1.11    2.87
## nopart          9.05     7.06  1.30     2.42    0.29
## disp           18.37    17.71  1.55     2.47    0.78
## GDP         61500.00 58900.00  0.90     0.70 2660.00
## demo           80.00    63.00  0.16    -1.80    5.17
## nogov          21.00    17.00  0.94     0.09    0.82
## nopartStand     3.37     4.69  1.30     2.42    0.19

Notice that normally we would write describe(EUNM)? Discuss Masking

Lets look at what is in that object we created. It is ALL THERE!

str(EUNMDesc)
## Classes 'psych', 'describe' and 'data.frame':    9 obs. of  13 variables:
##  $ vars    : int  1 2 3 4 5 6 7 8 9
##  $ n       : num  27 27 27 27 27 27 27 27 27
##  $ mean    : num  14.89 7.1 21.85 3.97 5.7 ...
##  $ sd      : num  8.12 6.92 14.91 1.51 4.06 ...
##  $ median  : num  15 5.91 15.06 3.54 5.21 ...
##  $ trimmed : num  14.96 5.78 20.8 3.83 5.12 ...
##  $ mad     : num  10.38 3.34 10.91 1.5 2.82 ...
##  $ min     : num  1 0.7 2.8 1.99 0.66 ...
##  $ max     : num  28 33.47 53.16 9.05 18.37 ...
##  $ range   : num  27 32.77 50.36 7.06 17.71 ...
##  $ skew    : num  -0.0583 2.3917 0.5832 1.2974 1.554 ...
##  $ kurtosis: num  -1.28 5.95 -1.11 2.42 2.47 ...
##  $ se      : num  1.562 1.332 2.869 0.29 0.781 ...

Want the standard errors? Sure

EUNMDesc$se
## [1]    1.5619527    1.3321965    2.8694346    0.2896896    0.7806087
## [6] 2659.9970010    5.1662761    0.8167550    0.1924501

Want the standard error of effective number of parties? Sure

EUNMDesc$se[4]
## [1] 0.2896896

Or you can write your own function :) (Can even make sure the function watches for missing)

se <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x))) #function for standard error

se(EUNM$nopart)
## [1] 0.2896896

want confidence intervals?

mean(EUNM$nopart, na.omit = T) # Check mean (I added to drop missing but we don't have any in this case)
## [1] 3.971111

Lower CI 95%

mean(EUNM$nopart, na.omit = T) - se(EUNM$nopart) * qt(0.025, length(na.omit(EUNM$nopart)) - 1, lower.tail = F)
## [1] 3.375646

Upper CI 95%

mean(EUNM$nopart, na.omit = T) + se(EUNM$nopart) * qt(0.025, length(na.omit(EUNM$nopart)) - 1, lower.tail = F)
## [1] 4.566577

Notice the qt command that gets us the critical t value for 0.025 probability in the tail. Degrees of freedom is the length of the thing that we calculate the CI for (I added to omit missing).# Lets not forget -1 because degrees of freedom is n-1 for this calculation.

Alright. We have chapters 1-8 and 12 covered. Now, for the t-test.

Now for Chapters 9-11 :)
Lets look at the t.test command’s help.

?t.test

Great. You can set confidence level with conf.level and you can set if it is a paired sample or not. Can do two tailed (default) or one tailed (never!!!). What’s this var.equal thing? Lets talk about the Welch test!

One Sample

Effective Number of Parties different than 3? (Can we say we have multiparty system in Europe?)

t.test(EUNM$nopart, mu=3)
## 
##  One Sample t-test
## 
## data:  EUNM$nopart
## t = 3.3522, df = 26, p-value = 0.002464
## alternative hypothesis: true mean is not equal to 3
## 95 percent confidence interval:
##  3.375646 4.566577
## sample estimates:
## mean of x 
##  3.971111

Two Sample

Let’s check if the wealthier countries differ in terms of political polarization from the poorer ones

EUNM$wealth <- 1
EUNM$wealth[EUNM$GDP < median(EUNM$GDP)] <- 0
(wealth_split <- split(EUNM$polar, EUNM$wealth))
## $`0`
##  [1]  6.43  9.97  8.10 11.66  5.75  8.57  2.50 23.00  0.70 11.75  4.59
## [12]  6.03  3.66
## 
## $`1`
##  [1]  1.40  2.97 33.47  3.74  7.28  7.29  4.07  5.91  6.59  3.87  1.85
## [12]  6.03  3.00  1.40
x <- as.vector(wealth_split[[1]])
y <- as.vector(wealth_split[[2]])
t.test(x,y)
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = 0.58255, df = 23.289, p-value = 0.5658
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.957715  7.063539
## sample estimates:
## mean of x mean of y 
##  7.900769  6.347857

Paired Sample

There is nothing in this data to do this with. But. How would you do it? Two ways.:

  1. Just feed it the vector of differences and set mu=0 (Actually it is 0 by default)
  2. Feed it two vectors (of same length AND ORDER) like it two sample t test and say paired = TRUE.

I prefer the first

Chi-square test

The package you want is called descr

#install.packages("descr", dependencies = T)
library(descr)

#Now let's restore the stability variable

EUNM$stability <- "Low"
EUNM$stability[(EUNM$polar > 5 & EUNM$polar <= 8) & (EUNM$nogov > 8 & EUNM$nogov <= 11)] <- "Medium"
EUNM$stability[(EUNM$polar <= 5) & (EUNM$nogov <= 8)] <- "High"  

EUNM$stabFactor <- factor(EUNM$stability, levels = c("Low", "Medium", "High"), ordered = T)

crosstab(EUNM$wealth, EUNM$stabFactor, chisq = T, plot = F) #by default, the crosstab function wants to create a pretty pointless chart
## Warning in chisq.test(tab, correct = FALSE, ...): Chi-squared approximation
## may be incorrect
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |-------------------------|
## 
## ==========================================
##                EUNM$stabFactor
## EUNM$wealth    Low   Medium   High   Total
## ------------------------------------------
## 0               11        2      0      13
## ------------------------------------------
## 1                8        0      6      14
## ------------------------------------------
## Total           19        2      6      27
## ==========================================
## 
## Statistics for All Table Factors
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 = 8.448236      d.f. = 2      p = 0.0146 
## 
##         Minimum expected frequency: 0.962963 
## Cells with Expected Frequency < 5: 4 of 6 (66.66667%)
anest04 <- read.csv("a02anes2004trimmed.csv")

So - quick review with a different dataset

summary(anest04)  # We see a few missing values.  Nothing substantial.  Can probably drop them.
##        id              pid           presvote          edu      
##  Min.   :  26.0   Min.   :0.000   Min.   :1.000   Min.   :1.00  
##  1st Qu.: 320.0   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:3.00  
##  Median : 624.0   Median :2.000   Median :1.000   Median :4.00  
##  Mean   : 620.0   Mean   :2.781   Mean   :1.477   Mean   :4.12  
##  3rd Qu.: 901.5   3rd Qu.:5.000   3rd Qu.:2.000   3rd Qu.:6.00  
##  Max.   :1206.0   Max.   :6.000   Max.   :2.000   Max.   :7.00  
##                   NA's   :5       NA's   :16                    
##     womnotru    
##  Min.   :1.000  
##  1st Qu.:3.000  
##  Median :4.000  
##  Mean   :3.794  
##  3rd Qu.:5.000  
##  Max.   :5.000  
##  NA's   :17

In cases when you have many variables, you want to be very careful with this.

anest04nomiss <- na.omit(anest04)

#install.packages("psych", dependencies = T)
library(psych)
describe(anest04nomiss)
##          vars   n   mean     sd median trimmed    mad min  max range  skew
## id          1 157 606.22 339.90    597  603.76 437.37  26 1206  1180  0.05
## pid         2 157   2.85   2.11      2    2.81   2.97   0    6     6  0.10
## presvote    3 157   1.48   0.50      1    1.48   0.00   1    2     1  0.06
## edu         4 157   4.24   1.61      4    4.20   1.48   1    7     6  0.26
## womnotru    5 157   3.87   1.13      4    3.99   1.48   1    5     4 -0.73
##          kurtosis    se
## id          -1.19 27.13
## pid         -1.42  0.17
## presvote    -2.01  0.04
## edu         -1.02  0.13
## womnotru    -0.51  0.09

So lets review the t test. Can we say that the US is not sexist?
(Middle value is 3 - we can cosnider anything below that sexist)

t.test(anest04nomiss$womnotru, mu=3)
## 
##  One Sample t-test
## 
## data:  anest04nomiss$womnotru
## t = 9.6697, df = 156, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 3
## 95 percent confidence interval:
##  3.694358 4.050865
## sample estimates:
## mean of x 
##  3.872611

Lets see if there is significant difference between Bush voters and Kerry voters? (I am going to do it differently from yesterday on purpose.)

First lets grab all the Kerry Voters (creates same dataframe with all the Kerry voters)

x <- subset(anest04nomiss, presvote == 1)

Then all the Bush voters (creates same dataframe with all the Bush voters)

y <- subset(anest04nomiss, presvote == 2)

t.test(x$womnotru,y$womnotru)
## 
##  Welch Two Sample t-test
## 
## data:  x$womnotru and y$womnotru
## t = 2.1816, df = 150.2, p-value = 0.03069
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.03683944 0.74451209
## sample estimates:
## mean of x mean of y 
##  4.061728  3.671053

Or if you prefer, we can do it this way too :)

t.test(anest04nomiss$womnotru ~ anest04nomiss$presvote)
## 
##  Welch Two Sample t-test
## 
## data:  anest04nomiss$womnotru by anest04nomiss$presvote
## t = 2.1816, df = 150.2, p-value = 0.03069
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.03683944 0.74451209
## sample estimates:
## mean in group 1 mean in group 2 
##        4.061728        3.671053

OK. So paired sample t test Lets load a different dataset. (Should look familiar.)
Note that other dataframe stays in memory. We will get back to it later.

imag <- read.csv("imaginationtherapy.csv")
summary(imag)
##  A          X15             X24      
##  B:1   Min.   : 5.00   Min.   :11.0  
##  C:1   1st Qu.: 7.50   1st Qu.:14.0  
##  D:1   Median : 9.50   Median :17.5  
##  E:1   Mean   :10.17   Mean   :18.0  
##  F:1   3rd Qu.:11.50   3rd Qu.:22.5  
##  G:1   Max.   :18.00   Max.   :25.0

Something went wrong.

imag <- read.csv("imaginationtherapy.csv", header = F)
summary(imag)
##  V1          V2              V3       
##  A:1   Min.   : 5.00   Min.   :11.00  
##  B:1   1st Qu.: 8.00   1st Qu.:14.00  
##  C:1   Median :10.00   Median :21.00  
##  D:1   Mean   :10.86   Mean   :18.86  
##  E:1   3rd Qu.:13.50   3rd Qu.:23.50  
##  F:1   Max.   :18.00   Max.   :25.00  
##  G:1

Looks much better

t.test(imag$V3 - imag$V2)
## 
##  One Sample t-test
## 
## data:  imag$V3 - imag$V2
## t = 7.0553, df = 6, p-value = 0.0004058
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   5.225463 10.774537
## sample estimates:
## mean of x 
##         8

Or. t.test command has it all built in

t.test(imag$V3, imag$V2, paired = T)
## 
##  Paired t-test
## 
## data:  imag$V3 and imag$V2
## t = 7.0553, df = 6, p-value = 0.0004058
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   5.225463 10.774537
## sample estimates:
## mean of the differences 
##                       8

Crosstabulation

(some of you wanted prettier tables - this may help)

install.packages("descr", dependencies = T)
library(descr)
crosstab(anest04nomiss$edu, anest04nomiss$presvote)

##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |-------------------------|
## 
## ====================================
##                      anest04nomiss$presvote
## anest04nomiss$edu     1    2   Total
## ------------------------------------
## 1                     1    3       4
## ------------------------------------
## 2                     8    2      10
## ------------------------------------
## 3                    29   26      55
## ------------------------------------
## 4                    13   13      26
## ------------------------------------
## 5                    10    7      17
## ------------------------------------
## 6                    13   15      28
## ------------------------------------
## 7                     7   10      17
## ------------------------------------
## Total                81   76     157
## ====================================

OK. So default may not be the most useful. Lets see our options

?crosstab

Lets try it this way.

crosstab(anest04nomiss$edu, anest04nomiss$presvote, chisq = T, plot = F)
## Warning in chisq.test(tab, correct = FALSE, ...): Chi-squared approximation
## may be incorrect
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |-------------------------|
## 
## ====================================
##                      anest04nomiss$presvote
## anest04nomiss$edu     1    2   Total
## ------------------------------------
## 1                     1    3       4
## ------------------------------------
## 2                     8    2      10
## ------------------------------------
## 3                    29   26      55
## ------------------------------------
## 4                    13   13      26
## ------------------------------------
## 5                    10    7      17
## ------------------------------------
## 6                    13   15      28
## ------------------------------------
## 7                     7   10      17
## ------------------------------------
## Total                81   76     157
## ====================================
## 
## Statistics for All Table Factors
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 = 5.811976      d.f. = 6      p = 0.445 
## 
##         Minimum expected frequency: 1.936306 
## Cells with Expected Frequency < 5: 3 of 14 (21.42857%)

We just did a chi-square tes of independence on this cross-tab. Note the Warning! Why do we get this?

It is possible to do chi-square test without the cross-tab too BTW

chisq.test(anest04nomiss$edu, anest04nomiss$presvote)
## Warning in chisq.test(anest04nomiss$edu, anest04nomiss$presvote): Chi-
## squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  anest04nomiss$edu and anest04nomiss$presvote
## X-squared = 5.812, df = 6, p-value = 0.4446

We can even di a chi-square test of indifference. Lets just use the coke example In class we had 3 Pepsi supperters and 17 Coke supporters.

cola <- c(3,17)
chisq.test(cola)
## 
##  Chi-squared test for given probabilities
## 
## data:  cola
## X-squared = 9.8, df = 1, p-value = 0.001745

Feel free to explore the options you have with ?chisq.test. But now lets get back to our cross-tabs We may also want more than frequencies. Lets discuss these options.

crosstab(anest04nomiss$edu, anest04nomiss$presvote, plot = F, prop.t = T)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |           Total Percent | 
## |-------------------------|
## 
## ==========================================
##                      anest04nomiss$presvote
## anest04nomiss$edu        1       2   Total
## ------------------------------------------
## 1                       1       3       4 
##                       0.6%    1.9%        
## ------------------------------------------
## 2                       8       2      10 
##                       5.1%    1.3%        
## ------------------------------------------
## 3                      29      26      55 
##                      18.5%   16.6%        
## ------------------------------------------
## 4                      13      13      26 
##                       8.3%    8.3%        
## ------------------------------------------
## 5                      10       7      17 
##                       6.4%    4.5%        
## ------------------------------------------
## 6                      13      15      28 
##                       8.3%    9.6%        
## ------------------------------------------
## 7                       7      10      17 
##                       4.5%    6.4%        
## ------------------------------------------
## Total                  81      76     157 
## ==========================================
crosstab(anest04nomiss$edu, anest04nomiss$presvote, plot = F, prop.c = T)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |          Column Percent | 
## |-------------------------|
## 
## ==========================================
##                      anest04nomiss$presvote
## anest04nomiss$edu        1       2   Total
## ------------------------------------------
## 1                       1       3       4 
##                       1.2%    3.9%        
## ------------------------------------------
## 2                       8       2      10 
##                       9.9%    2.6%        
## ------------------------------------------
## 3                      29      26      55 
##                      35.8%   34.2%        
## ------------------------------------------
## 4                      13      13      26 
##                      16.0%   17.1%        
## ------------------------------------------
## 5                      10       7      17 
##                      12.3%    9.2%        
## ------------------------------------------
## 6                      13      15      28 
##                      16.0%   19.7%        
## ------------------------------------------
## 7                       7      10      17 
##                       8.6%   13.2%        
## ------------------------------------------
## Total                  81      76     157 
##                      51.6%   48.4%        
## ==========================================
crosstab(anest04nomiss$edu, anest04nomiss$presvote, plot = F, prop.r = T)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |             Row Percent | 
## |-------------------------|
## 
## ==========================================
##                      anest04nomiss$presvote
## anest04nomiss$edu        1       2   Total
## ------------------------------------------
## 1                       1       3       4 
##                      25.0%   75.0%    2.5%
## ------------------------------------------
## 2                       8       2      10 
##                      80.0%   20.0%    6.4%
## ------------------------------------------
## 3                      29      26      55 
##                      52.7%   47.3%   35.0%
## ------------------------------------------
## 4                      13      13      26 
##                      50.0%   50.0%   16.6%
## ------------------------------------------
## 5                      10       7      17 
##                      58.8%   41.2%   10.8%
## ------------------------------------------
## 6                      13      15      28 
##                      46.4%   53.6%   17.8%
## ------------------------------------------
## 7                       7      10      17 
##                      41.2%   58.8%   10.8%
## ------------------------------------------
## Total                  81      76     157 
## ==========================================

Correlation

?cor
summary(anest04nomiss)
##        id              pid           presvote          edu       
##  Min.   :  26.0   Min.   :0.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.: 298.0   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:3.000  
##  Median : 597.0   Median :2.000   Median :1.000   Median :4.000  
##  Mean   : 606.2   Mean   :2.847   Mean   :1.484   Mean   :4.236  
##  3rd Qu.: 884.0   3rd Qu.:5.000   3rd Qu.:2.000   3rd Qu.:6.000  
##  Max.   :1206.0   Max.   :6.000   Max.   :2.000   Max.   :7.000  
##     womnotru    
##  Min.   :1.000  
##  1st Qu.:3.000  
##  Median :4.000  
##  Mean   :3.873  
##  3rd Qu.:5.000  
##  Max.   :5.000

Note also what I am doing with the data matrix!

cor(anest04nomiss[c(2,4:5)])
##                 pid       edu   womnotru
## pid       1.0000000 0.1878255 -0.1316365
## edu       0.1878255 1.0000000  0.1398544
## womnotru -0.1316365 0.1398544  1.0000000

Correlations with significance levels

install.packages("Hmisc", dependencies = T)
# install.packages("Hmisc", dependencies = T)
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## 
## Attaching package: 'Hmisc'
## The following object is masked from 'package:psych':
## 
##     describe
## The following objects are masked from 'package:base':
## 
##     format.pval, units
rcorr(as.matrix(anest04nomiss[c(2,4:5)]), type="pearson") 
##            pid  edu womnotru
## pid       1.00 0.19    -0.13
## edu       0.19 1.00     0.14
## womnotru -0.13 0.14     1.00
## 
## n= 157 
## 
## 
## P
##          pid    edu    womnotru
## pid             0.0185 0.1003  
## edu      0.0185        0.0806  
## womnotru 0.1003 0.0806

We can also do a spearman correlation here as well (may be more appropriate anyway)

rcorr(as.matrix(anest04nomiss[c(2,4:5)]), type="spearman") # type can be pearson or spearman
##            pid  edu womnotru
## pid       1.00 0.20    -0.12
## edu       0.20 1.00     0.12
## womnotru -0.12 0.12     1.00
## 
## n= 157 
## 
## 
## P
##          pid    edu    womnotru
## pid             0.0142 0.1494  
## edu      0.0142        0.1388  
## womnotru 0.1494 0.1388

For some reason it doesn’t like data frames. Don’t ask why.

Regressions

Finally. Running regressions. Lets start with a bivariate case. Use sexism as the dependent variable. Use party ID as IV. What do you expect to find?

lm(womnotru ~ pid, data = anest04nomiss)
## 
## Call:
## lm(formula = womnotru ~ pid, data = anest04nomiss)
## 
## Coefficients:
## (Intercept)          pid  
##     4.07320     -0.07045

That didn’t tell us much But appearances can be misleading

bivreg <- lm(womnotru ~ pid, data = anest04nomiss)

Now lets take a better look at this object

str(bivreg)
## List of 12
##  $ coefficients : Named num [1:2] 4.0732 -0.0705
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "pid"
##  $ residuals    : Named num [1:157] 0.997 1.209 -0.932 0.35 0.279 ...
##   ..- attr(*, "names")= chr [1:157] "1" "2" "4" "5" ...
##  $ effects      : Named num [1:157] -48.524 1.859 -0.977 0.141 0.111 ...
##   ..- attr(*, "names")= chr [1:157] "(Intercept)" "pid" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:157] 4 3.79 3.93 3.65 3.72 ...
##   ..- attr(*, "names")= chr [1:157] "1" "2" "4" "5" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:157, 1:2] -12.53 0.0798 0.0798 0.0798 0.0798 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:157] "1" "2" "4" "5" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "pid"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.08 1.05
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 155
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = womnotru ~ pid, data = anest04nomiss)
##  $ terms        :Classes 'terms', 'formula'  language womnotru ~ pid
##   .. ..- attr(*, "variables")= language list(womnotru, pid)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "womnotru" "pid"
##   .. .. .. ..$ : chr "pid"
##   .. ..- attr(*, "term.labels")= chr "pid"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(womnotru, pid)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "womnotru" "pid"
##  $ model        :'data.frame':   157 obs. of  2 variables:
##   ..$ womnotru: int [1:157] 5 5 3 4 4 4 4 5 5 5 ...
##   ..$ pid     : int [1:157] 1 4 2 6 5 6 0 1 2 0 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language womnotru ~ pid
##   .. .. ..- attr(*, "variables")= language list(womnotru, pid)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "womnotru" "pid"
##   .. .. .. .. ..$ : chr "pid"
##   .. .. ..- attr(*, "term.labels")= chr "pid"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(womnotru, pid)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "womnotru" "pid"
##  - attr(*, "class")= chr "lm"

In fact

summary(bivreg)
## 
## Call:
## lm(formula = womnotru ~ pid, data = anest04nomiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8618 -0.7209  0.1382  1.0677  1.3495 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.07320    0.15091  26.991   <2e-16 ***
## pid         -0.07045    0.04261  -1.653      0.1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.124 on 155 degrees of freedom
## Multiple R-squared:  0.01733,    Adjusted R-squared:  0.01099 
## F-statistic: 2.733 on 1 and 155 DF,  p-value: 0.1003

In fact, multivariate extension is pretty easy

multreg <- lm(womnotru ~ pid + edu + as.factor(presvote), data = anest04nomiss)
summary(multreg)
## 
## Call:
## lm(formula = womnotru ~ pid + edu + as.factor(presvote), data = anest04nomiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8257 -0.6019  0.1172  1.0231  1.6748 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.62571    0.26195  13.841   <2e-16 ***
## pid                  -0.01964    0.07112  -0.276   0.7828    
## edu                   0.11190    0.05661   1.977   0.0499 *  
## as.factor(presvote)2 -0.35354    0.29532  -1.197   0.2331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.11 on 153 degrees of freedom
## Multiple R-squared:  0.05426,    Adjusted R-squared:  0.03572 
## F-statistic: 2.926 on 3 and 153 DF,  p-value: 0.03569

Why as.factor? We could just use the dummy as well.
But this is not dummy So lets create a dummy

anest04nomiss$bush <- anest04nomiss$presvote - 1

Lets rerun the model. We get the same result

multreg2 <- lm(womnotru ~ pid + edu + bush, data = anest04nomiss)
summary(multreg2)
## 
## Call:
## lm(formula = womnotru ~ pid + edu + bush, data = anest04nomiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8257 -0.6019  0.1172  1.0231  1.6748 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.62571    0.26195  13.841   <2e-16 ***
## pid         -0.01964    0.07112  -0.276   0.7828    
## edu          0.11190    0.05661   1.977   0.0499 *  
## bush        -0.35354    0.29532  -1.197   0.2331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.11 on 153 degrees of freedom
## Multiple R-squared:  0.05426,    Adjusted R-squared:  0.03572 
## F-statistic: 2.926 on 3 and 153 DF,  p-value: 0.03569

But if we just do it this way, what would go wrong

multreg3 <- lm(womnotru ~ pid + edu + presvote, data = anest04nomiss)
summary(multreg3)
## 
## Call:
## lm(formula = womnotru ~ pid + edu + presvote, data = anest04nomiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8257 -0.6019  0.1172  1.0231  1.6748 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.97925    0.40172   9.906   <2e-16 ***
## pid         -0.01964    0.07112  -0.276   0.7828    
## edu          0.11190    0.05661   1.977   0.0499 *  
## presvote    -0.35354    0.29532  -1.197   0.2331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.11 on 153 degrees of freedom
## Multiple R-squared:  0.05426,    Adjusted R-squared:  0.03572 
## F-statistic: 2.926 on 3 and 153 DF,  p-value: 0.03569

So, we are done. I hope it wasn’t as painful as someone assigning homework for the weekend (because that is what I am going to do now).

PRACTICE!!

One good way to practice is using the package: swirl

#install.packages("swirl", dependencies = T)
library(swirl)
## 
## | Hi! I see that you have some variables saved in your workspace. To keep
## | things running smoothly, I recommend you clean up before starting swirl.
## 
## | Type ls() to see a list of the variables in your workspace. Then, type
## | rm(list=ls()) to clear your workspace.
## 
## | Type swirl() when you are ready to begin.
#To run swirl uncomment below and run this command
#swirl()