Math 189/289 Lab 3
Alex Hanbo Li
Dataframe
R has a working directory. This should be the main location to store external data files, as well as other code files.
getwd()One can also set the working directory to a desired location.
setwd() # desired directory as argument
getwd()If you are not sure where the source file is located, you can use “Session –> Set Working Directory –> To Source File Location”.
First, let’s load in the video game data for homework 2.
setwd("~/Dropbox/Math189Winter2017/lab3") # change to your own working directory
data <- read.table("videodata.txt", header=TRUE)
head(data)If a question was not answered or improperly answered, it was decoded as 99. There are total number of 65 such fields.
data[data == 99] <- NA
sum(is.na(data))[1] 65attach() and detach() functions
The $ notation for list components is not always very convenient. Sometimes it is more convenient to temporarily make the components of data frame visible as variables under their component name. We can use attach() function in R to do this.
For example, if we want to access the age variable, before attaching the data frame data, we can only access by:
data$age # returns an array [1] 19 18 19 19 19 19 20 19 19 19 20 19 19 18 18 19 21 20 18 19 20 24 19 19 21 20 22 18 19 20 19 19 19 19 19 18 19 20 19
[40] 20 19 19 19 19 19 19 19 19 20 19 18 20 19 18 19 21 18 20 19 19 19 21 20 19 18 19 18 20 19 19 19 20 19 23 19 20 19 19
[79] 25 19 20 19 19 21 18 19 19 20 33 19 19data['age'] # return a dataframeIf we directly type in age, an error will be returned.
ageError: object 'age' not foundAfter attaching the dataframe, we have an individual variable called age.
attach(data)The following object is masked _by_ .GlobalEnv:
ageage [1] 1000 18 19 19 19 19 20 19 19 19 20 19 19 18 18 19 21 20 18 19 20 24 19
[24] 19 21 20 22 18 19 20 19 19 19 19 19 18 19 20 19 20 19 19 19 19 19 19
[47] 19 19 20 19 18 20 19 18 19 21 18 20 19 19 19 21 20 19 18 19 18 20 19
[70] 19 19 20 19 23 19 20 19 19 25 19 20 19 19 21 18 19 19 20 33 19 19However, now any change on age will only affect the new variable age in the working directory, not the age attribute of the data frame data.
age[1] <- 1000
age [1] 1000 18 19 19 19 19 20 19 19 19 20 19 19 18 18 19 21 20 18 19 20 24 19
[24] 19 21 20 22 18 19 20 19 19 19 19 19 18 19 20 19 20 19 19 19 19 19 19
[47] 19 19 20 19 18 20 19 18 19 21 18 20 19 19 19 21 20 19 18 19 18 20 19
[70] 19 19 20 19 23 19 20 19 19 25 19 20 19 19 21 18 19 19 20 33 19 19data$age [1] 19 18 19 19 19 19 20 19 19 19 20 19 19 18 18 19 21 20 18 19 20 24 19 19 21 20 22 18 19 20 19 19 19 19 19 18 19 20 19
[40] 20 19 19 19 19 19 19 19 19 20 19 18 20 19 18 19 21 18 20 19 19 19 21 20 19 18 19 18 20 19 19 19 20 19 23 19 20 19 19
[79] 25 19 20 19 19 21 18 19 19 20 33 19 19To detach a data frame, use function detach().
search() [1] ".GlobalEnv" "data" "tools:rstudio" "package:stats" "package:graphics"
[6] "package:grDevices" "package:utils" "package:datasets" "package:methods" "Autoloads"
[11] "package:base" ls(2) [1] "age" "busy" "cdrom" "educ" "email" "freq" "grade" "home" "like" "math" "own" "sex" "time" "where"
[15] "work" detach(data)
search() [1] ".GlobalEnv" "tools:rstudio" "package:stats" "package:graphics" "package:grDevices"
[6] "package:utils" "package:datasets" "package:methods" "Autoloads" "package:base" More on estimating the distribution of data
Histogram
summary(data$age) Min. 1st Qu. Median Mean 3rd Qu. Max.
18.00 19.00 19.00 19.52 20.00 33.00 hist(data$age)
hist(data$age, breaks = 25, probability = TRUE, density = 20, col = 3, border = 3)
lines(density(data$age, bw = 0.25), col = 1)lines(density(data$age, bw = 0.5), col = 2)
rug(data$age)
Empirical cumulative distribution function (ecdf())
plot(ecdf(data$age), do.points = F, verticals = T)
x <- seq(min(data$age)-4, max(data$age), length.out = 1000)
lines(x, pnorm(x, mean=mean(data$age), sd=sd(data$age)), lty=3, col = 2)
Let’s compare it with a Normal distribution with the same sample mean and variance (red dotted line above).
Of course, one would check the Normality more carefully using Q-Q plot. We see a very long right tail.
par(pty = 's') # plot in a square
qqnorm(data$age)
qqline(data$age)
More formally, we can do a Shapiro-Wilk test. In this case, the p-value is significantly small and we can reject the Normality hypothesis \(H_0\).
shapiro.test(data$age)
Shapiro-Wilk normality test
data: data$age
W = 0.51412, p-value = 7.666e-16Sampling and Bootstrap
Sampling
Sampling (with or without replacement) can be done in R using the function sample().
sample(1:10, size = 10, replace = F) [1] 8 9 2 7 3 1 6 5 4 10sample(1:10, size = 10, replace = T) [1] 2 1 5 10 8 6 4 2 7 5Bootstrap
Suppose we want to find the point estimate of the proportion of male in the population and its confidence interval. Note that the sex variable is 1 if the student is male and 0 if female. Then the point estimation of the proportion of male is:
male.percentage <- mean(data$sex)
male.percentage[1] 0.5824176Now we also want to have a confidence interval of this estimator. However, clearly the distribution of sex variable is not Normal, it is a Bernoulli random variable. We know our data were drawn from a population with size \(N = 314\). Hence, we first create a bootstrap population of this size by repeating every sample for \(\frac{314}{91} = 3.45\) times. Here, we’ll just specify the parameter length.out to be 314.
boot.population <- rep(data$sex, length.out = 314)
length(boot.population)[1] 314Then we will choose \(n = 91\) samples from the Bootstrap population and call this a Bootstrap sample.
sample1 <- sample(boot.population, size = 91, replace = FALSE)Continue this procedure until we have 400 Bootstrap samples.
set.seed(189289)
B = 400 # the number of bootstrap samples we want
boot.sample <- array(dim = c(B, 91))
for (i in 1:B) {
boot.sample[i, ] <- sample(boot.population, size = 91, replace = FALSE)
}Then we can calculate the sample mean for each Bootstrap sample (i.e. each row of the Bootstrap sample matrix).
boot.mean <- apply(X = boot.sample, MARGIN = 1, FUN = mean)
head(boot.mean)[1] 0.5054945 0.6153846 0.6263736 0.6043956 0.4835165 0.5824176Let’s see the histogram of these Bootstrap sample means.
hist(boot.mean, breaks = 20, probability = TRUE, density = 20, col = 3, border = 3)
lines(density(boot.mean, adjust = 2), col = 2)
Check Normality by Q-Q plot and Shapiro-Wilk test.
par(pty = 's')
qqnorm(boot.mean)
qqline(boot.mean)
shapiro.test(boot.mean)
Shapiro-Wilk normality test
data: boot.mean
W = 0.99244, p-value = 0.04051So we can accept that the sample mean follows a Normal distribution. Then we can construct \(95\%\) confidence intervals.
boot.sd <- sd(boot.mean)
male.percentage + c(-1, 1)*1.96*boot.sd[1] 0.4969015 0.6679337