HW7: Creating Fake Data Sets To Explore Hypotheses

library(tidyverse) # To use ggplot2 and dplyr packages
library(kableExtra) # To make fancy tables
library(gridExtra)

##############################################################################################
# FUNCTION: sim_rnorm()
# Description: a function that creates a simulated gene expression dataset for different culture durations.
# Input Arguments: 'list' (a list containing the mean, variance and size of each 'Time' categories)
# Output: 'sim_df' (a dataframe containing a unique ID variable, Time and sim_AGER (simulated from normal distribution))
#___________________________________________________________________________________________

sim_rnorm <- function(list = input){

#Times = c("24h", "+4h", "+24h", "5d")

sim24h = rnorm(n = list[[1]][1], mean = list[[1]][2], sd = sqrt(list[[1]][3]))
sim48h = rnorm(n = list[[2]][1], mean = list[[2]][2], sd = sqrt(list[[2]][3]))
sim72h = rnorm(n = list[[3]][1], mean = list[[3]][2], sd = sqrt(list[[3]][3]))
sim144h = rnorm(n = list[[4]][1], mean = list[[4]][2], sd = sqrt(list[[4]][3]))

sim_AGER <- c(sim24h, sim48h, sim72h, sim144h)
ID <- 1:length(sim_AGER)
Time <- c (rep(24, length(sim24h)), rep(48, length(sim48h)), rep(72, length(sim72h)), rep(144, length(sim144h))) # Assigning Time categories to the rnorm() generated data

sim_df <- data.frame(ID, Time, sim_AGER) # Create a dataframe from the simulated Time and AGER data. 
  
return(sim_df)
}
##############################################################################################


##############################################################################################
# FUNCTION: Reg_Stat()
# Description: Performs linear regression on a dataset
# Input Arguments: 'df' (a dataframe containing the variable the linear regression will be ran on)
# Output: ‘a’ (an atomic vector containing the coefficient 'b' (Response = a + bPredictor), the p-value and conclusion (significant or insignificant))
#___________________________________________________________________________________________

Reg_Stat <- function(df = DF){
  
   z <- lm(sim_AGER ~ Time, data = df)
   coeff<-summary(z)$coefficients[2,1] # pulls out the coeffiecient of the relationship btw the response and predictor variables
   pVal <- summary(z)$coefficients[2,4] # pulls out the p-value of the relationship btw the response and predictor variables
  
   a <-c(signif(coeff, 4), signif(pVal, 4), ifelse(pVal < 0.05, "Significant", "Non-significant"))
   names(a) <- c("Coefficient", "P-Value", "Conclusion")
   return(a)
   }
##############################################################################################


##############################################################################################
# FUNCTION: Reg_Plot()
# Description: Plot the datset
# Input Arguments: ‘df’ (a dataframe containing the variable to be plotted)
# Output:'plot'
#___________________________________________________________________________________________

Reg_Plot <- function(df = DF){
  
   ### Plot of the data
   plot <- ggplot(df) +
           aes(x = Time, y = sim_AGER) +
           geom_point() +
           stat_smooth(method=lm,se=0.99) # default se=0.95 
   
  return(plot)
 }
##############################################################################################

df <- data.frame(read.csv("RawData.csv")) # read in the original csv dataset into a dataframe

a = lm(df$AGER ~ df$Time, data = df)
coeff<-summary(a)$coefficients[2,1]
pVal <- summary(a)$coefficients[2,4]
cat("The coefficient of the linear regression is: ", coeff, "\n",
    "The p-value is: ", pVal, "\n")

Subset <- df %>% group_by(Time) %>% 
                 summarize(Mean = mean(AGER, na.rm = T), Variance = sd(AGER, na.rm = T)^2, Size = n()) # Get the size, mean and variance of each Time group

kable(Subset, caption = "*Table 1 - Mean, Variance and Size for Each Time Group*", row.names = NA, align = "c")

## Create input variables needed to create simulated dataset

n_24h = 84; mean_24h = 0.8439886; var_24h = 0.0604894;
n_48h = 86; mean_48h = 2.2456403; var_48h = 3.0218083;
n_72h = 84; mean_72h = 2.7526489; var_72h = 7.8850666;
n_144h = 82; mean_144h = 2.4733090; var_144h = 9.2351024;

# Create a list with all input data

input <- list(sim24h = c(n_24h, mean_24h,var_24h), 
              sim48h = c(n_48h, mean_48h,var_48h),
              sim72h = c(n_72h, mean_72h,var_72h),
              sim144h = c(n_144h, mean_144h,var_144h))


df1 <- sim_rnorm(input) # create a randomly generated Time and AGER dataset
LR1 <- Reg_Stat(df1) # computes the linear regression analysis of AGER(response) and Time (predictor)

df2 <- sim_rnorm(input) # create a randomly generated Time and AGER dataset
LR2 <- Reg_Stat(df2) # computes the linear regression analysis of AGER(response) and Time (predictor)

df3 <- sim_rnorm(input) # create a randomly generated Time and AGER dataset
LR3 <- Reg_Stat(df3) # computes the linear regression analysis of AGER(response) and Time (predictor)

df4 <- sim_rnorm(input) # create a randomly generated Time and AGER dataset
LR4 <- Reg_Stat(df4) # computes the linear regression analysis of AGER(response) and Time (predictor)

Test <- c("1", "2", "3", "4")
RegStats1 <- bind_rows(LR1, LR2, LR3, LR4)
s <- cbind(Test,RegStats1)

kable(s, caption = "**Table 2 - Coefficient, p-Value and Conclusion from Linear Regression Analysis on 4 Different Randomly Generated Dataset**", align = "c")

plotList <-list(Reg_Plot(df1), Reg_Plot(df2), Reg_Plot(df3), Reg_Plot(df4)) # compress all the plots into a list

caption <- "Fig 1 - Plot of AGER (simulated from normal distribution with the same means, variance and sample size) and Time."

suppressMessages(grid.arrange(grobs = plotList, bottom = caption)) # plot all the plots on a grid

Reg_output <- NULL
trial_mean <- NULL

h <- seq(from = 1.5, to = 2.5, by = 0.05) # create a sequence of means for one group

for(i in seq_along(h)){
   
    n_24h = 84; mean_24h = h[i]; var_24h = 0.0604894;
    n_48h = 86; mean_48h = 2.5; var_48h = 3.0218083;
    n_72h = 84; mean_72h = 2.5; var_72h = 7.8850666;
    n_144h = 82; mean_144h = 2.5; var_144h = 9.2351024;
    
    # Create a list with all input data
    input2 <- list(sim24h = c(n_24h, mean_24h,var_24h), 
                  sim48h = c(n_48h, mean_48h,var_48h),
                  sim72h = c(n_72h, mean_72h,var_72h),
                  sim144h = c(n_144h, mean_144h,var_144h))
    
    mean_diff <- mean_48h - mean_24h # compute mean difference
    set.seed(99)
    dfz <- sim_rnorm(input2) # create a randomly generated Time and AGER dataset
    
    RegOutput <- data.frame(t(Reg_Stat(dfz))) # computes the linear regression analysis of AGER(response) and Time
    tf <- data.frame(t(c(i, mean_diff))) # create a dataframe with the trial no. and the mean diff. for the current trial
    
    trial_mean <- rbind(trial_mean, tf)    # update the 'trial_mean' dataframe
    Reg_output <-rbind(Reg_output, RegOutput) # update the 'Reg_output' dataframe
    
}

names(trial_mean) <- c("Trial No", "Mean Difference")
RegResults <- cbind(trial_mean,Reg_output)
output <- RegResults %>% select(-Coefficient)

kable(output, caption = "**Table 3 - LR Output with Varying Group Means Difference **", align = "c")

Reg_output <- NULL
trial_size <- NULL

h <- seq(from = 0, to = 84, by = 2) # create a sequence of means for one group

for(i in seq_along(h)){
   
    n_24h = h[i]; mean_24h = 0.8439886; var_24h = 0.0604894;
    n_48h = 84; mean_48h = 2.2456403; var_48h = 3.0218083;
    n_72h = 84; mean_72h = 2.7526489; var_72h = 7.8850666;
    n_144h = 84; mean_144h = 2.4733090; var_144h = 9.2351024;
        
    # Create a list with all input data
    input2 <- list(sim24h = c(n_24h, mean_24h,var_24h), 
                  sim48h = c(n_48h, mean_48h,var_48h),
                  sim72h = c(n_72h, mean_72h,var_72h),
                  sim144h = c(n_144h, mean_144h,var_144h))
    
    size_diff <- n_48h - n_24h # compute mean difference
    set.seed(99)
    dfz <- sim_rnorm(input2) # create a randomly generated Time and AGER dataset
    
    RegOutput <- data.frame(t(Reg_Stat(dfz))) # computes the linear regression analysis of AGER(response) and Time
    tf <- data.frame(t(c(i, size_diff))) # create a dataframe with the trial no. and the mean diff. for the current trial
    
    trial_size <- rbind(trial_size, tf)    # update the 'trial_mean' dataframe
    Reg_output <-rbind(Reg_output, RegOutput) # update the 'Reg_output' dataframe
    
}

names(trial_size) <- c("Trial No", "Size Difference")
RegResults <- cbind(trial_size,Reg_output)
output <- RegResults %>% select(-Coefficient)

kable(output, caption = "**Table 4 - LR Output with Varying Group Size Difference (with seed set to 99)**", align = "c")

Reg_output <- NULL
trial_size <- NULL

h <- seq(from = 0, to = 84, by = 2) # create a sequence of means for one group

for(i in seq_along(h)){
   
    n_24h = h[i]; mean_24h = 0.8439886; var_24h = 0.0604894;
    n_48h = 84; mean_48h = 2.2456403; var_48h = 3.0218083;
    n_72h = 84; mean_72h = 2.7526489; var_72h = 7.8850666;
    n_144h = 84; mean_144h = 2.4733090; var_144h = 9.2351024;
        
    # Create a list with all input data
    input2 <- list(sim24h = c(n_24h, mean_24h,var_24h), 
                  sim48h = c(n_48h, mean_48h,var_48h),
                  sim72h = c(n_72h, mean_72h,var_72h),
                  sim144h = c(n_144h, mean_144h,var_144h))
    
    size_diff <- n_48h - n_24h # compute mean difference
    #set.seed(99)
    dfz <- sim_rnorm(input2) # create a randomly generated Time and AGER dataset
    
    RegOutput <- data.frame(t(Reg_Stat(dfz))) # computes the linear regression analysis of AGER(response) and Time
    tf <- data.frame(t(c(i, size_diff))) # create a dataframe with the trial no. and the mean diff. for the current trial
    
    trial_size <- rbind(trial_size, tf)    # update the 'trial_mean' dataframe
    Reg_output <-rbind(Reg_output, RegOutput) # update the 'Reg_output' dataframe
    
}

names(trial_size) <- c("Trial No", "Size Difference")
RegResults <- cbind(trial_size,Reg_output)
output <- RegResults %>% select(-Coefficient)

kable(output, caption = "**Table 5 - LR Output with Varying Group Size Difference (without seed set) **", align = "c")