1Stat 511 HW#7 Spring 2009 (not to be collected) 1. An R function that will create bootstrap samples (taken from page 399 of Efron and Tibshirani) is > bootstrap<-function(x,nboot,theta) + {data<-matrix(sample(x,size=length(x)*nboot,replace=T),nrow=nboot) + return(apply(data,1,theta)) + } (A fancier version that will compute BCa and ABC intervals can be had by downloading the R contributed package "bootstrap" that has been written to accompany the Efron and Tibshirani book.) Load the MASS package. Below are times to failure (in millions of cycles) for high speed turbine engine bearings made from two different compounds (taken from "Analysis of Single classification Experiments Based on Censored Samples from the Two-Parameter Weibull Distribution: by J.I. McCool, Journal of Statistical Planning and Inference, 1979). (These are pretty small samples to be trying out a bootstrap technology, but for sake of exercise we will ignore this fact.) Compound 1 3.03 5.53 5.60 9.30 9.92 12.51 12.95 15.21 16.04 16.84 Compound 2 3.19 4.26 4.47 4.53 4.67 4.69 5.78 6.79 9.37 12.75 Enter these two sets of values into R as > compound1<-c(3.03,5.53,5.60,9.30,9.92,12.51,12.95,15.21,16.04,16.84) > compound2<-c(3.19,4.26,4.47,4.53,4.67,4.69,5.78,6.79,9.37,12.75) a) Make normal plots for the two samples as > qqnorm(compound1) > qqnorm(compound2) Do you expect "constant variance/normal distribution" ordinary statistical methods to be reliable in the analysis of these data? Why? b) Make normal plots of the log-lifetimes here as > qqnorm(log(compound1)) > qqnorm(log(compound2)) Are you any more optimistic about the reliability of "constant variance/normal distribution" ordinary statistical methods in the analysis of the log lifetimes? c) Consider inference for the median of a distribution F that one might assume to be generating the transplant Compound 2 lifetimes. The sample median Compound 2 lifetime can be had by typing > median(compound2) As a way of assessing the precision of the sample median (as an estimator of the median of F , you may generate, say, 10,000B = bootstrap sample medians by typing2> B<-10000 > comp2boot.non<-bootstrap(compound2,B,"median") Based on these values, what number might you report as a standard error for the sample median? b) Use the bootstrapped sample medians to produce a 95% (unadjusted) percentile bootstrap confidence interval for the median of F . Some R code for doing this is > kl<-floor((B+1)*.025) > ku<-B+1-kl > sortcomp2boot.non<-sort(comp2boot.non) > sortcomp2boot.non[kl] > sortcomp2boot.non[ku] c) A parametric bootstrap approach here might be to use a Weibull distribution to describe transplant survival times. It's easy to use maximum likelihood to fit a Weibull to these data by typing > fit2<-fitdistr(compound2,"weibull") > fit2 Then a set of 10,000 simulated sample medians from the fitted distribution can be simulated by typing > Wboot<-function(samp,nboot,theta,shape,scale) + {data<-scale*matrix(rweibull(samp*nboot,shape),nrow=nboot) + return(apply(data,1,theta)) + } > comp2boot.Wei<-Wboot(10,B,"median",fit2$estimate[1],fit2$estimate[2]) Use these to produce a parametric bootstrap standard error for the sample median and a parametric bootstrap 95% (unadjusted) percentile bootstrap confidence interval for the median of F . How do these compare to what you found above? d) As Koheler's Assignment #9 from Spring 2002 points out, large sample theory says that iffis the density function forF and()1.5Fθ−= (the population median) the standard deviation of the sample median is approximately ()12fnθ How does a (Weibull) parametric bootstrap standard error for the sample median Compound 2 lifetime compare to a plug-in version of the above (where the sample median is used in place of θ and the fitted Weibull density is used forf)? e) If one models the Compound 1 survival times as iid from second distribution,G , it might be of interest to estimate the difference in medians()()11.5 .5FG−−−. Modify the above to create 10,000B = (nonparametric) bootstrap sample medians from the Compound 1 cases. Then compute 10,000 differences between a bootstrap Compound 2 sample median and a bootstrap Compound 1 sample median. Based on these 10,000 differences, find 95% percentile confidence limits for the difference in underlying median lifetimes. Is this difference clearly non-zero?32. Problem 3 of the Spring 2003 Stat 511 Final Exam is a bootstrap problem. Do the R computing required to get your own copy of the printout provided on that Exam. Then answer the Exam questions. 3. Problem 3 of the Spring 2004 Stat 511 Final Exam is a bootstrap problem. Do the R computing required to get your own copy of the printout provided on that Exam. Then answer the Exam questions. 4. The table below contains information on 23 (out of 24) pre-Challenger space shuttle flights. (On one flight, the solid rocket motors were lost at sea and so no data are available.) Provided are launch temperatures,t (in ° F), and a 0-1 response,y, indicating whether there was post-launch evidence of a field joint primary O-ring incident. (O-ring failure was apparently responsible for the tragedy.) 1y= indicates that at least one of the 6 primary O-rings showed evidence of erosion. Temperature O-ring Incident? Temperature O-ring Incident?66 0 67 0 70 1 53 1 69 0 67 0 68 0 75 0 67 0 70 0 72 0 81 0 73 0 76 0 70 0 79 0 57 1 75 1 63 1 76 0 70 1 58 1 78 0 Treat the response variables,iy, as Bernoulli distributed (binomial()1,ip ) and independent launch to launch. Note that Eypμ==here. We'll model ()01log log11pnhtpnμμββμ⎛⎞⎜⎟⎛⎞===+⎜⎟⎜⎟−⎝⎠⎜⎟−⎝⎠ This is a generalized linear model, with binomial (Bernoulli) response, and the "logit" link … otherwise known as a logistic regression model. We may use R's glm function to analyze these data. (You may learn about the function in the usual way, by typing ?glm .) Load the MASS package. Read in the data for this problem by typing > temp<-c(66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81, 76,79,75,76,58) > incidents<-c(0,1,0,0,0,0,0,0,1,1,1,0,0,1,0,0,0,0,0,0,1,0,1) And create and view a logical variable equivalent to the 0-1 response variable by typing > indicate<-(incidents>0) > indicate a) We may then fit and summarize a generalized linear model here by typing4 >