2 == 1+1[1] TRUE1/2 == .5[1] TRUE2 == 3[1] FALSE2/3 == .6666[1] FALSE3 R Crash Course
Course textbook for F2026 under active development
📄 Download the R code from this chapter
This course is designed for people who already have taken PSCI1800 or have a basic familiarity with R. This chapter serves as a refresher on specific use cases that will be important for the way we will use R in this course.
The textbook for my intro to R class, PSCI 1800, is available here if you feel you need additional reminders about anything I don’t cover here.
There are a few main things that we will be doing in R that are particular to this class that I want to refresh: working with Boolean variables, subsetting, generating simulations, and writing our own functions.
3.1 Boolean Variables
A very helpful feature of R is how it makes use of True or False, or “Boolean”, variables. R can easily test a large number of objects for some sort of logical condition and return a bunch of trues and falses. Further, R very helpfully treats “True” as the number 1, and “False” as the number 0, which will be critical for using R to simulate probability.
If we give R some expressions that can be True or False, it will evaluate them:
A couple of things to note here. First, we use a double equals sign for conditional logic statements. A single equals sign is used for assignment (we use<- as our assignment operator, but you can also use =. I like using the arrow so I don’t get confused about what I’m doing). Second, while TRUE and FALSE are words, they are not in quotation marks. This is not the words “True” and “False”, R is treating these as special values, not just characters.
When we have a vector of information that we create with concatenate (c()) we can also use the double equals sign (or any other boolean operator):
vec <- c(1,2,3,4,5)
vec == 3[1] FALSE FALSE TRUE FALSE FALSEWe can see that R went through and evaluated each of the items in vec as true or false.
Other common boolean operators are greater than and less than:
2 > 3[1] FALSE2 < 3[1] TRUEvec > 3[1] FALSE FALSE FALSE TRUE TRUEvec >= 3[1] FALSE FALSE TRUE TRUE TRUEvec <= 3[1] TRUE TRUE TRUE FALSE FALSEAs I mentioned above, something that ends up being very helpful down the road for us is that R treats Trues as being equal to 1 and Falses be equal to 0. That means that we can quickly figure out how many items in a vector are true:
sum(vec >= 3)[1] 3Moreso, it lets us take advantage of the fact that the mean of a vector of 1’s and 0’s is equal to the probability that the vector is equal to 1.
mean(vec >= 3)[1] 0.6This is going to keep coming up so let’s take a second to prove to ourselves that taking the mean of a variable that is 0 or 1 returns the proportion of 1s.
Consider this vector:
\[ v = \lbrace 0,1,1\rbrace \]
We want to take the mean of this vector. The formula for the mean is:
\[ \bar{v} = \frac{\sum v_i}{n} \]
Which just means: add up all of the values of the vector and divide by the number of values. So for this variable:
\[ \begin{aligned} \bar{v} = \frac{0 + 1+ 1}{3}\\ \bar{v} = \frac{2}{3} \end{aligned} \]
We can see really clearly that what is going to happen here is that the numerator will equal the number of 1s and the denominator will be the total length of the vector. So the mean will return the proportion of entries that equal 1. Or, if we have a boolean variable (which treats T=1 and F=0) the proportion of the entries that are true. We are going to make use of this fact a lot.
Another key boolean operator is the exclamation point, which we use in two distinct ways. First, we use it simply to say “not equal”:
2 == 3 #FALSE[1] FALSE2 != 3 #TRUE[1] TRUEWe can also use it in front of any logical statement to negate the whole thing, turning all TRUE to FALSE and vice versa:
!(2 == 3) #TRUE. essentially the same as the one above[1] TRUEvec == 3[1] FALSE FALSE TRUE FALSE FALSE!(vec == 3)[1] TRUE TRUE FALSE TRUE TRUEWhy is this helpful? Imagine a situation where you have data on all counties in the US and you want to look at all the values except those in Alaska. It is way easier to say “not Alaska” than it is to say “Arizona, Alabama, Arkansas, Connecticut….”.
An important extension to this is the ability to have more than one condition. For this we have both an AND and OR function.
For the AND function we use the ampersand:
(2 == 2) & (3 == 3) #TRUE[1] TRUE(2 == 2) & (3 == 2) #FALSE[1] FALSEThe ‘OR’ operator is denoted by the pipe:
(2 == 2) | (3 == 2) #TRUE[1] TRUE(2 == 8) | (3 == 2) #FALSE[1] FALSEA helpful shortcut for separating AND and OR functions is:
- With an AND function the whole statement is true if and only if all component statements are true.
- With an OR function the whole statement is true if any of the component statements are true.
The final boolean operator we will make use of is %in%, which is maybe the most helpful one. This will test whether things on the left hand side are present in the right hand side. The logical test is being performed on the things on the left hand side.
So consider the confusing and overlapping band memberships of Buffalo Springfield, Crosby Stills and Nash, and Crosby Stills Nash and Young.
buffalo.springfield <- c("stills", "martin","palmer","furay","young")
csn <- c("crosby","stills","nash")
csny <- c("crosby","stills","nash","young")The %in% operator will tell me: which members of Buffalo Springfield were in CSN?
buffalo.springfield %in% csn[1] TRUE FALSE FALSE FALSE FALSEbuffalo.springfield %in% csny[1] TRUE FALSE FALSE FALSE TRUE3.2 Subsetting with Boolean Variables
How do we use these boolean variables to subset? The key insight is that if we put a vector of TRUE and FALSE inside square brackets, R will return the index positions of the TRUE values and not return the index positions of the FALSE values.
Let’s pull in some data on all US counties to demonstrate:
acs <- rio::import("https://github.com/marctrussler/IIS-Data/raw/main/ACSCountyData.csv")
head(acs[,c("county.name","state.abbr","median.income","percent.college")]) county.name state.abbr median.income percent.college
1 Etowah County AL 44023 17.73318
2 Winston County AL 38504 13.53585
3 Escambia County AL 35000 12.65935
4 Autauga County AL 58786 27.68929
5 Baldwin County AL 55962 31.34588
6 Barbour County AL 34186 12.21592If we want the median incomes of counties where more than 40% of the population has a college degree:
head(acs$median.income[acs$percent.college > 40])[1] 63417 75761 84196 116178 78041 93712I would read that expression as: give me the values of median income where percent college is greater than 40.
We can similarly use conditional logic to return rows of the full data frame that meet a certain condition. If I want to return the rows of counties in Pennsylvania:
head(acs[acs$state.abbr == "PA", ]) V1 county.fips county.name state.full state.abbr state.alpha
2245 2245 42015 Bradford County Pennsylvania PA 38
2246 2246 42017 Bucks County Pennsylvania PA 38
2247 2247 42013 Blair County Pennsylvania PA 38
2248 2248 42083 McKean County Pennsylvania PA 38
2249 2249 42085 Mercer County Pennsylvania PA 38
2250 2250 42087 Mifflin County Pennsylvania PA 38
state.icpsr census.region population population.density percent.senior
2245 14 northeast 61304 53.42781 20.26295
2246 14 northeast 626370 1036.30200 17.61275
2247 14 northeast 123842 235.53000 19.79781
2248 14 northeast 41806 42.69485 18.51409
2249 14 northeast 112630 167.45910 20.70674
2250 14 northeast 46362 112.79050 20.77779
percent.white percent.black percent.asian percent.amerindian
2245 96.92353 0.4404280 0.6769542 0.06361738
2246 88.01555 4.0017881 4.6973833 0.16475885
2247 95.50799 1.6214208 0.6992781 0.03391418
2248 94.48405 2.4876812 0.4449122 0.06697603
2249 91.22969 5.6698926 0.6889816 0.10920714
2250 96.88322 0.7915966 0.1898106 0.15745654
percent.less.hs percent.college unemployment.rate median.income gini
2245 11.085675 17.84021 4.808069 51457 0.4310
2246 6.144487 40.46496 4.524462 86055 0.4480
2247 9.070472 20.80933 4.961147 47969 0.4448
2248 8.908716 18.14864 6.829542 46953 0.4475
2249 10.161785 22.55270 5.271703 48768 0.4381
2250 14.792334 12.47843 3.598866 47526 0.4144
median.rent percent.child.poverty percent.adult.poverty
2245 711 17.012323 11.132679
2246 1207 7.358998 5.743422
2247 711 19.575034 14.600937
2248 644 27.077236 16.563244
2249 679 24.827619 13.507259
2250 665 23.594624 13.347636
percent.car.commute percent.transit.commute percent.bicycle.commute
2245 90.12223 0.2475153 0.27036290
2246 88.75924 3.4039931 0.16052502
2247 91.58161 1.1913318 0.18328182
2248 91.10273 0.9152430 0.17750166
2249 90.04547 0.4175279 0.08887970
2250 90.67850 0.3038048 0.05786758
percent.walk.commute average.commute.time percent.no.insurance
2245 3.625148 23 6.692875
2246 1.717094 30 4.212366
2247 3.056494 20 5.240548
2248 3.372532 22 4.568722
2249 3.577925 21 6.158217
2250 2.956069 24 12.339847Remember the use of the comma here! What we are subsetting is a data frame with 2 dimensions, so we have to give information on which columns we want to return. By putting a comma and leaving the second part blank, we are saying “give me all the columns”.
We can string together booleans with & and | to subset on multiple conditions:
#Counties in Pennsylvania with over 40% college educated:
head(acs[acs$state.abbr == "PA" & acs$percent.college > 40, ]) V1 county.fips county.name state.full state.abbr state.alpha
2246 2246 42017 Bucks County Pennsylvania PA 38
2252 2252 42091 Montgomery County Pennsylvania PA 38
2260 2260 42027 Centre County Pennsylvania PA 38
2261 2261 42029 Chester County Pennsylvania PA 38
2289 2289 42003 Allegheny County Pennsylvania PA 38
state.icpsr census.region population population.density percent.senior
2246 14 northeast 626370 1036.3020 17.61275
2252 14 northeast 821301 1700.5640 17.00095
2260 14 northeast 161443 145.4622 13.26970
2261 14 northeast 517156 689.0596 15.44138
2289 14 northeast 1225561 1678.5810 18.11171
percent.white percent.black percent.asian percent.amerindian
2246 88.01555 4.001788 4.697383 0.1647588
2252 79.41607 9.079375 7.476552 0.1271154
2260 87.51200 3.961770 6.345893 0.0953897
2261 85.33711 5.937473 5.177741 0.1077044
2289 80.11409 12.883161 3.654734 0.1104800
percent.less.hs percent.college unemployment.rate median.income gini
2246 6.144487 40.46496 4.524462 86055 0.4480
2252 5.781227 48.70617 4.619781 88166 0.4667
2260 5.673250 44.73959 4.253275 58055 0.4696
2261 6.749068 51.81124 4.356723 96726 0.4572
2289 5.712555 40.70367 5.286755 58383 0.4869
median.rent percent.child.poverty percent.adult.poverty
2246 1207 7.358998 5.743422
2252 1253 7.293578 6.084560
2260 966 11.803944 22.887424
2261 1287 8.224267 6.776429
2289 865 16.381297 11.940936
percent.car.commute percent.transit.commute percent.bicycle.commute
2246 88.75924 3.403993 0.1605250
2252 85.07758 5.457351 0.2106967
2260 77.86436 5.472074 2.0904255
2261 86.12955 2.732424 0.1857280
2289 79.94947 9.523755 0.5728954
percent.walk.commute average.commute.time percent.no.insurance
2246 1.717094 30 4.212366
2252 2.276138 29 3.767316
2260 8.864362 20 5.295987
2261 2.747870 28 5.890486
2289 4.109690 27 4.278286The %in% operator is often the cleanest way to subset to multiple categories. If I want counties in any of the northeastern states:
northeast <- c("PA","NY","NJ","CT","MA","VT","NH","ME","RI")
head(acs[acs$state.abbr %in% northeast, ]) V1 county.fips county.name state.full state.abbr state.alpha
309 309 9005 Litchfield County Connecticut CT 7
310 310 9007 Middlesex County Connecticut CT 7
311 311 9003 Hartford County Connecticut CT 7
312 312 9001 Fairfield County Connecticut CT 7
313 313 9013 Tolland County Connecticut CT 7
314 314 9015 Windham County Connecticut CT 7
state.icpsr census.region population population.density percent.senior
309 1 northeast 183031 198.8359 19.91794
310 1 northeast 163368 442.3669 19.00739
311 1 northeast 894730 1217.3820 16.45860
312 1 northeast 944348 1511.0420 15.18709
313 1 northeast 151269 368.6320 14.81136
314 1 northeast 116538 227.2013 15.66613
percent.white percent.black percent.asian percent.amerindian
309 92.92579 1.927542 1.854331 0.19450257
310 88.66118 5.297243 3.072817 0.16159835
311 71.47721 13.646910 5.205705 0.32914958
312 73.24016 11.353336 5.286928 0.24757822
313 88.35650 3.054823 4.651978 0.04429196
314 88.86972 2.268788 1.384956 0.61010143
percent.less.hs percent.college unemployment.rate median.income gini
309 7.410651 35.07926 5.234069 78314 0.4479
310 5.823624 41.52384 4.693489 84761 0.4427
311 10.469985 37.75315 6.543712 72321 0.4716
312 10.110784 47.35474 6.955179 92969 0.5432
313 5.439658 41.15365 5.248147 84916 0.4315
314 11.350250 24.14521 6.283623 64774 0.4151
median.rent percent.child.poverty percent.adult.poverty percent.car.commute
309 1036 7.701885 7.096809 89.94674
310 1162 8.338651 7.000711 88.49398
311 1076 15.441200 10.286984 89.51368
312 1470 10.954530 8.405622 80.34434
313 1119 4.802281 8.191013 86.55236
314 895 13.548696 10.178451 91.43660
percent.transit.commute percent.bicycle.commute percent.walk.commute
309 1.3545533 0.13105198 2.250949
310 1.4582524 0.06971269 2.621654
311 3.2843624 0.22519651 1.936371
312 10.2409910 0.18379396 2.528293
313 1.9894721 0.22062239 4.602106
314 0.5409203 0.21250439 2.736214
average.commute.time percent.no.insurance
309 28 4.172517
310 26 3.284609
311 23 4.285092
312 31 8.589842
313 26 3.009870
314 26 3.872557Now all of these examples have been super obvious. Things we could eyeball. Try to think about how helpful these tools are going to be when we have a dataset with thousands of rows.
3.3 For Loops
One of the most flexible tools in all coding languages is for() loops. In short: for() loops allow us to repeat the same code many, many times. There are two main use cases in this course. The first is when we need to perform the same operation on different clusters in the data (below, for example, we’ll compute the mean population for the counties in each state). The second is when we want to use R to simulate probability, which we’ll get to in the next section.
Previously we saw that we can use square brackets to give us conditional means. For example we can easily get the mean population for all counties:
mean(acs$population)[1] 102769.9And if we want the mean population for just the counties in Alabama:
mean(acs$population[acs$state.abbr=="AL"])[1] 72607.16OK, but what if we wanted to do the same thing for all 50 states?
We could think about brute forcing it, subbing in the right two letter abbreviation and repeating the code 50 times:
mean(acs$population[acs$state.abbr=="AK"])[1] 25466.07mean(acs$population[acs$state.abbr=="AZ"])[1] 463112.3mean(acs$population[acs$state.abbr=="AR"])[1] 39875.61mean(acs$population[acs$state.abbr=="CA"])[1] 674978.6#
#
#
mean(acs$population[acs$state.abbr=="WI"])[1] 80255.47mean(acs$population[acs$state.abbr=="WY"])[1] 25297.22This is unwieldy! And I didn’t even do all of the states!
But notice something about the above block of code: most of the command on each line is exactly the same. All that we are changing in every line is what is in the quotations. We are going from AL, to AZ, to AR…. all the way to WY.
So all we need is a method that allows us to run mean(acs$population[acs$state.abbr=="XX"]) 50 times, each time through changing XX to the next two-letter state abbreviation.
The way that we will accomplish this is with a for() loop.
Here is the basic setup for a for() loop:
for(i in 1:3){
#CODE HERE
}We can look at the for loop as having two parts. Everything from the start to the opening squiggly bracket are the instructions for the loop. It basically tells the loop how many times to repeat the same thing. “i” in this case is the variable that is going to change its value each time through the loop. Having this variable is what allows the loop to do something slightly different each time through the loop. The second part is whatever you put in between the squiggly brackets. This is the code that will be repeated again, and again, the specified number of times.
Let’s put a simple bit of code in the loop to see how these work. I’m just going to put the command print(i) in the loop so that all the loop does is print the current value of i.
for(i in 1:3){
print(i)
}[1] 1
[1] 2
[1] 3What a loop does is the following: on the first time through the loop it sets our variable (i) to the number 1. Literally. Anywhere we have the variable i, R will replace it with the integer 1. Once all the code has been executed and R reaches the closing squiggly bracket it goes back up to the start. Now it will set i equal to 2. Anywhere in our code where there is currently the variable “i”, R will literally treat this as the integer 2. After it reaches the end of the code again, it will go back to the beginning and run the code again, this time making “i” equal to 3.
So the result of running this code is to simply print out the numbers 1, then 2, then 3.
To make this crystal clear, here is what code is literally being executed by R here. It is just the print command 3 times, the first time changing the variable i to the integer 1, the second time changing it to 2, and the third time changing it to 3.
print(1)[1] 1print(2)[1] 2print(3)[1] 3There is nothing special about “i” — we can name the variable that changes every time through the script anything we want. So we could edit this same code to be for(batman in 1:3), and have the variable be called batman. And if we run this we get the exact same result:
for(batman in 1:3){
print(batman)
}[1] 1
[1] 2
[1] 3While we can use anything, it is convention to use the letter “i” as the indexing variable in a for loop. This is true across coding platforms so it’s a good idea to just use that.
Now let’s use our knowledge of for() loops to complete the task we set out: calculating the mean population for the counties in every state.
The way we are going to do that is to first create and save a vector of all the state abbreviations. Using the unique() command, I can create a vector of all the state abbreviations. If we look at this, we now have a vector of all the possible state abbreviations.
states <- unique(acs$state.abbr)
head(states)[1] "AL" "AK" "AZ" "AR" "CA" "CO"length(states)[1] 51Remember that by using square brackets we can select a certain value out of a vector. So the first entry in states is AL, the 13th entry is ID, and the 27th entry is Montana.
states[1][1] "AL"states[13][1] "ID"states[27][1] "MT"So instead of subsetting to a particular two letter abbreviation, now we can insert this vector and have the loop sequentially select each of the items in it each time through the loop.
for(i in 1:3){
mean(acs$population[acs$state.abbr==states[i]])
}So instead of acs$state.abbr=="AL", we now have acs$state.abbr==states[i].
Again, let’s think about what R is going to do in this code. The first time through the loop it will change i to the integer 1, then the integer 2, then the integer 3. So this successively subsets to the first state, the second state, then the third state.
#What it will literally do
mean(acs$population[acs$state.abbr==states[1]])[1] 72607.16mean(acs$population[acs$state.abbr==states[2]])[1] 25466.07mean(acs$population[acs$state.abbr==states[3]])[1] 463112.3But hold on, right now we are only going from 1 to 3, but there are 50 states + DC. We could just write in 1:51, but a better coding practice is to set the end point of the loop to be the length of the vector states, regardless of what the length is. That way if we ever change the vector in the future our loop still works.
for(i in 1:length(states)){
mean(acs$population[acs$state.abbr==states[i]])
}Now let’s run it and see what happens: nothing!
R, in fact, did run the command in the brackets 51 times, each time changing i from 1 to 2 to 3, effectively looping over the 50 states + DC. But unless we explicitly tell R to save the results of what it is doing, it will do all of this in the background when running code inside a loop.
So the last step in this process is to create an empty vector to store the results. We need to create a vector of NAs that is the same length as states.
state.pop.means <- rep(NA, length(states))
for(i in 1:length(states)){
state.pop.means <- mean(acs$population[acs$state.abbr==states[i]])
}
state.pop.means[1] 25297.22Huh. Well that didn’t work. We were expecting to get 51 state means in this vector, but instead we got just one. Why did that happen?
Well, the first time through the loop we told R: assign the mean of Alabama to the object state.pop.means, and it did. Then we said, OK, nevermind, assign the mean of Alaska to state.pop.means. Then we said, OK, nevermind, assign the mean of Arkansas to state.pop.means. So this last number is just the last mean we calculated (the average population of counties in Wyoming).
The last step here is to tell R to save the results in the right “slot” of state.pop.means. We do this by using our variable i to index the results vector as well. Each time through the loop we save the results to the first slot in state.pop.means, and then the next time in the second slot, etc.
state.pop.means <- rep(NA, length(states))
for(i in 1:length(states)){
state.pop.means[i] <- mean(acs$population[acs$state.abbr==states[i]])
}
head(cbind(states, state.pop.means)) states state.pop.means
[1,] "AL" "72607.1641791045"
[2,] "AK" "25466.0689655172"
[3,] "AZ" "463112.333333333"
[4,] "AR" "39875.6133333333"
[5,] "CA" "674978.620689655"
[6,] "CO" "86424.078125" Each time through the loop our variable i gets changed to a new value. The first time through the loop, changing i to the integer 1 accomplishes two tasks. First, it tells R which slot in the empty vector state.pop.means to save the answer. Second, it tells R which entry in the vector states to access in order to subset our data.
So always try to remember that nothing more complicated than this is happening inside of a loop. And you can always just enter the integers yourself to get a sense of what is going on.
3.4 R for Simulations
The other major use case for for() loops is to simulate probability. A helpful definition of the probability of an event is the frequency that an outcome occurs if the event is repeated a large number of times.
So, for example, the probability of a coin coming up heads on a fair coin is 50%. That does not mean that if we flip a tail we will get heads next. It means that if we flip a coin a large number of times, the number of heads will approach 50%.
We can write code that simulates one coin flip:
coin <- c("H","T")
sample(coin, 1)[1] "T"The sample() command is a random command. We say what we want to sample in the first argument, and how many times we want to sample it in the second argument. Every time we run the command we get something different. (Well, here we may not because there are only two things we are sampling. But every time we run it is a new random draw.)
sample(coin, 1)[1] "H"sample(coin, 1)[1] "H"sample(coin, 1)[1] "T"sample(coin, 1)[1] "H"sample(coin, 1)[1] "H"Because sample() produces a different result each time we run it, we’ll frequently want to make our code reproducible. Setting a seed at the top of a simulation makes sure everyone (including future-you) gets the same “random” numbers each time the code is run. Throughout this course I set the seed to Penn’s zip code:
set.seed(19104)
sample(coin, 1)[1] "T"Above we defined probability as the frequency of an outcome when we repeat an event a large number of times. With R we can just repeat an event a large number of times using a loop:
for(i in 1:1000){
sample(coin, 1)
}Like with all loops, this code did indeed run that sample() command 1000 times, but if we don’t save what we are doing we don’t see any result. Instead, let’s save the outcome of each loop to a results vector. The best practice for this is to pre-allocate a vector of NAs of the length we want, and then fill in each position as we go:
result <- rep(NA, 1000)
for(i in 1:1000){
result[i] <- sample(coin, 1)
}
head(result)[1] "T" "T" "H" "T" "H" "T"What percentage of these results are Heads? We can figure that out, but I’m going to propose that we edit our code slightly to make better use of what we just learned about boolean variables. All we want to know is if the coin is heads or not, so let’s evaluate that each time we sample:
result <- rep(NA, 1000)
for(i in 1:1000){
result[i] <- sample(coin, 1) == "H"
}
head(result)[1] FALSE TRUE FALSE FALSE FALSE FALSEWe now have 1000 trues or falses. Because R treats True as 1 and False as 0, we can determine the proportion of Trues by taking the mean of this variable:
mean(result)[1] 0.503As expected it’s approximately 50%.
Now that was a pretty obvious example. The nice thing about this method is that it can answer probability questions as long as we can generate one answer. We just need to put that code in a loop and determine how frequently a certain answer occurs in the long run.
So: if I roll 3 dice and flip 2 coins, how often does the dice add up to a number greater than 12 and we get two heads?
Let’s write code to do this event once:
dice <- c(1,2,3,4,5,6)
#Using the replace=T option because each time we sample from the dice
#we want all sides to be available
sum(sample(dice, 3, replace=T))[1] 9#Same with the coin
sample(coin, 2, replace=T)[1] "H" "H"This code accomplished this task, but we want to know if the conditions we set have been met.
#This one is easy: is the sum greater than 12
sum(sample(dice, 3, replace=T))>=12[1] TRUE#For the coin we want to see if each entry is equal to "H", but only return one true
#if both of them are H. We can do that with all()
sample(coin,2, replace=T)=="H"[1] TRUE TRUEall(sample(coin,2, replace=T)=="H")[1] FALSE#And then we want to know if both conditions are met:
sum(sample(dice, 3, replace=T))>=12 & all(sample(coin,2, replace=T)=="H")[1] FALSEOK, we’ve written code to do this once. All we need to do is repeat this a large number of times to determine the probability of this compound event:
result <- rep(NA, 1000)
for(i in 1:1000){
result[i] <- sum(sample(dice, 3, replace=T))>=12 & all(sample(coin,2, replace=T)=="H")
}
mean(result)[1] 0.103Around 10% of the time do we get a sum over 12 and two heads. Not sure what to do with that information tbh.
This method can be applied to any sort of probability problem, but where we are particularly interested in it will be helping us to understand the principles of sampling. When we take a survey we get one sample of many possible samples we could get. Each of those samples will be a bit different, and statistics helps us understand just how different those samples will be. Much more on this to come!
One principle of sampling that we will learn is that confidence intervals contain the true population parameter 95% of the time. If the true proportion we are measuring with our sample is 30%, if we repeatedly take samples and form confidence intervals, 95% of them will contain 30%.
Here is simulating that once:
coin <- c(0,1)
samp <- sample(coin, 1000, replace=T, prob=c(.7,.3))
#Lower bound less than .3 and upper bound greater than .3?
t.test(samp)$conf.int[1] < .3 & t.test(samp)$conf.int[2] > .3[1] TRUETo prove this axiom of probability, we can repeat this process many times and store the results:
result <- rep(NA, 1000)
for(i in 1:1000){
samp <- sample(coin, 1000, replace=T, prob=c(.7,.3))
#Lower bound less than .3 and upper bound greater than .3?
result[i] <- t.test(samp)$conf.int[1] < .3 & t.test(samp)$conf.int[2] > .3
}
mean(result)[1] 0.95795% of confidence intervals contained the truth (30%). This is where we are headed! For right now, don’t worry too much about what a confidence interval is or how t.test() works. The point is the pattern — simulate an event once, put tha code in a loop, and take the mean of the boolean result to get a probability.
3.5 table() and prop.table()
When we look at categorical variables we frequently want to know the frequency and proportion of each category. This is what table() and prop.table() are for.
Let’s load some data from the American National Election Study (ANES):
anes <- read.csv("https://raw.githubusercontent.com/marctrussler/IIS-Data/main/ANES2020Clean.csv")If we want the raw counts of a variable we use table():
table(anes$race)
Asian/Hawaiian/Pacific-Islander Black, non-Hispanic
284 726
Hispanic Multiple races, non-Hispanic
762 271
Native American White, non-Hispanic
172 5963 Wrapping that in prop.table() gives us relative frequencies (proportions) instead of counts:
prop.table(table(anes$race))
Asian/Hawaiian/Pacific-Islander Black, non-Hispanic
0.03472732 0.08877476
Hispanic Multiple races, non-Hispanic
0.09317682 0.03313769
Native American White, non-Hispanic
0.02103204 0.72915138 We can also make two-way tables by giving table() two variables. I like to label the variables in the command to make things more clear:
table(race = anes$race, gender = anes$gender) gender
race Female Male
Asian/Hawaiian/Pacific-Islander 133 151
Black, non-Hispanic 466 255
Hispanic 404 357
Multiple races, non-Hispanic 156 114
Native American 78 90
White, non-Hispanic 3177 2763And we can pass that to prop.table() too. By default, this gives us the percentage of the total sample in each cell:
prop.table(table(race = anes$race, gender = anes$gender)) gender
race Female Male
Asian/Hawaiian/Pacific-Islander 0.016331041 0.018541257
Black, non-Hispanic 0.057220039 0.031311395
Hispanic 0.049607073 0.043835953
Multiple races, non-Hispanic 0.019155206 0.013998035
Native American 0.009577603 0.011051081
White, non-Hispanic 0.390103143 0.339268173We can add an option to prop.table to calculate either row or column percentages. 1 will calculate row percentages and 2 will calculate column percentages. (This is just like square brackets where the first spot is for row information and the second spot is for column information!)
Row percentages — of people of each race, what percent are Female vs Male?
prop.table(table(race = anes$race, gender = anes$gender), 1) gender
race Female Male
Asian/Hawaiian/Pacific-Islander 0.4683099 0.5316901
Black, non-Hispanic 0.6463245 0.3536755
Hispanic 0.5308804 0.4691196
Multiple races, non-Hispanic 0.5777778 0.4222222
Native American 0.4642857 0.5357143
White, non-Hispanic 0.5348485 0.4651515Column percentages — of people of each gender, what percent are each race?
prop.table(table(race = anes$race, gender = anes$gender), 2) gender
race Female Male
Asian/Hawaiian/Pacific-Islander 0.03013140 0.04048257
Black, non-Hispanic 0.10557318 0.06836461
Hispanic 0.09152696 0.09571046
Multiple races, non-Hispanic 0.03534209 0.03056300
Native American 0.01767105 0.02412869
White, non-Hispanic 0.71975532 0.74075067These three tables answer three genuinely different questions. When we get to conditional probability in the next chapter you will see that these are actually the joint, and two different conditional, probabilities.
3.6 Writing Functions
Functions take a specific set of inputs and produce a specific set of outputs.
That might sound confusing, but you have already been using functions constantly. Every time you write mean(), sum(), length(), table(), or head() you are calling a function. Someone, at some point, wrote the code that makes mean() work, and packaged it up so that all you have to do is give it a vector (an input) and it gives you back the mean (the output).
In this section we are going to learn how to write our own functions. Here’s the basic idea: if you find yourself copy-pasting the same chunk of code and just changing one or two things each time, you should probably write a function instead.
A function in R has three parts: a name, arguments (the inputs), and a body (the code that runs). Here is the simplest possible function:
say.hello <- function(){
print("Hello!")
}Let’s break this down. say.hello is the name – this is what we will type to call the function, just like we type mean to call the mean function. The function() part tells R “I am defining a function.” The curly braces {} contain the body – the code that runs when you call the function. In this case, all it does is print “Hello!”.
To call the function:
say.hello()[1] "Hello!"Note the parentheses. If you forget them, R doesn’t run the function – it just tells you that the function exists:
say.hellofunction ()
{
print("Hello!")
}That’s not what we want. Always include the parentheses when you want to actually run a function.
A function that does the same thing every time isn’t very useful. The real power comes from giving functions arguments – inputs that change what the function does.
To continue our example, let’s write our own function for calculating the mean:
new.mean <- function(vec){
sum(vec)/length(vec)
}The input into this function is vec, which will be a vector of numbers. The function sums that vector and divides by its length, giving us the mean.
Now we can set vec equal to any vector of numbers and it will calculate the mean:
x <- 1:175
new.mean(vec=x)[1] 88R will automatically treat what you put in parentheses as the first argument, so we can shortcut this to:
new.mean(x)[1] 88#To show it works with other things
y <- 4:17
new.mean(y)[1] 10.5We can put as many arguments as we want into a function. Let’s do a more useful example. I have a locker in Pottruck and the code has 3 numbers picked from a dial of 40. Not every locker is like this though: some dials have 30 numbers, some have 50, and codes can be shorter or longer. Let’s write a function that generates a random combination for any size dial and any length of code.
random.combo <- function(dial, entries){
sample(1:dial, entries, replace=F)
}The function takes two arguments: dial is the number of positions on the dial, and entries is how many numbers are in the combination. Inside, we just use sample() to draw entries numbers from 1:dial, without replacement (because a real lock code doesn’t repeat numbers).
Now we can generate a combination for any lock we like:
set.seed(19104)
#My Pottruck locker
random.combo(dial=40, entries=3)[1] 28 21 37#Some fancy safe
random.combo(dial=100, entries=5)[1] 51 78 20 97 86#The tiniest lock in the world
random.combo(dial=10, entries=2)[1] 7 2Every time we call the function we get a new random combination, because sample() is a random command.
We can use this function to estimate the probability of correctly guessing a locker combination by pure chance. What are the odds that someone who has one guess at my Pottruck locker gets in?
We can first set a true combination:
set.seed(19104)
truth <- random.combo(dial=40, entries=3)
truth[1] 28 21 37And then generate guesses to see whether each position matches the truth:
guess <- random.combo(dial=40, entries=3)
guess == truth[1] FALSE FALSE FALSEThe all() function returns TRUE only if every entry of a vector is true. That is exactly the condition we want – every position of the guess has to match every position of the truth:
all(guess == truth)[1] FALSENow we do this a large number of times, saving whether each attempt matches, and take the mean of the result vector to get the probability:
set.seed(19104)
truth <- random.combo(dial=40, entries=3)
result <- rep(NA, 100000)
for(i in 1:100000){
guess <- random.combo(dial=40, entries=3)
result[i] <- all(guess == truth)
}
mean(result)[1] 1e-05We get a tiny number, meaning correctly guessing a lock code is very unlikely. Reassuring! We can also confirm this with math: there are \(40 \times 39 \times 38 = 59{,}280\) possible combinations, and only one is correct, so \(P(\text{correct guess}) = 1/59{,}280 \approx 0.0000169\). Note that even doing 100,000 simulations, we only expect around 1 or 2 matches, so simulation values here are going to jump around a bit.
And because we wrote the function generally, we can immediately answer the same question for a lock of any size, or a code of any length, just by changing the arguments.
3.7 if() Statements
Sometimes we want to say: if a certain condition holds, then run this code. We can do this via if statements. Similar to when we do conditional logic for a single line of code, we start with a logical statement – but in this case a statement that produces a single T or F, and then R will only run the code if the statement is True.
If statements look a lot like for loops, but in the opening bracket we put a conditional statement that is either T or F. If the statement is T R will run what is in the squiggly brackets. If it is false it won’t:
if(2+2==4){
print("I ran code chunk 1")
}[1] "I ran code chunk 1"if(2+2==5){
print("I ran code chunk 2")
}When we run the first bit of code it runs what is in the squiggly brackets, when we run the second chunk of code it does not run what is in the squiggly brackets.
Right now our code just says: if the statement is true run this code, if it is not, move on. But many times, we want to say: if this statement is true run this code, if it is not, run this other code. In that case we can add an else to our code:
if(2+2==5){
print("I ran code chunk 2")
} else {
print("I ran code chunk 3")
}[1] "I ran code chunk 3"Because the logical statement at the start is false, R runs the code following the else.
We aren’t limited to just two conditions, we can use else if to add a second logical statement:
if(2+2==5){
print("I ran code chunk 1")
} else if(2+2==4){
print("I ran code chunk 2")
} else {
print("I ran code chunk 3")
}[1] "I ran code chunk 2"Where this becomes useful in this course is inside a loop, when the code we want to run depends on the outcome of the current iteration. For example, when we simulate the Monty Hall problem in the probability chapter, what a player gets from “switching” depends on whether they initially picked a car or a goat. An if()/else inside the loop lets us handle those two cases differently. We’ll return to this then.