Non-Parametric Statistics

Instructions

This exam is due May 4, 2019 at 11:55pm. There are 100 points available for students in both 388/488. Students enrolled in 388 will be scored out of 80 points. Don’t cheat. Seriously, don’t cheat. Have a great summer. Try not to do anything stupid.

(1) [10 points]

What is the difference between parametric and non-parametric statistics (2-3 lines should suffice for an answer)?

Parametric statistics makes assumptions about the probability distributions of data. These assumptions often involve parametric distributions or the central limit theorem is applied as samples get larger. Non-parametric statistics, however, does not make these distributional assumptions about the data; usually non-parametric statistical methods make other, more flexible assumptions about the data. In general, assuming more about data corresponds with greater power in statistical testing, but not always. Some advantages of non-parametric statistical methods over parametric statistical methods include the following:

• Fewer/less stringent assumptions
• Exactness of p-values, coverage probabilities & error rates
• Typically easier to apply than their normal-theory counterparts
• Procedures are typically easier to understand
• Procedures are only slightly less efficient than their normal theory competitors when underlying populations are normal, and they are at least mildly, if not wildly, more efficient when underlying populations are not normal
• Methods are relatively insensitve to outliers
• Widely applicable since many procedures simply rely on the ranks of observations rather than their true magnitudes

(Taken directly from class notes)

(2) [10 points]

Give an example of a hypothesis test that can be tested using both parametric and non-parametric tests. State the null and alternative hypotheses and describe (or give the names of) the parametric and nonparametric tests used to test the hypothesis that you have given as an example. (HINT: Pick an example that is easy for yourself!)

Randomly generated data

set.seed(2634)
x <- rnorm(100, mean = 0, sd = 1)
kable(head(x), caption = "100 Randomly Generated Variables from Standard Normal")

Parametric hypothesis test: One Sample t-Test
\(H_{0}: \mu = 0\)
\(H_{A}: \mu \neq 0\)

t.test(x)

    One Sample t-test

data:  x
t = -0.037259, df = 99, p-value = 0.9704
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.2193076  0.2112232
sample estimates:
   mean of x 
-0.004042188 

Assumptions:
* Data are IID * Approximately normally distributed * No outliers

Non-parametric alternative test: Sign Test (let \(\theta\) be the median)
\(H_{0}: \theta = 0\)
\(H_{A}: \theta \neq 0\)

table(sign(x))
binom.test(x = 47, n = 100, p = 0.5)

-1  1 
47 53 

    Exact binomial test

data:  47 and 100
number of successes = 47, number of trials = 100, p-value = 0.6173
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.3694052 0.5724185
sample estimates:
probability of success 
                  0.47 

Assumptions:
* data are IID

(3) [10 points]

A random forest consists of many trees, but there is a difference between the trees in a random forest and a regular tree model. Describe the differences between these trees.

The major differences between trees in CART models and Random Forests is the way trees are “grown.” When trees are being modeled in CART, all observations start in a single node and all binary splits for all independent variables are searched until the one that reduces within leaf variance the most is found. This process is repeated until the reduction in variance meets some threshold or a minimum number of observations in each node is reached.

Random forests are an ensemble method combining CART and bootstrapping. When a random forest model is being trained, a bootstrap sample is taken from a random sample. At each node a random subset of the independent variables is taken, all possible binary splits are considered only for that subset, and the split that reduces variance most is selected. Since each tree is constructed using only a sample of the data, there is an inherent holdout dataset; predictions can be made on the data which was not used to train the model and predictions are averaged across trees.

In summary, the major difference between trees in CART models and random forests is the number of independent variables considered at each split and the sample of observations that the models are trained on.

(4) Answer the following question with the “USArrest” dataset in R.

rm(list = ls())
data("USArrests")
df <- USArrests
kable(head(df), caption = "USA Arrests Data from R")

a. [10 points] Using the “USArrests” data set that is pre-loaded in R, construct a kernel density estimator of the number of assault arrests and show the density in a figure. (hint: ggplot2 makes this very easy).

set.seed(2634)

ggplot(df, aes(x = Assault, y = ..density..)) + geom_histogram(bins = 10, closed = "right") + 
    geom_density(bw = 5, col = "red") + labs(title = "Kernel Density Estimator", 
    subtitle = "Bandwidth = 55", y = "Density", x = "Number of Assault Arrests")

ggplot(df, aes(x = Assault, y = ..density..)) + geom_histogram(bins = 10, closed = "right") + 
    geom_density(bw = 10, col = "red") + labs(title = "Kernel Density Estimator", 
    subtitle = "Bandwidth = 10", y = "Density", x = "Number of Assault Arrests")

ggplot(df, aes(x = Assault, y = ..density..)) + geom_histogram(bins = 10, closed = "right") + 
    geom_density(bw = 15, col = "red") + labs(title = "Kernel Density Estimator", 
    subtitle = "Bandwidth = 15", y = "Density", x = "Number of Assault Arrests")

b. [10 points] Using the “USArrests” data set that is pre-loaded in R, fit a loess regression model with UrbanPop as a predictor and Assault as the response variable. (hint: ggplot2 makes this very easy).

# Fitting Loess regression models with spans between .1 and 1
fits <- data_frame(span = seq(from = 0.1, to = 2, by = 0.1)) %>% group_by(span) %>% 
    do(augment(loess(Assault ~ UrbanPop, df, degree = 2, span = .$span)))

# Plotting smoothing parameters
ggplot(fits, aes(x = UrbanPop, y = Assault, frame = span)) + geom_point() + 
    geom_line(aes(y = .fitted), color = "red") + transition_states(span, wrap = TRUE, 
    transition_length = 3) + labs(title = "Loess Regression: Assault Arrests as a function of Population", 
    subtitle = "Span: {closest_state}", x = "Urban Population", y = "Assault Arrests")

(5) Answer the following questions with the data sets in the dslabs pacakge. We will focus on measles. (I did some data cleaning for you. Also, assume that the 0’s in the early years are correct.)

rm(list = ls())
measles <- subset(us_contagious_diseases, disease == "Measles")
measles <- melt(measles, id = c("disease", "state", "year"), measure = "count")
measles <- dcast(measles, formula = year ~ state)
# Remove the space in some of the state names.
names(measles)[10] <- "DC"
names(measles)[31] <- "NH"
names(measles)[32] <- "NJ"
names(measles)[33] <- "NM"
names(measles)[34] <- "NY"
names(measles)[35] <- "NC"
names(measles)[36] <- "ND"
names(measles)[41] <- "RI"
names(measles)[42] <- "SC"
names(measles)[43] <- "SD"
names(measles)[50] <- "WV"
kable(head(measles), caption = "Measles Data")

a. [10 points] Construct a CART model to predict number of measles cases in Illinois using the number of cases in other states as potential predictors. Display the best tree that you find and try to offer some interpretation.

# Set seed and partition training and test data
set.seed(2634)
ind <- sample(1:dim(measles)[1], 50, replace = FALSE)
df.train <- measles[ind, ]
df.test <- measles[-ind, ]


m3 <- rpart(formula = Illinois ~ ., data = df.train, method = "class", control = list(minsplit = 10, 
    maxdepth = 12, xval = 10))


hyper_grid <- expand.grid(minsplit = seq(5, 20, 1), maxdepth = seq(8, 15, 1))


models <- list()
for (i in 1:nrow(hyper_grid)) {
    # get minsplit, maxdepth values at row i
    minsplit <- hyper_grid$minsplit[i]
    maxdepth <- hyper_grid$maxdepth[i]
    # train a model and store in the list
    models[[i]] <- rpart(formula = Illinois ~ ., data = df.train, method = "class", 
        control = list(minsplit = minsplit, maxdepth = maxdepth))
}

get_cp <- function(x) {
    min <- which.min(x$cptable[, "xerror"])
    cp <- x$cptable[min, "CP"]
}

# function to get minimum error
get_min_error <- function(x) {
    min <- which.min(x$cptable[, "xerror"])
    xerror <- x$cptable[min, "xerror"]
}

hyper_grid %>% mutate(cp = purrr::map_dbl(models, get_cp), error = purrr::map_dbl(models, 
    get_min_error)) %>% arrange(error) %>% top_n(-5, wt = error)


optimal_tree <- rpart(formula = Illinois ~ ., data = df.train, method = "class", 
    control = list(minsplit = 13, maxdepth = 11, cp = 0.01))

pred <- predict(optimal_tree, newdata = df.test, type = "prob")
RMSE(pred = pred, obs = df.test$Illinois)
pred.v <- predict(optimal_tree, newdata = df.test, type = "vector")
RMSE(pred = pred.v, obs = df.test$Illinois)

# MSE
mse <- mean((pred - df.train$Illinois)^2)
mse
mse.v <- mean((pred.v - df.train$Illinois)^2)
mse.v

rpart.plot(optimal_tree)

  minsplit maxdepth         cp    error
1       18        9 0.02083333 1.020833
2       16       10 0.02083333 1.020833
3       20       10 0.02083333 1.020833
4        5       11 0.02083333 1.020833
5        8       14 0.02083333 1.020833
[1] 20994.23
[1] 20976.59
[1] 367342596
[1] 371806390

b. [10 points] Construct a random forest model to predict number of measles cases in Illinois using the number of cases in other states as potential predictors. Is the random forest model better than the CART model at predicting measles cases?

set.seed(2634)
IL.allForests <- randomForest(Illinois ~ ., data = df.train, importance = TRUE, 
    ntree = 1000)
plot(IL.allForests)

varImpPlot(IL.allForests)

kable(IL.allForests$importance)

# IL.Forests <-randomForest(Illinois ~., data = df.train, importance = TRUE,
# ntree = 1000)
predicted_IL_all <- predict(IL.allForests, newdata = df.test, OOB = TRUE)
mean((predicted_IL_all - df.train$Illinois)^2)

[1] 475351741

c. [10 points] Compute Kendall’s \(\tau\) between the Illinois and Massachusetts measles cases data. Compute the standard error of this estimator using the bootstrap procedure. Also, estimate the bias of this estimator.

d. [10 points] Perform PCA on the the measles data with states as variables. (Remember to exclude year as a predictor) How many components are needed to account for 85% of the variability? Also, try to interpret the first two principal components.

e. [10 points] Perform hierarchical clustering on the the measles data with states as observations. (Remember to exclude year data). How many clusters do there appear to be? Try to interpret each of the clusters.