Many of you have probably seen Fight Club and this course or working in R might be similar:
Some things make perfect sense in the end
Some things were very elegant, when knowing the story in the end
It does not hurt if you re-watch it once or twice
lots of things under the hood (first rule of FC)
It might get you satisfaction (or even a good job), but you might take a punch or two on the road
IDE (Integrated Development Environments). There are others: emacs and vi and sublime text (even steeper learning curve)
Scripting:
Never code in the console! Seriously don’t
Running code ctrl+r on Windows, cmd+enter on mac
In console, up and down arrow are useful
saving (not the console, the script) Use .R at the end.
line breaks (+s) in the console (not in the script)
probably too early (but it is never too early) See here for how to script well clicky1, clicky2
We insert comments in R using #
So:
#writewhateverhere
Does nothing
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
26 == (3 * 13) - 13
## [1] TRUE
26 == (3 * 13) - 12
## [1] FALSE
Get using getw:
getwd()
## [1] "/Users/levi/data"
Set using setwd()
setwd('~/data/') # on Windows, it is 'c:/data/'
getwd()
## [1] "/Users/levi/data"
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)
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
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
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
a <- c(1:25) # Lets make some data
sd(a) # Lets check the standard deviation
## [1] 7.359801
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
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
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"
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
(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)
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 ...
setwd("~/data/") #or whatever you want to use
2.Put the data from the elearning in that folder.** (All of them.)
## 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)
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 <- 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"
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!)
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:
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
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"
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
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
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!
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
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
There is nothing in this data to do this with. But. How would you do it? Two ways.:
I prefer the first
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
(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
## ==========================================
?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.
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()