# Function to count number of times C.I. captures # population proportion in a simulation binom.ci.answer<-function(data.x,p.pop){ x<-data.x n<-length(x) prop.obs<-mean(x) LL<-prop.obs-((prop.obs*(1-prop.obs)/n)^.5)*1.96 UL<-prop.obs+((prop.obs*(1-prop.obs)/n)^.5)*1.96 if (LL<=p.pop && UL>=p.pop) ci.answer<-1 else ci.answer<-0 return(ci.answer) } # Initialize population and monte carlo parameters p<-.20 n<-100 mc.samples<-2000 # Draw a single sample of size n binary.obs.vector<-rbinom(n, 1, p) prop.obs<-mean(binary.obs.vector) cat("Observed Proportion") cat("\n") prop.obs # Normal theory 95% C.I. for an observed proportion # with a given sample size n cat("Normal theory confidence interval for Observed Proportion") cat("\n") prop.obs-((prop.obs*(1-prop.obs)/n)^.5)*1.96 prop.obs+((prop.obs*(1-prop.obs)/n)^.5)*1.96 # Monte Carlo simulation of binomial sampling distribution of sample size n x.mat<-matrix(rbinom(n*mc.samples,1,p),ncol=n) x.mean<-apply(x.mat, 1, mean) x.ci<-apply(x.mat, 1, binom.ci.answer,p.pop=p) hist(x.mean) prop.captured<-mean(x.ci)*100 cat("Percentage that true Population Value is captured with Monte Carlo Sim.") cat("\n") prop.captured cat("Average Percentage over all Monte Carlo Sim.") cat("\n") mean(x.mean) cat("Confidence intervals calculated from Monte Carlo Sim.") cat("\n") quantile(x.mean, probs=c(.025, .975)) # Margin of error for a specified p with 95% C.I. p<-.50 margin<-((p*(1-p)/n)^.5)*1.96 cat("Margin of error for a specified p with 95% C.I.") cat("\n") cat("margin of error\n") margin cat("for p") cat("\n") p # Sample size needed for a specified p and a specified # margin of error with 95% C.I. p<-.50 margin<-.01 n<-(p*(1-p))*((1.96/margin)^2) cat("Sample size needed for a specified p and a specified margin of error with 95% C.I.") cat("\n") cat("for p") cat("\n") p cat("with margin of error") cat("\n") margin cat("one needs a sample size of") cat("\n") n # Bootstrap simulation p<-.20 bootstrap.samples<-2000 n<-100 x.mat<-matrix(sample(binary.obs.vector, n*bootstrap.samples,replace=T),ncol=n) x.mean<-apply(x.mat, 1, mean, trim=.20) x.ci<-apply(x.mat, 1, binom.ci.answer,p.pop=p) hist(x.mean) prop.captured<-mean(x.ci)*100 cat("Percentage of true Pop. values captured by Bootstrap Sim.") cat("\n") cat("with bootstrap samples\n") bootstrap.samples cat("\n") cat("proportion captured") cat("\n") prop.captured cat("Average Bootstrap proportion") cat("\n") mean(x.mean) cat("Bootstrap C.I.") cat("\n") quantile(x.mean, probs=c(.025, .975))