Math 189/289 Lab 4

Create Data

Let’s say today we are going to work with some simulated data. We generate \(N\) numbers from bernoulli distribution with success rate 0.3, and take them as our population.

N <- 1000
data.population <- rbinom(n=N, size=1, prob=0.3)

One can give a story line to this data. If put in the same context as the video games dataset that we have for homework, for example, this could be the response of the students in the whole school whether they played games in the past week or not. Thus, the school has 1000 students in total. 1 indicates the student played, and 0 indicates did not play.

Then in that sense, we might not have observed every one of these responses as in for the survey data in the video games dataset. Thus, we sample, say \(n\), observations from the population, and take them as our sample units.

n <- 300
ind.sample <- sample.int(n=N, size=n)
data.sample <- data.population[ind.sample]

Sample Statistics

Now we will stick with the story line. Once we have the sample data, we can compute a point estimate, as well as an interval estimate for the fraction of the students who played games in the past week or not.

A point estimate in this case is just the mean of the sample,

mean.sample <- mean(data.sample)
mean.sample

[1] 0.3033333

To get an interval estimate for the fraction, we follow the derivation in the lecture slides, the interval estimate is then given by, \[ \left(\bar x - 1.96 \sqrt{\frac{\bar x (1 - \bar x)}{n-1} \frac{N - n}{N}}, \bar x + 1.96 \sqrt{\frac{\bar x (1 - \bar x)}{n-1} \frac{N - n}{N}} \right), \] where \(\bar x\) indicates the sample mean, \(N\) indicates the population size, and \(n\) indicates the sample size. Thus, an interval estimate is then given as the following.

width <- 1.96 * sqrt(mean.sample*(1-mean.sample)*(N-n)/((n-1)*N))
int.sample <- c(mean.sample - width, mean.sample + width)
int.sample

[1] 0.2597378 0.3469289

Bootstrap Estimate

Another popular method is using the idea of bootstrap. We need a bootstrap population to start with. With the population size of \(N=1000\), and sample size of \(n=300\), approximately each sample occurs about 3 times in the bootstrap population. One can simply duplicate each observation in the sample 3 times, and treat the resulting sample as the bootstrap population.

One can also sample from the sample with replacement.

ind.boot <- sample.int(n, size=N, replace=TRUE)
data.boot <- data.sample[ind.boot]

With the bootstrap population, we are now ready to generate bootstrap sample means. We will take, say 2000, random samples of size \(n\).

B <- 2000
boot.sample.mean <- rep(NA, B)
for(i in 1:B){
  ind <- sample.int(N, size=n, replace=FALSE)
  boot <- data.boot[ind]
  boot.sample.mean[i] <- mean(boot)
}

Let’s take a look at the distribution of bootstrap sample means.

hist(boot.sample.mean)

The point estimate is the mean of bootstrap sample means.

mean.boot <- mean(boot.sample.mean)
mean.boot

[1] 0.300755

There are two approaches one can take from here. We can derive a confidence interval using the bootstrap sample means, as given by \[ \left( \bar x - 1.96 s, \bar x + 1.96s \right) \]

s <- sd(boot.sample.mean)
int.boot <- c(mean.boot - 1.96*s, mean.boot + 1.96*s)
int.boot

[1] 0.2580095 0.3435005

Another approach to derive an interval estimate using the bootstrap sample means, one can simply extract the 0.025-quantile and 0.975 quantile of the bootstrap sample means and arrive at an interval estimate.

int.boot <- c(quantile(boot.sample.mean, 0.025), quantile(boot.sample.mean, 0.975))
int.boot

     2.5%     97.5% 
0.2566667 0.3433333 

Comparing Two Distributions

To compare two distributions, the Kolmogorov-Smirnov (KS) statistics is a helpful measure. For example, say if we take another sample from the population data earlier, and compare it with the earlier sample we took.

ind.sample2 <- sample.int(n=N, size=n)
data.sample2 <- data.population[ind.sample2]
ks.test(data.sample, data.sample2)

p-value will be approximate in the presence of ties

    Two-sample Kolmogorov-Smirnov test

data:  data.sample and data.sample2
D = 0.043333, p-value = 0.9408
alternative hypothesis: two-sided

Now let’s take another sample that is completely different.

data.sample2 <- rbinom(n=n, size=1, prob=0.7)
ks.test(data.sample, data.sample2)

p-value will be approximate in the presence of ties

    Two-sample Kolmogorov-Smirnov test

data:  data.sample and data.sample2
D = 0.39, p-value < 2.2e-16
alternative hypothesis: two-sided

Important Features

Often we woule like to pick out important reasons resulting an outcome. For example, we would like to figure out what are the import reasons why students like/dislike video games.

We first load the data.

data <- read.table("videodata.txt", header=TRUE)

The package we are going to introduce today is named ‘tree’.

install.packages("tree")

Error in install.packages : Updating loaded packages

library(tree)

Since our task only ask for like or dislike, we need to congregate the column ‘like’ into 2 categories from 5, as a new column named ‘dis_like’. We will denote 1=never played, 4=not really, and 5=not at all of the column ‘like’ as 0=dislike in the new column ‘dis/like’, and the rest as 1=like.

data['dis_like'] <- rep(NA, dim(data)[1])
for(i in 1:dim(data)[1]){
  like <- data[i, 'like']
  if(like==0 || like==4 || like==5){
    data[i, 'dis_like'] = 0
  }else{
    data[i, 'dis_like'] = 1
  }
}

We then fit a tree on the data frame with appropriate columns.

data.tree <- tree(dis_like~educ+sex+age+home+math+work+own+cdrom+grade, data=data)
plot(data.tree, type="uniform")
text(data.tree)