+
+
+
+
+
+
+
+
8 Monte Carlo Simulation Exercise
+
Question a and b
+
set.seed(9252568)
+
+# n.sample is number of samples that we wanna sample
+# size is the sample size of each sample
+# mean is the sample mean of each sample that we sampled form
+# sd is the sample sd of each sample that we sampled from
+
+sample_samples <- function(n.sample, size, mean, sd){
+
+ # create the matrix to store each values in each sample
+ samples <- matrix(0, nrow = size, ncol = n.sample)
+ colnames(samples) <- c(1:n.sample)
+
+ # create the matrix to store means, absolute bias, SE, 95% CI, sample bias for each sample
+ info_sample <- data.frame(matrix(0, nrow = n.sample, ncol = 6))
+ colnames(info_sample) <- c("means", "abs_bias", "SE", "low_95_CI", "high_95_CI",
+ "sample_bias")
+
+ # sample samples
+ for (i in 1:n.sample) {
+ samples[,i] <- rnorm(n = size, mean = mean, sd = sd)
+ }
+
+ # display the info of each sample
+ for (i in 1:n.sample) {
+ info_sample[i,1] <- mean(samples[,i])
+ info_sample[i,2] <- abs(info_sample[i,1] - mean)
+ info_sample[i,3] <- sd(samples[,i]) / sqrt(size)
+ info_sample[i,4] <- mean(samples[,i]) - 1.96*info_sample[i,3]
+ info_sample[i,5] <- mean(samples[,i]) + 1.96*info_sample[i,3]
+ info_sample[i,6] <- ifelse((info_sample[i,4] > mean) | (info_sample[i,5] < mean),
+ "Biased", "Correct")
+ }
+
+ return(info_sample)
+}
+
+# draw the samples
+my_draw <- sample_samples(100, 1000, 2, 0.8)
+
Question c
+
library(tidyverse)
+
ggplot(data = my_draw, aes(x = 1:100, y = means)) +
+ geom_point(color = "#585858") +
+ geom_errorbar(aes(ymin = low_95_CI, ymax = high_95_CI, colour = sample_bias)) +
+ scale_color_manual(values=c("Blue", "#383838"), name = "Biased Sample", ) +
+ geom_hline(yintercept = 2, color = "Red", linetype = "dashed") +
+ labs(x="Simulations", y = "95% Confident Interval") +
+ theme_classic()
+
+
Question d
+
my_draw %>% filter(sample_bias == "Biased")
+
## means abs_bias SE low_95_CI high_95_CI sample_bias
+## 1 1.946501 0.05349950 0.02553148 1.896459 1.996542 Biased
+## 2 1.936868 0.06313216 0.02524686 1.887384 1.986352 Biased
+## 3 2.057306 0.05730614 0.02516567 2.007981 2.106631 Biased
+## 4 2.055018 0.05501774 0.02533916 2.005353 2.104682 Biased
+## 5 1.947215 0.05278480 0.02577351 1.896699 1.997731 Biased
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Introduction
+
The study is concerned with whether the number of rat sighting is
+tied to the population of a city and the waste of food thrown on the
+streets.
+
The data analyzed are simulated with 4 variables, they are:
+
city_id: ID number of each city popu: population of
+a city in millions. was_food: the amount of food thrown on the
+street an year, in tons rat_sight: the number of rat sighting
+an year
+
library(truncnorm) # for function rtruncnorm() that allow to set the max and min for normal dist
+options(scipen = 999)
+
Simulate the data
+
set.seed(9252568)
+city_id <- c(1:2000)
+popu <- rtruncnorm(n = 2000, a = 0.05, mean = 1, sd = 1)
+was_food <- rtruncnorm(2000, a = 10, mean = 30, sd = 20)
+
+rat_sight <- round(700.5 + 8.7*popu + 12.7*was_food + rnorm(2000, mean = 0.5),1)
+
+rat <- data.frame(city_id, popu, was_food, rat_sight)
+
Analysis
+
# operate the linear regression for analyzing
+rat_model <- lm(rat_sight~ popu + was_food, data = rat)
+
+summary(rat_model)
+
##
+## Call:
+## lm(formula = rat_sight ~ popu + was_food, data = rat)
+##
+## Residuals:
+## Min 1Q Median 3Q Max
+## -2.8992 -0.6873 0.0055 0.6723 3.5858
+##
+## Coefficients:
+## Estimate Std. Error t value Pr(>|t|)
+## (Intercept) 701.046782 0.067300 10416.7 <0.0000000000000002 ***
+## popu 8.675431 0.028409 305.4 <0.0000000000000002 ***
+## was_food 12.699526 0.001381 9194.8 <0.0000000000000002 ***
+## ---
+## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+##
+## Residual standard error: 0.9905 on 1997 degrees of freedom
+## Multiple R-squared: 1, Adjusted R-squared: 1
+## F-statistic: 4.229e+07 on 2 and 1997 DF, p-value: < 0.00000000000000022
+
+
+
Replication
+
Simulation
+
set.seed(92525568)
+
+city_id <- c(1:2000)
+popu <- rtruncnorm(n = 2000, a = 0.05, mean = 1, sd = 1)
+was_food <- rtruncnorm(2000, a = 10, mean = 30, sd = 20)
+
+rat_sight <- round(700.5 + 8.7*popu + 12.7*was_food + rnorm(2000, mean = 0.5),1)
+
+rat_rep <- data.frame(city_id, popu, was_food, rat_sight)
+
Analysis
+
rep <- lm(rat_sight~ popu + was_food, data = rat_rep)
+summary(rep)
+
##
+## Call:
+## lm(formula = rat_sight ~ popu + was_food, data = rat_rep)
+##
+## Residuals:
+## Min 1Q Median 3Q Max
+## -3.8009 -0.6931 0.0073 0.6500 3.3009
+##
+## Coefficients:
+## Estimate Std. Error t value Pr(>|t|)
+## (Intercept) 701.024240 0.066896 10479 <0.0000000000000002 ***
+## popu 8.688244 0.028120 309 <0.0000000000000002 ***
+## was_food 12.698665 0.001419 8949 <0.0000000000000002 ***
+## ---
+## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+##
+## Residual standard error: 0.9992 on 1997 degrees of freedom
+## Multiple R-squared: 1, Adjusted R-squared: 1
+## F-statistic: 4.005e+07 on 2 and 1997 DF, p-value: < 0.00000000000000022
+
The results are quite similar after changing the seed.
+
sessionInfo()
+
## R version 4.2.1 (2022-06-23 ucrt)
+## Platform: x86_64-w64-mingw32/x64 (64-bit)
+## Running under: Windows 10 x64 (build 22000)
+##
+## Matrix products: default
+##
+## locale:
+## [1] LC_COLLATE=English_United States.utf8
+## [2] LC_CTYPE=English_United States.utf8
+## [3] LC_MONETARY=English_United States.utf8
+## [4] LC_NUMERIC=C
+## [5] LC_TIME=English_United States.utf8
+##
+## attached base packages:
+## [1] stats graphics grDevices utils datasets methods base
+##
+## other attached packages:
+## [1] truncnorm_1.0-8
+##
+## loaded via a namespace (and not attached):
+## [1] digest_0.6.29 R6_2.5.1 jsonlite_1.8.0 magrittr_2.0.3
+## [5] evaluate_0.16 stringi_1.7.8 cachem_1.0.6 rlang_1.0.4
+## [9] cli_3.3.0 rstudioapi_0.14 jquerylib_0.1.4 bslib_0.4.0
+## [13] rmarkdown_2.16 tools_4.2.1 stringr_1.4.1 xfun_0.31
+## [17] yaml_2.3.5 fastmap_1.1.0 compiler_4.2.1 htmltools_0.5.3
+## [21] knitr_1.40 sass_0.4.2
+
+
+
+
+
+
+