Chapter 5 and 6 DJM 21 February 2017 Exam admonitions 0. You must compile the .Rmd to .html 1. No code in the html. 2. If you put plots or tables in, you must talk about them. Rule of thumb: if you don’t have anything good to say about a number, don’t give the number (or plot) at all. 3. MOST IMPORTANT you must explain your results. Simply providing them is not likely to get you credit. 4. Look over the model solutions for the last assignments.

Project progress report • Due 8 March at 11:59 pm ( just over 2 weeks from today) • Your report should have 3 components: 1. A list of teammate names and an explanation of what is interesting to you about this data. 2. A short introductory paragraph introducing the data and describing some potential questions you might investigate. 3. A lengthy exploratory data analysis. • The third part is a big deal. • You need to provide evidence that you have explored the data carefully and meaningfully. • Code must be integrated. • Think of this like HW 2. • Just like with HW 2 and HW 3, much of what you do on the midterm will end up in the final report, so spend the time to do a good job.

Simulation Why Simulation? • Up until now, when we do linear models, we used t-statistics, p-values, CIs ˆ if the model is true and we • These things are based on the sampling distribution of the estimators (β) don’t do any model selection. • What if we do model selection, use Kernels, think the model is wrong? • None of those formulas work. And analogous formulas can be impossible (or painfully annoying) to derive.

1

Some simulation basics set.seed(2018-02-20) sample(1:10, replace=TRUE, prob=1:10/10) ##



6

8 10

9

7

6

3

9 10

8

sample(letters[1:10], replace=TRUE, prob=1:10/10) ##

 "i" "i" "a" "h" "e" "e" "h" "f" "d" "j"

sample(letters[1:10], replace=TRUE) ##

 "c" "e" "h" "j" "i" "b" "e" "a" "a" "j"

sample(letters[1:10]) ##

 "g" "i" "b" "d" "c" "h" "e" "j" "a" "f"

Resampling data set.seed(2018-02-20) n = 100; x = runif(n) df = data.frame(x=x, y=3+2*x+rnorm(n)) ggplot(df, aes(x,y)) + geom_point(color=blue,shape=17)

2

6

5

y

4

3

2

1 0.00

0.25

0.50

0.75

x

A sample (with replacement), and a new draw from the same distribution resample <- function(df) { stopifnot(is.data.frame(df) | is.matrix(df)) df[sample(1:nrow(df), replace = TRUE),] } df2 = resample(df) xn = runif(n) df3 = data.frame(x=xn, y=3+2*xn+rnorm(n)) df = rbind(df,df2,df3) df\$grp = rep(c('original','resample','new draw'), each=n) p <- ggplot(df, aes(x,y,color=grp)) + geom_point(aes(shape=grp)) + scale_color_manual(values=c(red,blue,orange)) + theme(legend.title = element_blank(),legend.position = 'bottom') p

3

1.00

y

6

4

2

0.00

0.25

0.50

0.75

x new draw

original

Add some lines p + geom_smooth(method='lm', se = FALSE)

4

resample

1.00

y

6

4

2

0.00

0.25

0.50

0.75

x new draw

original

resample

Using simulations to check modeling assumptions x = runif(n) - 0.5; y = 3+2*x + rnorm(n)*x^2 dfHetero = data.frame(x=x, y=y) ggplot(dfHetero, aes(x,y)) + geom_point(color=blue,shape=17) + geom_smooth(method='lm', se=FALSE,color=blue) + geom_abline(intercept = 3, slope=2, color=green)

5

1.00

4.0

y

3.5

3.0

2.5

2.0

−0.50

−0.25

0.00

0.25

x

If the noise is homoskedastic. . . • The red and blue points should have the same distribution heteromod = lm(y~x,data=dfHetero) dfHetero\$resids = residuals(heteromod) dfHetero\$resample = sample(residuals(heteromod), replace = TRUE) dfHetero %>% gather(key='type', value='resids',-c(y,x)) %>% ggplot(aes(x,resids,color=type,shape=type)) + geom_point() + scale_color_manual(values=c(red,blue)) + theme(legend.title = element_blank(),legend.position = 'bottom') + geom_hline(yintercept=0, color=blue)

6

0.50

resids

0.2

0.0

−0.2

−0.4 −0.50

−0.25

0.00

0.25

x resample

resids

That one was easy x = runif(n)-0.5 y = 3+2*x + c(arima.sim(list(ar=.8), n, rand.gen = function(n) 0.1* rt(n, df=5))) dfTS = data.frame(x=x, y=y) ggplot(dfTS, aes(x,y)) + geom_point(color=blue) + geom_smooth(method='lm',se=FALSE, color=red) + geom_abline(intercept = 3, slope = 2, color=blue)

7

0.50

4.0

y

3.5

3.0

2.5

2.0

−0.50

−0.25

0.00

0.25

x

If the noise is homoskedastic. . . • The red and blue points should have the same distribution tsMod = lm(y~x, data=dfTS) dfTS\$resids = residuals(tsMod) dfTS\$resample = sample(residuals(tsMod), replace = TRUE) dfTS %>% gather(key='type', value='resids',-c(y,x)) %>% ggplot(aes(x,resids,color=type,shape=type)) + geom_point() + scale_color_manual(values=c(red,blue)) + theme(legend.title = element_blank(),legend.position = 'bottom') + geom_hline(yintercept=0, color=blue)

8

0.50

0.25

resids

0.00

−0.25

−0.50

−0.50

−0.25

0.00

0.25

x resample

resids

But. . . lag.resids = with(dfTS, data.frame(lag.resids = resids[-n], resids = resids[-1])) lag.resids\$resample = sample(lag.resids\$resids, replace = TRUE) lag.resids %>% gather(key='type', value='resids',-lag.resids) %>% ggplot(aes(lag.resids,resids,color=type,shape=type)) + geom_point() + scale_color_manual(values=c(red,blue)) + theme(legend.title = element_blank(),legend.position = 'bottom') + geom_smooth(method='lm',se=FALSE)

9

0.50

0.25

resids

0.00

−0.25

−0.50

−0.50

−0.25

0.00

lag.resids resample

Another useful command sample.int(10) ##

 10

5

8

1

4

9

7

3

2

6

10

resids

0.25

What’s the deal with this Bootstrap?

What’s the deal with this Bootstrap? • Suppose I want to estimate something and get a CI. • But I don’t know how to calculate the CI (or maybe I do, but it’s hard) • Then what?

Example 1 • Let Xi ∼ χ24 . ¯ then by the CLT (if n is big), • I know if I estimate the mean with X, √

¯ − E [X]) n(X ≈ N (0, 1). s 11

• This gives me a 95% confidence interval like √ ¯ ± 2 ∗ s/ n X • But I don’t want to estimate the mean, I want to estimate the median. ggplot(data.frame(x=c(0,12)), aes(x)) + stat_function(fun=function(x) dchisq(x, df=4), color=blue) + geom_vline(xintercept = 4, color=blue) + # mean geom_vline(xintercept = qchisq(.5,4), color=red) # median

0.15

y

0.10

0.05

0.00 0.0

2.5

5.0

7.5

10.0

x

Now what • I give you a sample of size 50, you give me the sample median. • How do you get a CI? • You can use the bootstrap! set.seed(2018-02-20) x = rchisq(n, 4) (med = median(x)) ##  2.570231 B = 100 alpha = 0.05 bootMed <- function(x) median(sample(x, replace=TRUE)) bootDist = replicate(B, bootMed(x)) bootCI = 2* med - quantile(bootDist, probs = c(1-alpha/2, alpha/2)) ggplot(data.frame(bootDist), aes(bootDist)) + geom_density(color=blue) + geom_vline(xintercept = bootCI, col=blue, linetype=2) + 12

12.5

geom_vline(xintercept = med, col=blue) + geom_vline(xintercept = qchisq(.5, 4), col=red) # truth

1.5

density

1.0

0.5

0.0 2.0

2.5

3.0

3.5

bootDist

An alternative • In that bootstrap, I didn’t use any information about the data-generating process. • What if I told you that the data came from a χ2 , but I didn’t tell you the degrees of freedom? • You could try a “parametric” bootstrap: xbar = mean(x) s = sd(x) ParaBootSamp <- function(B, xbar, s){ means = rnorm(B, mean=xbar, sd=s/sqrt(n)) meds = qchisq(.5, means) return(meds) } ParaBootDist = ParaBootSamp(B, xbar, s) ParaBootCI = 2* med - quantile(ParaBootDist, probs = c(1-alpha/2, alpha/2)) ggplot(data.frame(bootDist), aes(ParaBootDist)) + geom_density(color=blue) + geom_vline(xintercept = ParaBootCI, col=blue, linetype=2) + geom_vline(xintercept = med, col=blue) + geom_vline(xintercept = qchisq(.5, 4), col=red) # truth

13

density

0.75

0.50

0.25

0.00 2

3

ParaBootDist

In truth • Let’s compare these intervals • The nonparametric bootstrap (first one) had a width of bootCI - bootCI ## 2.5% ## 1.156143 • The parametric bootstrap (second one) had a width of ParaBootCI - ParaBootCI ## 2.5% ## 1.479937 • Using theory, we could find the exact CI. In this case, it has a width of 1.76.

14

4

Bootstrap diagram

15

Bootstrap intuition

Bootstrap error sources (From the bottom up on the last slide) 1. Simulation error: using only B samples to estimate F with Fˆ . 2. Statistical error: our data depended on a sample from the population. We don’t have the whole population so we make an error by using a sample (Note: this part is what always happens with data, and what the science of statistics analyzes.) 3. Specification error: If we use the model based bootstrap, and our model is wrong, then we think we are badly overconfident in our assessment of error.

Recap • There are essentially 2 types of bootstrap 1. Parametric 2. Nonparametric • If you really believe your model, use the first • If not, use the second 16

• Both are valid

Example 2 (using code from Book on new data) ## ## Attaching package: 'MASS' ## The following object is masked from 'package:dplyr': ## ## select library(MASS) ggplot(fatcats, aes(Bwt, Hwt)) + geom_point(color=blue)

Hwt

16

12

8

2.0

2.5

3.0

Bwt

A model ## Running the model on the original data cats.lm = lm(Hwt ~ 0+Bwt,data=fatcats) summary(cats.lm) ## ## Call: 17

3.5

## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

lm(formula = Hwt ~ 0 + Bwt, data = fatcats) Residuals: Min 1Q Median -5.4713 -0.6530 -0.0427

3Q 0.8189

Max 4.9289

Coefficients: Estimate Std. Error t value Pr(>|t|) Bwt 3.88959 0.04385 88.7 <2e-16 *** --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.456 on 143 degrees of freedom Multiple R-squared: 0.9821, Adjusted R-squared: 0.982 F-statistic: 7868 on 1 and 143 DF, p-value: < 2.2e-16

confint(cats.lm) ## 2.5 % 97.5 % ## Bwt 3.802907 3.976266

I think that that CI is wrong. . . qqnorm(residuals(cats.lm)) qqline(residuals(cats.lm))

2 0 −2 −4

Sample Quantiles

4

Normal Q−Q Plot

−2

−1

0 Theoretical Quantiles

18

1

2

Functions to borrow resample <- function(x) { sample(x, replace = TRUE) } resample.data.frame <- function(data) { sample.rows <- resample(1:nrow(data)) return(data[sample.rows, ]) } rboot <- function(statistic, simulator, B) { tboots <- replicate(B, statistic(simulator())) if (is.null(dim(tboots))) { tboots <- array(tboots, dim = c(1, B)) } return(tboots) } bootstrap <- function(tboots, summarizer, ...) { summaries <- apply(tboots, 1, summarizer, ...) return(t(summaries)) } equitails <- function(x, alpha) { lower <- quantile(x, alpha/2) upper <- quantile(x, 1 - alpha/2) return(c(lower, upper)) } bootstrap.ci <- function(statistic = NULL, simulator = NULL, tboots = NULL, B = if (!is.null(tboots)) { ncol(tboots) }, t.hat, level) { if (is.null(tboots)) { stopifnot(!is.null(statistic)) stopifnot(!is.null(simulator)) stopifnot(!is.null(B)) tboots <- rboot(statistic, simulator, B) } alpha <- 1 - level intervals <- bootstrap(tboots, summarizer = equitails, alpha = alpha) upper <- t.hat + (t.hat - intervals[, 1]) lower <- t.hat + (t.hat - intervals[, 2]) CIs <- cbind(lower = lower, upper = upper) return(CIs) }

Model ## Simulator resamp.resids.cats <- function(){ resids = residuals(cats.lm) newResids = sample(resids, replace=TRUE) # resample the residuals from the original model newCats = data.frame(Bwt = fatcats\$Bwt, Hwt=fitted(cats.lm) + newResids) # create a new dataframe

19

# with the original x's but new y's return(newCats) } ## Estimator fitCats <- function(newCats) coef(lm(Hwt~0+Bwt, data=newCats)) # get the coef from OLS fitCats(fatcats) # test the above on original data, should give same coef ## Bwt ## 3.889586

Model-based bootstrap cisPara = bootstrap.ci(statistic = fitCats, simulator = resamp.resids.cats, B = 1000, t.hat = fitCats(fatcats), level = 0.95)

Nonparametric bootstrap resamp.cats <- function() resample.data.frame(fatcats) cisNonPara = bootstrap.ci(statistic = fitCats, simulator = resamp.cats, B = 1000, t.hat = fitCats(fatcats), level = 0.95) # use the prev func to # bootstrap on resampled data cisPara ## lower upper ## Bwt 3.80835 3.982101 cisNonPara ## lower upper ## Bwt 3.798535 3.975114

Bootstrapping with nonparametric regression • This is a bit harder • The reason is that we use CV to choose the bandwidth • So we have to repeat that step in the bootstrapping • That is: 1. 2. 3. 4. 5.

Input data Use CV to choose a smoothing parameter Use the chosen parameter to estimate the smooth function Resample the data Using this new data, repeat 2 and 3

20

## Chapter 5 and 6 - GitHub

Mar 8, 2018 - These things are based on the sampling distribution of the estimators (ËÎ²) if the model is true and we don't do any model selection. â¢ What if we do model selection, use Kernels, think the model is wrong? â¢ None of those formulas work. And analogous formulas can be impossible (or painfully annoying) to.

#### Recommend Documents

AIFFD Chapter 6 - Mortality - GitHub
6.5 Adjusting Catch-at-Age Data for Unequal Recruitment . . . . . . . . . . . . . . . . . . . . . . ...... System: Windows, i386-w64-mingw32/i386 (32-bit). Base Packages: base ...

AIFFD Chapter 5 - Age and Growth - GitHub
May 13, 2015 - The following additional packages are required to complete all of the examples (with ... R must be set to where these files are located on your computer. ...... If older or younger age-classes are not well represented in the ... as the

Chapter 5
not in the domain. The only critical point is x = 0. As x moves away from 0 on either side, the values of y decrease. The function has a local maximum value at (0, ...... (b) Since. ,. dV. dV dr dt dr dt. = we have. 2 . dV dr rh dt dt Ï. = (c). 2. 2

Chapter 4 - GitHub
The mathematics: A kernel is any function K such that for any u, K(u) â¥ 0, â« duK(u)=1 and â« uK(u)du = 0. â¢ The idea: a kernel is a nice way to take weighted averages. The kernel function gives the .... The âbig-Ohâ notation means we have

Chapter 2 - GitHub
Jan 30, 2018 - More intuitively, this notation means that the remainder (all the higher order terms) are about the size of the distance between ... We don't know Âµ, so we try to use the data (the Zi's) to estimate it. â¢ I propose 3 ... Asymptotica

Chapter 3 - GitHub
N(0, 1). The CLT tells us about the shape of the âpilingâ, when appropriately normalized. Evaluation. Once I choose some way to âlearnâ a statistical model, I need to decide if I'm doing a good job. How do I decide if I'm doing anything good?

chapter 5-6 Mr Popper's Penguins.pdf
Summer 2011. Chapter 5- Troubles with a Penguin. Language Arts. Vocabulary â. Look up the following words and write the definitions. icebox â ...

Chapter 1 - GitHub
Jan 23, 2017 - 1. What are all these things? 2. What is the mean of yi? 3. What is the distribution of Ïµi? 4. What is the notation X or Y ? Drawing a sample yi = xi Î² + Ïµi. Write code which draws a sample form the population given by this model. p

Chapter 6 - Home
write each new term and its definition in your notebook. â¡ communication (p. 174). â¡ verbal communication. (p. 175). â¡ nonverbal communication ..... Page 9 ...

Chapter Number 6 Chapter Title Incentives and ...
Mar 30, 2015 - well as by technology. .....  Casper, B. A., and S. A. Tyson. .... cisions under incomplete information, whereas we abstract from collective action ...

Chapter 6
when a lot of features are irrelevant and drown out the relevant features' signal in the distance calculations. Notice that the nearest-neighbor method can easily be applied to regression problems with a real-valued target variable. In fact, the meth

Chapter 5 - DLSCRIB
Three different washing solutions are being compared to study their ... Plot the mean tensile strengths observed for each chemical type in Problem 4.3 and ...... np y p y .... h... n-1. Treatment x Squares. Squares. Treatments .... h.j.. SS. SS np y

Chapter 5
Every Document object has: â¢forms - an array of references to the forms of the document. â¢Each forms object has an elements array, which has references to the form's elements. Document also has property arrays for anchors, links, & images. JavaSc

AIFFD Chapter 12 - Bioenergetics - GitHub
The authors fit a power function to the maximum consumption versus weight variables for the 22.4 and ... The linear model for the 6.9 group is then fit with lm() using a formula of the form ..... PhD thesis, University of Maryland, College Park. 10.

AIFFD Chapter 10 - Condition - GitHub
May 14, 2015 - 32. This document contains R versions of the boxed examples from Chapter 10 of the âAnalysis and Interpretation of Freshwater Fisheries Dataâ ...

chapter iv: the adventure - GitHub
referee has made changes in the rules and/or tables, simply make a note of the changes in pencil (you never kno, ,hen the rules ,ill ... the Game Host's rulebook until they take on the mantle of GH. The excitement and mystery of ...... onto the chara

CHAPTER 5: Graphs and Trees - DAINF
Page 166 Mathematical Structures for Computer Science Gersting. CHAPTER ... Not isomorphic; graph in (b) has a node of degree 5, graph in (a) does not. 14. f:.

Chapter 5 Lab 5-1, Configure and Verify Path Control
Note: This lab uses Cisco 1841 routers with Cisco IOS Release 12.4(24)T1, ... Cisco IOS Software versions if they have comparable capabilities and features.

AIFFD Chapter 4 - Recruitment - GitHub
Some sections build on descriptions from previous sections, so each ... setwd("c:/aaaWork/web/fishR/BookVignettes/AIFFD/") ... fact is best illustrated with a two-way frequency table constructed from the two group factor variables with ..... 10. Year

CHAPTER 6 Clean.pdf
time by approved providers. Approved providers include: (A) The American Physical Therapy Association (APTA),. including any sections, credentialed ...

Chapter 6 Friction.pdf
R. T1. R. A 1g. 1a. R. R. 1g. 1a.. 0.5g. 0.5g. a. A B. =0.2. 0.5kg. 1kg 1kg. =0.2. 2Ã0.5 16N. R1. 4g. 16N=T. 2R. 4Ã0.5. 30Â°. a. 1. 0.5 m/s2. 2kg. 4kg. 2. Page 3 of 10. Chapter 6 Friction.pdf. Chapter 6 Friction.pdf. Open. Extract. Open with. Sign

Chapter 6 Review.pdf
The Worthington Pottery Company manufactures beer mugs in batches of 120. and the overall rate of defects is 5%. Find the probability of having more than 6. defects in a batch. 8. A bank's loan officer rates applicants for credit. The ratings are nor

chapter 6.pdf
There was a problem previewing this document. Retrying... Download. Connect more apps... Try one of the apps below to open or edit this item. chapter 6.pdf.