R Notebook

library(tidyverse) library(moderndive) library(ggplot2) library(skimr) library(gapminder) library(infer)

set.seed(4869)

#bringing in dataframe, has 4969 entries with 7 columns movie_boxoffice.1 <- read.csv(“C:/Users/zahir/OneDrive/Desktop/Rutgers Semester Folders/2024 Fall Semester/Inference Data Science/HW 4 Lab/movie_boxoffice-1.csv”) View(movie_boxoffice.1)

#cleaning up dataframe of duplicates, has 4869 entries with 7 columns movie_boxoffice_distinct <- movie_boxoffice.1 %>% distinct() view(movie_boxoffice_distinct)

#0. How many duplicated records were removed from boxoffice data? What is your random seed number? #Answer: 100 dublicated records(entries) were removed, and my random seed number is set to 4869

#************************************************************************************************ #Part I: Sampling #1.Treat the box office data as population. Have a histogram of the global box office earning. Describe shape of distribution.

ggplot(movie_boxoffice_distinct, aes(x = Worldwide_Gross)) + geom_histogram(fill = "lightblue", color = "black") + 
labs(title = "Worldwide Box Office Earnings Distribution", x = "Worldwide Gross", y = "Count")

#Answer: Seems to be an exponential distribution with mode occuring at 0.

#2. In population, what is average global box earning? What is the std? What is the proportion of movies who global box earning exceeds budget?

average_globalbox_earning <- mean(movie_boxoffice_distinct$Worldwide_Gross)
print(average_globalbox_earning)

[1] 95.83149

std_globalbox_earning <- sd(movie_boxoffice_distinct$Worldwide_Gross)
print(std_globalbox_earning)

[1] 177.4594

global_earning_exceeds_budget <- movie_boxoffice_distinct %>%
  filter(Worldwide_Gross > Budget)
proportion_exceed <- nrow(global_earning_exceeds_budget) / nrow(movie_boxoffice_distinct)
print(proportion_exceed)

[1] 0.6461286

#Answer:Average global box earning is $95.83 million. Standard deviation is $177.46 million. Proportion of global box earning exceeding budget is 64.61%

#3. Take a random sample of 200 movies from the pop. Get histogram of global box earning, describe shape of distribution. set.seed(4869) movies_sample <- movie_boxoffice_distinct %>% sample_n(size = 200)

ggplot(movies_sample, aes(x = Worldwide_Gross)) + geom_histogram(fill = "lightblue", color = "black") + 
labs(title = "Worldwide Box Office Earnings Distribution", x = "Worldwide Gross", y = "Count")

#Answer: Seems to be an exponential distribution with mode occuring at 0. Also seems to be similar to the overall population histogram.

#4. In sample, what {is avg global box office earning? is std? is proportion of global box movie earning exceeding budget?} Are these summary stats from sample close to population paramets?

average_sample_earning <- mean(movies_sample$Worldwide_Gross)
print(average_sample_earning)

[1] 136.5527

std_sample_earning <- sd(movies_sample$Worldwide_Gross)
print(std_sample_earning)

[1] 231.5467

global_sample_earning_exceeds_budget <- movies_sample %>%
  filter(Worldwide_Gross > Budget)
proportion_sample_exceed <- nrow(global_sample_earning_exceeds_budget) / nrow(movie_boxoffice_distinct)
print(proportion_sample_exceed)

[1] 0.02813719

#Answer:Average global box earning is $136.55 million. Standard deviation is $231.55 million. Proportion of global box earning exceeding budget is 2.81%. These summary stats differ quite a bit from the population summary stats. Specifically about $40 million in avg, $54 million in std, and about 62% difference in proportion. This sample may contain more movies with lesser global box earnings than their budgets unlike the overall population explaining such low proportion.

#5.Take rand sample of n movies from the pop, calculate the avg global box office earning, and proportion of global box earning exceeding budget.

set.seed(4869)
sample_stats <- function(n, reps = 500) {
  movie_boxoffice_distinct %>%
    rep_sample_n(size = n, reps = reps) %>%
    group_by(replicate) %>%
    summarize(
      gross_avg = mean(Worldwide_Gross),
      proportion_ex = mean(Worldwide_Gross > Budget),
    )
}

n_val <- sample(50:500, 1)

output_ex <- sample_stats(n_val)
print(output_ex)

overall_avg <- mean(output_ex$gross_avg)
cat(
  "Sample Size =", n_val, 
  ": AVG Global Box Office Earnings =", round(overall_avg, 2),
  ", Proportion Exceeding Budget =", round(mean(output_ex$proportion_ex, na.rm = TRUE), 4)
  )

Sample Size = 78 : AVG Global Box Office Earnings = 95.06 , Proportion Exceeding Budget = 0.6478

#6 For each n=20, 50, 100, and 200, get histrogram of the avg global box earning using face_wrap(), also get mean and std error

set.seed(4869)
get_repetition_stats <- function(n, reps = 500) {
  movie_boxoffice_distinct %>%
    rep_sample_n(size = n, reps = reps) %>%
    group_by(replicate) %>%
      summarize(
        gross_mean = mean(Worldwide_Gross),
        proportion_exceed = mean(Worldwide_Gross > Budget)
      )
}

n_20 <- get_repetition_stats(20)
n_50 <- get_repetition_stats(50)
n_100 <- get_repetition_stats(100)
n_200 <- get_repetition_stats(200)

n_20$Sample_Size <- 20
n_50$Sample_Size <- 50
n_100$Sample_Size <- 100
n_200$Sample_Size <- 200
  
allSamplesDF <- bind_rows(n_20, n_50, n_100, n_200)

allSamplesDF$Sample_Size <- rep(c(20, 50, 100, 200), each = 500)

ggplot(allSamplesDF, aes(x = gross_mean)) + geom_histogram(fill = "lightblue", color = "black") + facet_wrap(~ Sample_Size, scales = "free") + labs(title = "Average Global Box Office Earnings Distribution by n Size", x = "Worldwide Gross Mean", y ="Count")

summary_stats <- allSamplesDF %>%
  group_by(Sample_Size) %>%
  summarize(
    mean_allgross = mean(gross_mean),
    std_err_allgross = sd(gross_mean) / sqrt(n())
  )
print(summary_stats)

#7. Compare distributions among the different n values, and also compare them to the distribution from the population in Q1 #Answer: For distrubtions among the different n values, starting at n = 20, more closely aligned with the distribution from the population in Q1, but mode no longer being at 0, but a value around $70 million. Going up the scale on the n values to 50, 100, and 200, the distribtuion was approaching normality, especially with n being 200. Although mode count decreased, the concisesness of it went up to around the value of $90 million.

#8. For each n = {20, 50, 100, 200}, get the histogram of the proportion of movies whose global box office earning exceeds budget, using face_wrap, also get mean and std error

set.seed(4869)
get_repetition_prop <- function(n, reps = 500) {
  movie_boxoffice_distinct %>%
    rep_sample_n(size = n, reps = reps) %>%
    group_by(replicate) %>%
      summarize(
        proportion_exceed = mean(Worldwide_Gross > Budget)
      )
}

pn_20 <- get_repetition_prop(20)
pn_50 <- get_repetition_prop(50)
pn_100 <- get_repetition_prop(100)
pn_200 <- get_repetition_prop(200)

pn_20$Sample_Size <- 20
pn_50$Sample_Size <- 50
pn_100$Sample_Size <- 100
pn_200$Sample_Size <- 200
  
allSamDF <- bind_rows(n_20, n_50, n_100, n_200)

allSamDF$Sample_Size <- rep(c(20, 50, 100, 200), each = 500)

ggplot(allSamDF, aes(x = proportion_exceed)) + geom_histogram(fill = "lightblue", color = "black") + facet_wrap(~ Sample_Size, scales = "free") + labs(title = "Proportion of Successful Movies by n Size", x = "Proportion of Worldwide Gross Exceeding Budget", y ="Count")

prop_summary_stats <- allSamDF %>%
  group_by(Sample_Size) %>%
  summarize(
    mean_proportion = mean(proportion_exceed),
    std_err_allgross = sd(proportion_exceed) / sqrt(n())
  )
print(prop_summary_stats)

#9. Compare the distributions among different n values. #Answer: For n = 20 describles a somewhat left-skewed distribtuion with mode being at 0.7. For n = 50, distribtuion appears to be normal, while at n = 100, repeating similarities with n = 20, but now there are two modes making it a bimodal left-skewed distribution. For N = 200, distribution appears to be more normal similar to n = 50.

#************************************************************************************************ #Part II: Bootstrapping

#10. Use data from Q3 as the initial sample, use bootstrapping method to resample once with 200 movies. Are there any duplicated movies, expected or not?

set.seed(4869)

bootstrap_samp <- movies_sample %>%
  rep_sample_n(size = 200, replace = TRUE) %>%
  summarize(gross_avg = mean(Worldwide_Gross))
print(bootstrap_samp)

set.seed(4869)
duplicateNum <- movies_sample %>%
  rep_sample_n(size = 200, replace = TRUE) %>%
  count(Movie) %>%
  filter(n > 1)
print(duplicateNum)

#Answer:Yes, there were duplicated movies, was this expected? Of course it was expected because we set replace to be TRUE allowing picked movies to be picked again.

#11.Get the histrogram of global box office earning in the bootstrap sample. Describe the shape of distribution and compare it to Q3

bootstrap2_samp <- movies_sample %>%
  rep_sample_n(size = 200, replace = TRUE)

ggplot(bootstrap2_samp, aes(x = Worldwide_Gross)) + geom_histogram(fill = "lightblue", color = "black") + labs(title = "Boostrap Sample Histogram of Worldwide Box Earnings", x = "Worldwide Gross", y ="Count")

#Answer: The shape of bootstrap sample histogram is closer to the exponential distribution, with mode still occuring at 0, and more variance with higher counts of outliers than as seen in Q3.

#12.In the bootstrap sample from Q10, what {is the avg global box earning, is the std, is the proportion of succesful movies?} Are they close enough to those in the intial sample?

avg_gross <- mean(bootstrap2_samp$Worldwide_Gross)
print(avg_gross)

[1] 127.9658

std_gross <- sd(bootstrap2_samp$Worldwide_Gross)
print(std_gross)

[1] 221.3411

prop_ex_bu <- mean(bootstrap2_samp$Worldwide_Gross > bootstrap2_samp$Budget)
print(prop_ex_bu)

[1] 0.69

#Answer: Average gross given is $127.97 million, Std dev is $221.34 million, and proportion of successful movies is 69%. These are pretty close to the intial values of avg gross being $136.55 million, std dev as $231.55 million, with the biggest defference being proportion from 2.81%. However, this is now closer to the overall population proportion of 64.61.

#13. Get the boostrapping distr of the avg global box earning, and proportion of movies whose global box earning exceeded budget, by resampling 500 times with bootstrapping method. Get the mean and std error of the average global box earning and mean and std error for proportion.

bootstrap_samples <- movies_sample %>%
  rep_sample_n(size = nrow(movies_sample), reps = 500, replace = TRUE)

bootstrap_stats <- bootstrap_samples %>%
  summarize(
    bootstrap_grossmean = mean(Worldwide_Gross),
    bootstrap_prop = mean(Worldwide_Gross > Budget)
  )
bootstrap_grossavg <- mean(bootstrap_stats$bootstrap_grossmean)
print(bootstrap_grossavg)

[1] 120.3005

bootstrap_std_grossavg <- sd(bootstrap_stats$bootstrap_grossmean)
print(bootstrap_std_grossavg)

[1] 14.59414

mean_prop <- mean(bootstrap_stats$bootstrap_prop)
print(mean_prop)

[1] 0.65442

se_prop <- sd(bootstrap_stats$bootstrap_prop) / sqrt(nrow(bootstrap_stats))
print(se_prop)

[1] 0.001605641

ggplot(bootstrap_stats, aes(x = bootstrap_grossmean)) + geom_histogram(fill = "lightblue", color = "black") + labs(title = "Average Global Box Office Earnings Distribution", x = "AVG Worldwide Gross", y = "Count")

ggplot(bootstrap_stats, aes(x = bootstrap_prop)) + geom_histogram(fill = "lightblue", color = "black") + labs(title = "Proportion of Successful Movies", x = "Successful Movies", y = "Count")

#Answer: Given mean of Average global office earning is $120.30 million, with std dev being $14.59 million. Proportion of succesful movies given is 65.44%, and the standard error for it is 0.16%. The distribution for Average Global Box Earnings is right leaning with occuring around $112 million, and for the poportion of succesful moves, distribution is more normal with a slight left skew and mode occuring around 67.5%.