Introduction to Simulations in R

1 Introduction to Simulations in RCharles DiMaggio, PhD, MPH, PA-CNew york university department of Surgery and population HealthNYU-Bellevue Division of Trauma and Surgical Critical CareJune 10, DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20151 / 48 Outline1sampling in R2simulating risk ratios3simulation for statistical inference4simulation to summarize and predict regression resultssimulating predictive uncertainty in complex models5simulation for model checking and fitPoisson exampleCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20152 / 48 This material has been shamelessly and read this book!

2 Gellman and Hill, Data Analysis Using Regression andMulitlevel/Hierarchical Models , Cambridge university Press, 2007.(mostly chapters 7 and 8). ~gelman/arm/Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20153 / 48sampling in ROutline1sampling in R2simulating risk ratios3simulation for statistical inference4simulation to summarize and predict regression resultssimulating predictive uncertainty in complex models5simulation for model checking and fitPoisson exampleCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20154 / 48sampling in Rsample()simple random samplesample(c("H","T"), size = 8, replace = TRUE) # fair coinsample(1.)

3 6, size = 2, replace = TRUE, prob=c(3,3,3,4,4,4)) #loaded dicereplace=TRUEto over ride the default sample without replacementprob=to sample elements with different probabilities, oversample based on some ()function allow you to make a reproducible set ofrandom DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20155 / 48sampling in Rprobability distributions in Rbeta(shape1, shape2, ncp)binom(size, prob)chisq(df, ncp)exp(rate)gamma(shape, scale)logis(location, scale)norm(mean, sd)pois(lambda)t(df, ncp)unif(min, max)Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20156 / 48sampling in Rconvention for using probability functions in Rdxxx(x,)returns the density or the value on the y-axis of a probabilitydistribution for a discrete value of xpxxx(q,)returns the cumulative density function (CDF) or the areaunder the curve to the left of an x value on a probability distributioncurveqxxx(p,)returns the quantile value, the standardized z value for xrxxx(n,) returns a random simulation of size nqnorm( )qnorm( )Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20157 / 48sampling in Rsampling from probability distributionsrnorm(6) # 6 std nrml distribution valuesrnorm(10, mean = 50, sd = 19) # set parametersrunif(n = 10, min = 0, max = 1)

4 #uniform distributionrpois(n = 10, lambda = 15) # Poisson distribution# toss coin 8 times using binomial distributionrbinom(n = 8, size = 1, p = )rbinom(8,1,.5) # args correct order# 18 trials, sample size 10, prob success =.2 rbinom(18, 10, )Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20158 / 48sampling in RExercise 1: Sampling and SimulationsCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 20159 / 48simulating risk ratiosOutline1sampling in R2simulating risk ratios3simulation for statistical inference4simulation to summarize and predict regression resultssimulating predictive uncertainty in complex models5simulation for model checking and fitPoisson exampleCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201510 / 48simulating risk ratiosbootstrapping a relative riskUse Simulations to approximate results when no direct or Risk ()to simulate (many times)

5 Rates of disease in exposed andunexposed populationsDivide results by the number of Simulations and use the mean tails for the point estimate and confidence DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201511 / 48simulating risk ratiosapproach to bootstrapping a relative risk1simulate 5000 replicate bernoulli trials in sample sizen1= exposed2divide those results byn1to get 5000 simulated risk estimates for theexposed group3repeat that process for the unexposed groupn24divide 5000 simulated risks in exposed by 5000 simulated risks inunexposed to get 5000 simulate relative risks5calcualate mean and tails from that populationCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201512 / 48simulating risk ratiosExample.

6 ASA and MIRR calculation usingepitab()Hennekens, 1987 study protective benefits of myocardial infarctions (fatal and non-fatal) among 11,037 people inthe treatment group189 MI s among 11,034 people in the placebo groupCalculate RR and CI using with log-approximated approach(epitools::epitab())library(epit ools) <- matrix(c(104,11037,189,11034),2,2)epitab ( , method="riskratio")Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201513 / 48simulating risk ratiossimulate RR estimateuserbinom()to repeat 5,000 times an experiment where we count thenumber of outcomes (MI s) in two populationsprobability of the outcome in a population defined by the results of theHennekens studyfor each replicate, divide the number outcomes by number of people ineach population to get 5,000 risk estimates for each group (treatmentand placebo)calculate the RR for each simulationcollect and describe (151)tx <- rbinom(5000, 11037, 104/11037)plac <- rbinom(5000, 11034, 189/11034) <- <- <- ( )quantile( , c( , ))sd( )

7 Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201514 / 48simulating risk ratiosTry writing a function to calculate bootstrap estimates of relative the function using the aspirin DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201515 / 48simulation for statistical inferenceOutline1sampling in R2simulating risk ratios3simulation for statistical inference4simulation to summarize and predict regression resultssimulating predictive uncertainty in complex models5simulation for model checking and fitPoisson exampleCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201516 / 48simulation for statistical inferencesome simple simulationsbirth genderpredict number of girls in 400 births, where prob of female birth <-rbinom(1,400.)

8 488) get distribution of the Simulations , repeat the simulation many < <-rbinom( , 400, .488)hist( )can do same thing with a loopvectorized operation preferred in R, but loops useful in < <-rep(NA, ) # create vector to store simulationsfor (i in 1 ){ [i]<-rbinom(1,400, )}hist( )Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201517 / 48simulation for statistical inferencemore complex simulationsaccount for twins1/125 chance fraternal twins, each with chance being girl1/300 chance identical twins, chance of both being <- sample (c("fraternal twin","identical twin","single birth"),size=400, replace=TRUE, prob=c(1/125, 1/300, 1 - 1/125 - 1/300))girls <- rep (NA, 400)for (i in 1:400){if ( [i]=="single birth"){girls[i] <- rbinom (1, 1.}}

9 488)}else if ( [i]=="identical twin"){girls[i] <- 2*rbinom (1, 1, .495)}else if ( [i]=="fraternal twin"){girls[i] <- rbinom (1, 2, .495)}} <- sum (girls)vectorized version of the loopgirls <- ifelse ( "single birth", rbinom (400, 1, .488),ifelse ( "identical twins", 2*rbinom (400, 1, .495),rbinom (400, 2, .495)))Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201518 / 48simulation for statistical inferenceusing replicate()repeat the simulation many <-function(x){ <- sample (c("fraternal twin","identical twin","single birth"),size=x, replace=TRUE, prob=c(1/125, 1/300, 1 - 1/125 - 1/300))girls <- ifelse ( "single birth", rbinom (400, 1, .488),ifelse ( "identical twins", 2*rbinom (400, 1, .495),rbinom (400, 2, .495)))return(sum(girls))} (400) <-replicate(1000, (400))hist( )Charles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201519 / 48simulation for statistical inferenceExercise 2.

10 Using simulation to Draw Statistical InferencesCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201520 / 48simulation for statistical inferenceconfidence ratio of two proportionsSurvey 1000 people, 500 men and 500 men support death penalty, 65% of men to women is = a standard error for this ratio can be the p qnCharles DiMaggio, PhD, MPH, PA-C (New york university department of Surgery and population health NYU-Bellevue Division of Trauma and Surgical Critical Care) Introduction to Simulations in RJune 10, 201521 / 48simulation for statistical inferenceconfidence intervalsdoing the <- <- <- sqrt ( *( ) ) <- <- <- sqrt ( *( ) )# Run 10,000 normal Simulations for each <- <- rnorm ( , , ) <- rnorm ( , , )# ratio of the simulationratio <- # 95% CI of the <- quantile (ratio, c(.))