diff --git a/Figures/Figure_1.png b/Figures/Figure_1.png new file mode 100755 index 0000000..3066bc9 Binary files /dev/null and b/Figures/Figure_1.png differ diff --git a/Figures/Figure_2.png b/Figures/Figure_2.png new file mode 100755 index 0000000..9df9d0a Binary files /dev/null and b/Figures/Figure_2.png differ diff --git a/Figures/Figure_3.png b/Figures/Figure_3.png new file mode 100644 index 0000000..1fcd878 Binary files /dev/null and b/Figures/Figure_3.png differ diff --git a/Figures/Figure_4.png b/Figures/Figure_4.png new file mode 100644 index 0000000..1a264a4 Binary files /dev/null and b/Figures/Figure_4.png differ diff --git a/Figures/Figure_5.pdf b/Figures/Figure_5.pdf new file mode 100644 index 0000000..670693e Binary files /dev/null and b/Figures/Figure_5.pdf differ diff --git a/Final Writeup/BIOS60318-FinalProject-Writeup.pdf b/Final Writeup/BIOS60318-FinalProject-Writeup.pdf new file mode 100644 index 0000000..cd646ae Binary files /dev/null and b/Final Writeup/BIOS60318-FinalProject-Writeup.pdf differ diff --git a/PeerReviewComments.txt b/PeerReviewComments.txt new file mode 100644 index 0000000..e244d8f --- /dev/null +++ b/PeerReviewComments.txt @@ -0,0 +1,5 @@ +Comments for other group: +-Add more to the ggplot statements to make the graph more appealing +-When running ANOVA, "x" is not defined as antibiotic treatment +-Comment out that code so the whole script runs +-Consider building matrices using functions instead of manually inputting zeros and 1's diff --git a/RCode/10000SimulationEnvironment.RData b/RCode/10000SimulationEnvironment.RData new file mode 100755 index 0000000..1ccbdc8 Binary files /dev/null and b/RCode/10000SimulationEnvironment.RData differ diff --git a/RCode/BIOS60318_FinalProject.R b/RCode/BIOS60318_FinalProject.R new file mode 100755 index 0000000..aae7939 --- /dev/null +++ b/RCode/BIOS60318_FinalProject.R @@ -0,0 +1,383 @@ +########################################################## +########################################################## +################ BIOS 60318 Final Project ################ +########################################################## +########################################################## +rm(list=ls()) # Clear workspace +if(!is.null(dev.list())) dev.off() # Clear plots +par(mfrow=c(1,1)) # Setup plot parameters +par(ps = 12, font.lab = 1) # Plot parameters +set.seed(1) # set random seed generator for reproducibility +src = dirname(rstudioapi::getActiveDocumentContext()$path) # Get the path to where this script is located +setwd(src) # Set working directory to that path + +# If you don't have the packages necessary, let's make sure you do +list.of.packages <- c("ggplot2", "CorReg") +new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])] +if(length(new.packages)) install.packages(new.packages) + +# Let's load relevant libraries +# NOTE: Please make sure extractPVal.R and superNll.R are within your working directory +library(ggplot2) # Allow plotting capabilities +library(CorReg) # Allows use of modified BoxPlot() +source(paste0(getwd(),"/extractPVal.R"), encoding = 'UTF-8') # Custom function to extract a Pvalue from an lm() model +source(paste0(getwd(),"/superNll.R"), encoding = 'UTF-8') # Custom function to run MLLE on a data set in either regression or ANOVA format +cat("\014") # Clear console + +# Let's load the data +sugarData = read.table(paste0(getwd(),"/sugar.csv"), header=TRUE, sep=",") +antibioticsData = read.table(paste0(getwd(),"/antibiotics.csv"), header=TRUE, sep=",") + +########################################################## +########################################################## +##################### Part I: ANOVA ###################### +########################################################## +########################################################## +# Relevel the data to have the control as the reference +antibioticsData$trt = as.numeric(relevel(antibioticsData$trt, ref=4)) + +# Use a custom nll function to obtain parameters and pValue for log likelihood ratio test +antibioticNll = superNll(antibioticsData$trt, antibioticsData$growth, 4, length(antibioticsData$growth), anova = TRUE) + +# Fit using aov() that passes each group through lm() as a comparison +aov.fit = aov(growth ~ as.factor(trt), data=antibioticsData) + +# PValue from aov() using F statistic +aovPVal=summary(aov.fit)[[1]][["Pr(>F)"]][1] + +# Parameters found from nll +antibioticParam_nll = c(antibioticNll$coefficients, antibioticNll$sigma, antibioticNll$pValue) + +# Parameters found from lm() and aov() +antibioticParam_aov = c(aov.fit$coefficients, sigma(aov.fit), aovPVal) + +# Comparison of nll versus lm() +comparison_antibiotics = cbind(antibioticParam_nll, antibioticParam_aov) +dimnames(comparison_antibiotics)[[1]]=c("Control", "Treatment 1", "Treatment 2", "Treatment 3", "Residual Std Error", "P Value") +comparison_antibiotics +#write.table(comparison_antibiotics, "clipboard", sep="\t", row.names=FALSE) # Transport comparison to excel + +# 95% Confidence interval of the parameters using LSE +antibiotic_CI = confint(aov.fit, level = 0.95) +dimnames(antibiotic_CI)[[1]]=c("Control", "Treatment 1", "Treatment 2", "Treatment 3") +antibiotic_CI + +# 95% Confidence interval of the parameters using Hessian matrix and MLLE +invHessian_antibiotics = solve(antibioticNll$hessian) +parameter.std.error_antibiotics = sqrt(diag(invHessian_antibiotics)) +lowerBound_antibiotics = matrix(0, 4, 1) +upperBound_antibiotics = matrix(0, 4, 1) +plusMinus_antibiotics = matrix(0, 4, 1) +crit.value_antibiotics = qt(p = 0.025, df = length(antibioticsData$growth)) +for (i in 1:4){ + plusMinus_antibiotics[i] = crit.value_antibiotics*parameter.std.error_antibiotics[i] + lowerBound_antibiotics[i] = antibioticNll$coefficients[i]-plusMinus_antibiotics[i] + upperBound_antibiotics[i] = antibioticNll$coefficients[i]+plusMinus_antibiotics[i] +} +antibiotic_CI_MLLE = cbind(lowerBound_antibiotics, upperBound_antibiotics) +dimnames(antibiotic_CI_MLLE)[[1]] = c("Control", "Treatment 1", "Treatment 2", "Treatment 3") +antibiotic_CI_MLLE + + +# Boxplot of the data comparing control versus treatment with 95% CI of the means for each group in red +BoxPlot(antibioticsData$growth, + as.factor(antibioticsData$trt), + AnoVa = TRUE, ylab="Growth of Bacteria", + names=c("Control", "Treatment 1", "Treatment 2", "Treatment 3"), + verbose=FALSE) + +# Comparison of each treatment to the control for significance +# NOTE: Familywise error rate is not corrected for given the low number of tests performed +controlVSTrt1 = t.test(antibioticsData$growth[1:4], antibioticsData$growth[5:8]) +controlVSTrt2 = t.test(antibioticsData$growth[1:4], antibioticsData$growth[9:12]) +controlVSTrt3 = t.test(antibioticsData$growth[1:4], antibioticsData$growth[13:16]) + +########################################################## +########################################################## +############### Part II: Linear Regression ############### +########################################################## +########################################################## +colnames(sugarData)[2] = "growthSugar" +sugarNll = superNll(sugarData$sugar,sugarData$growthSugar,2,length(sugarData$sugar),anova=FALSE) + +# Combine all parameters for nll +sugarParam_nll = c(sugarNll$coefficients, sugarNll$sigma, sugarNll$pValue) + +# Linear regression via lm() +linear.mod = lm(growthSugar ~ sugar, data=sugarData) + +# All parameters from lm() +sugarParam_lm = c(linear.mod$coefficients, sigma(linear.mod), extractPVal(linear.mod)) + +# Comparison of the parameters +comparison_sugar = cbind(sugarParam_nll, sugarParam_lm) +dimnames(comparison_sugar)[[1]]=c("Beta 0", "Beta 1", "Residual Std Error", "P Value") +comparison_sugar + +# 95% confidence interval of the parameters b0 and b1 +CI_LM = confint(linear.mod) +dimnames(CI_LM)[[1]]=c("Beta 0", "Beta 1") +CI_LM + +# 95% Confidence interval of the parameters using Hessian matrix and MLLE +invHessian_sugar = solve(sugarNll$hessian) +parameter.std.error_sugar = sqrt(diag(invHessian_sugar)) +lowerBound_sugar = matrix(0, 2, 1) +upperBound_sugar = matrix(0, 2, 1) +plusMinus_sugar = matrix(0, 2, 1) +crit.value_sugar = qt(p = 0.025, df = length(sugarData$sugar)) +for (i in 1:2){ + plusMinus_sugar[i] = crit.value_sugar*parameter.std.error_sugar[i] + lowerBound_sugar[i] = sugarNll$coefficients[i]-plusMinus_sugar[i] + upperBound_sugar[i] = sugarNll$coefficients[i]+plusMinus_sugar[i] +} +sugar_CI_MLLE = cbind(lowerBound_sugar, upperBound_sugar) +dimnames(sugar_CI_MLLE)[[1]] = c("Intercept", "Slope") +sugar_CI_MLLE + +# Plot the data and fit +ggplot(sugarData,aes(sugarData$sugar,sugarData$growthSugar))+ + geom_point()+ + geom_abline(slope=linear.mod$coefficients[2],intercept=linear.mod$coefficients[1])+ + geom_abline(slope=0,intercept=3.33,color="red")+ + theme_classic()+ + xlab("Sugar Concentration")+ + ylab("Growth of E.Coli")+ + geom_smooth(method="lm") + +########################################################## +########################################################## +############### Part III: Power Analysis ################# +########################################################## +########################################################## + +# NOTE TO CHECK FUNCTION CODE CHANGE Nsim to 100 or below +Nsim=10000 # Number of simulations to run +beta0=10 # The intercept of the line +beta1=0.4 # The slope of the line +sigmaVals = c(1,2,4,6,8,12,16,24) # Standard Deviations for the error term + +N=24 # Number of experimental units +maxX = 50 # Maximum value of experimental data +minX = 0 # Minimum value of experimental data +nLevelsVals = c(2,4,8) # Number of levels in ANOVA + +# Define NULL vectors to hold data later in the code +paramComparison = NULL +paramComparison_nll = NULL +pMatrixAll = NULL +pMatrixRun = matrix(0, Nsim, 6) + +# Loop over all levels and all sigmas +for (nLevels in nLevelsVals){ + for (sigma in sigmaVals){ + + # Linear Regression + X=seq(minX,maxX,(maxX+1)/N) # Sequenced X values + + # Declare Linear Regression Matrices for lm() + Y=matrix(rep(N*Nsim),nrow=N,ncol=Nsim) # Create our Y matrix + coeff_matrix_lm=matrix(rep(2*Nsim),nrow=2,ncol=Nsim) # Create our coefficent matrix + sigma_matrix_lm=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our sigma matrix + p_matrix_lm=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our p value matrix + + # Declare Anova matrices for aov() + coeff_matrix_aov=matrix(rep(nLevels*Nsim),nrow=nLevels,ncol=Nsim) # Create our coefficent matrix for aov + sigma_matrix_aov=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our sigma matrix for aov + p_matrix_aov=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our p value matrix for aov + + # running the simulations using lm() and aov() + for (i in 1:Nsim){ + epsilon=rnorm(N,mean=0,sd=sigma) # Normalized error + Y[,i]=beta0+beta1*X+epsilon # Generate data points + + mod=lm(Y[,i]~X) # Fit the data using lm() + coeff_matrix_lm[,i]=mod$coefficients # Extract coefficients + sigma_matrix_lm[i]=summary(mod)$sigma # Extract sigmas + p_matrix_lm[i]=extractPVal(mod) # Extract P Value + + anovaResults=data.frame(levels=rep(seq(0,nLevels-1,1),each=N/nLevels),YVal=Y[,i]) + aovMod = aov(anovaResults$YVal~as.factor(anovaResults$levels)) + coeff_matrix_aov[,i]=aovMod$coefficients + sigma_matrix_aov[i]=sigma(aovMod) + p_matrix_aov[i]=summary(aovMod)[[1]][["Pr(>F)"]][1] + } + + # Calculate the averages of the simulations + averageb0_reg = mean(coeff_matrix_lm[1]) + averageb1_reg = mean(coeff_matrix_lm[2]) + averagep_reg = mean(p_matrix_lm) + significantP_reg = length(which(p_matrix_lm<0.05)) + + averageb0_aov = mean(coeff_matrix_aov[1]) + averageb1_aov = mean(coeff_matrix_aov[2]) + averagep_aov = mean(p_matrix_aov) + significantP_aov = length(which(p_matrix_aov<0.05)) + + # Compare the two parameters and P Values between lm() and aov() + linearParameters = c(averageb0_reg, averageb1_reg,averagep_reg,significantP_reg) + anovaParameters = c(averageb0_aov, averageb1_aov,averagep_aov,significantP_aov) + comparison_III = cbind(linearParameters,anovaParameters) + dimnames(comparison_III)[[1]]=c("Beta 0", "Beta 1", "P Value","Num Significant P") + paramComparison = cbind(paramComparison, sigma, comparison_III) + + ########################################################################################################## + + # Declare Linear Regression Matrices for custom nll function + Y_nll=matrix(rep(N*Nsim),nrow=N,ncol=Nsim) # Create our Y matrix + coeff_matrix_lm_nll=matrix(rep(2*Nsim),nrow=2,ncol=Nsim) # Create our coefficent matrix + sigma_matrix_lm_nll=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our sigma matrix + p_matrix_lm_nll=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our p value matrix + + # Declare Anova matrices for custom nll function + coeff_matrix_aov_nll=matrix(rep(nLevels*Nsim),nrow=nLevels,ncol=Nsim) # Create our coefficent matrix for aov + sigma_matrix_aov_nll=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our sigma matrix for aov + p_matrix_aov_nll=matrix(rep(Nsim),nrow=1,ncol=Nsim) # Create our p value matrix for aov + + # running the simulations using custom nll() functions + for (i in 1:Nsim){ + epsilon=rnorm(N,mean=0,sd=sigma) # Normalized error + Y_nll[,i]=beta0+beta1*X+epsilon # Generate data points + + mod_nll=superNll(X,Y[,i],nLevels = 2,nExpUnits = N,anova = FALSE) # Use nll to obtain parameter estimates for linear regression + coeff_matrix_lm_nll[,i]=mod_nll$coefficients # Extract coefficients + sigma_matrix_lm_nll[i]=mod_nll$sigma # Extract sigmas + p_matrix_lm_nll[i]=mod_nll$pValue # Extract P Value + + anovaResults=data.frame(levels=rep(seq(0,nLevels-1,1),each=N/nLevels),YVal=Y[,i]) # set up anova data frame and bin data + aovMod_nll = superNll(anovaResults$levels, anovaResults$YVal, nLevels, N,anova = TRUE) # use nll to obtain parameter estimates from anova regression + coeff_matrix_aov_nll[,i]=aovMod_nll$coefficients # Extract coefficients + sigma_matrix_aov_nll[i]=aovMod_nll$sigma # Extract sigmas + p_matrix_aov_nll[i]=aovMod_nll$pValue # Extract p values + } + + # Calculate the averages of the simulations + averageb0_reg_nll = mean(coeff_matrix_lm_nll[1]) + averageb1_reg_nll = mean(coeff_matrix_lm_nll[2]) + averagep_reg_nll = mean(p_matrix_lm_nll) + significantP_reg_nll = length(which(p_matrix_lm_nll<0.05)) + + averageb0_aov_nll = mean(coeff_matrix_aov_nll[1]) + averageb1_aov_nll = mean(coeff_matrix_aov_nll[2]) + averagep_aov_nll = mean(p_matrix_aov_nll) + significantP_aov_nll = length(which(p_matrix_aov_nll<0.05)) + + # Compare the two parameters and P Values + linearParameters_nll = c(averageb0_reg_nll, averageb1_reg_nll,averagep_reg_nll,significantP_reg_nll) + anovaParameters_nll = c(averageb0_aov_nll, averageb1_aov_nll,averagep_aov_nll,significantP_aov_nll) + comparison_III_nll = cbind(linearParameters_nll,anovaParameters_nll) + dimnames(comparison_III_nll)[[1]]=c("Beta 0", "Beta 1", "P Value","Num Significant P") + paramComparison_nll = cbind(paramComparison_nll, sigma, comparison_III_nll) + + # Combine all P Values for each MLLE and LSE run + pMatrixRun[,1] = p_matrix_lm + pMatrixRun[,2] = p_matrix_aov + pMatrixRun[,3] = p_matrix_lm_nll + pMatrixRun[,4] = p_matrix_aov_nll + pMatrixRun[,5] = sigma + pMatrixRun[,6] = nLevels + pMatrixAll = cbind(pMatrixAll, pMatrixRun) + } +} + +# Extract the two level comparison data for lm() vs. aov() +twoLevelANOVA = cbind(paramComparison[,1:24]) +# Extract the four level comparison data for lm() vs. aov() +fourLevelANOVA = cbind(paramComparison[,25:48]) +# Extract the eight level comparison data for lm() vs. aov() +eightLevelANOVA = cbind(paramComparison[,49:72]) + +# Extract the two level comparison data for linear regression nll vs. anova nll +twoLevelANOVA_nll = cbind(paramComparison_nll[,1:24]) +# Extract the four level comparison data for linear regression nll vs. anova nll +fourLevelANOVA_nll = cbind(paramComparison_nll[,25:48]) +# Extract the eight level comparison data for linear regression nll vs. anova nll +eightLevelANOVA_nll = cbind(paramComparison_nll[,49:72]) + +# Plot the histogram of the p-values for lm() vs. aov() for 8 levels and sigma = 24 (MSE) +# LSE For linear regression p distribution +ggplot(data.frame(pVal=t(p_matrix_lm)), aes(x=pVal))+ + geom_histogram(binwidth = 0.1, color="white",fill="blue")+ + theme_classic()+ + xlab("p-Values")+ + ylab("Frequency")+ + ggtitle(label="p-Value distribution of lm()")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2) + +# LSE For ANOVA p distribution +ggplot(data.frame(pVal=t(p_matrix_aov)), aes(x=pVal))+ + geom_histogram(binwidth = 0.1, color="white",fill="blue")+ + theme_classic()+ + xlab("p-Values")+ + ylab("Frequency")+ + ggtitle(label="p-Value distribution of aov()")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2) + +# Plot the histogram of the p-values for linear nll vs. anova nll for 8 levels and sigma = 24 (MLLE) +# MLLE For linear regression p distribution +ggplot(data.frame(pVal=t(p_matrix_lm_nll)), aes(x=pVal))+ + geom_histogram(binwidth = 0.1, color="white",fill="blue")+ + theme_classic()+ + xlab("p-Values")+ + ylab("Frequency")+ + ggtitle(label="p-Value distribution of linear regresssion nll()")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2) + +# MLLE For ANOVA p distribution +ggplot(data.frame(pVal=t(p_matrix_aov_nll)), aes(x=pVal))+ + geom_histogram(binwidth = 0.1, color="white",fill="blue")+ + theme_classic()+ + xlab("p-Values")+ + ylab("Frequency")+ + ggtitle(label="p-Value distribution of anova nll()")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2) + +# Plot the histogram of the p-values for lm() vs. aov() for 2 levels and sigma = 1 (LSE) +# LSE For linear regression p distribution +ggplot(data.frame(pVal=-log10(pMatrixAll[,1])), aes(x=pVal))+ + geom_histogram(binwidth = 0.1, color="white",fill="blue")+ + theme_classic()+ + xlab("-log10(p-Values)")+ + ylab("Frequency")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2)+ + xlim(min(-log10(pMatrixAll[,1])),max(-log10(pMatrixAll[,1]))) + +# LSE For ANOVA p distribution +ggplot(data.frame(pVal=-log10(pMatrixAll[,2])), aes(x=pVal))+ + geom_histogram(binwidth = 0.05, color="white",fill="blue")+ + theme_classic()+ + xlab("-log10(p-Values)")+ + ylab("Frequency")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2)+ + xlim(min(-log10(pMatrixAll[,2])),max(-log10(pMatrixAll[,2]))) + +# Plot the histogram of the p-values for custom nll for Linear regression and ANOVA for 2 levels and sigma = 1 (MLLE) +# MLLE For linear regression p distribution +# NOTE: 1e24 was added to all values due to the values being below the machine computable epsilon value +ggplot(data.frame(pVal=-log10(pMatrixAll[,3]+1e-24)), aes(x=pVal))+ + geom_histogram(binwidth = 0.1, color="white",fill="blue")+ + theme_classic()+ + xlab("-log10(p-Values)")+ + ylab("Frequency")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2)+ + xlim(23,25) + +# MLLE For ANOVA p distribution +ggplot(data.frame(pVal=-log10(pMatrixAll[,4])), aes(x=pVal))+ + geom_histogram(binwidth = 0.05, color="white",fill="blue")+ + theme_classic()+ + xlab("-log10(p-Values)")+ + ylab("Frequency")+ + theme(plot.title=element_text(hjust=0.5))+ + geom_vline(aes(xintercept=mean(pVal)),color="red",size=1.2)+ + xlim(min(-log10(pMatrixAll[,4])),max(-log10(pMatrixAll[,4]))) + +#write.table(eightLevelANOVA, "clipboard", sep="\t", row.names=TRUE) # Copy data to clipboard +#write.table(eightLevelANOVA_nll, "clipboard", sep="\t", row.names=TRUE) # copy data to clipboard diff --git a/antibiotics.csv b/RCode/antibiotics.csv old mode 100644 new mode 100755 similarity index 100% rename from antibiotics.csv rename to RCode/antibiotics.csv diff --git a/RCode/extractPVal.R b/RCode/extractPVal.R new file mode 100755 index 0000000..c703275 --- /dev/null +++ b/RCode/extractPVal.R @@ -0,0 +1,6 @@ +extractPVal <- function (linearFit) { + fStat = summary(linearFit)$fstatistic + pVal = pf(fStat[1],fStat[2],fStat[3],lower.tail=F) + attributes(pVal) = NULL + return(pVal) +} diff --git a/sugar.csv b/RCode/sugar.csv old mode 100644 new mode 100755 similarity index 100% rename from sugar.csv rename to RCode/sugar.csv diff --git a/RCode/superNll.R b/RCode/superNll.R new file mode 100755 index 0000000..f0f75d7 --- /dev/null +++ b/RCode/superNll.R @@ -0,0 +1,172 @@ +# Custom function to optimize parameters using +# negative log likelihood for anova and linear regression +# x := The x data we wish to observe a relationship for +# y := The y data we wish to observe a relationship for +# nLevels := number of levels to be used; Default is 2 for simple linear regression +# nExpUnits := Number of data points in the data set +# anova := Will we be using linear regression or anvoa? Default is false + +superNll = function(x,y,nLevels=2,nExpUnits, anova=FALSE){ + +# If we are using anova, use a more complex matrice solution for nll +if (anova == TRUE){ + # Define the function to optimize + nll=function(p,x,y){ + # Declare empty matrices to hold data + B = matrix(0, 1, nLevels) + expected = 0 + XX = matrix(0, nExpUnits, nLevels) + + # Fill in a matrix of size [nExpUnits x nLevels] + # This will hold binary values for whether or not + # a specific beta parameter is used + for (j in 1:nLevels){ + for (i in ((nExpUnits*j/nLevels)-((nExpUnits/nLevels)-1)):(nExpUnits*j/nLevels)){ + XX[i,j] = 1 + } + } + + #Define the values for the beta parameters + B[1] = p[1] + for (i in 2:nLevels){ + B[i] = p[i] + expected = expected + B[i]*XX[,i] + } + expected = expected + B[1] + + # Define the value for the residual standard error parameters + sigma=exp(p[nLevels+1]) + + # Compute and return the negative log likelihood + nll=-sum(dnorm(x=y,mean=expected,sd=sigma,log=TRUE)) + return(nll) + } + + # Initial guess for optimization + # Declare empty matrix to be filled in + initialGuess = matrix(0, 1, nLevels+1) + # First guess is the mean of the first subset of data + initialGuess[1] = mean(y[((nExpUnits/nLevels)-((nExpUnits/nLevels)-1)):(nExpUnits/nLevels)]) + # Initial guess for residual standard error is 1 + initialGuess[nLevels+1] = 1 + # Subsequent guesses after the first are differences in + # means of the current subset of the data compared to the first subset + for (i in 2:nLevels){ + initialGuess[i] = mean(y[((nExpUnits*i/nLevels)-((nExpUnits/nLevels)-1)):(i*nExpUnits/nLevels)])-initialGuess[1] + } + + # Optimize the nll function + myNll=optim(par=initialGuess,fn=nll,x=x,y=y, hessian = TRUE) + + # Simple model for nll function + nllSimple_lin<-function(p,x,y){ + B0=p[1] + sigma=exp(p[2]) + expected=B0 + nll=-sum(dnorm(x=y,mean=expected,sd=sigma,log=TRUE)) + return(nll) + } + + # Initial guess for simple model + simpleGuess=c(1,1) + + # Optimize the simple model + simpleFit=optim(par=simpleGuess,fn=nllSimple_lin,x=x,y=y) + + # Compute the LR statistic + teststat_lin=2*(simpleFit$value-myNll$value) + + # Compute the degrees of freedom + df_lin=length(myNll$par)-length(simpleFit$par) + + # P value from nll function + pValue=1-pchisq(teststat_lin,df_lin) + + # Coefficients from nll function + coefficients = myNll$par[1:nLevels] + + # Sigma from nll function + sigma = exp(myNll$par[nLevels+1]) + + # Return the pValue, coefficients, and sigma + nllReturns = NULL + attributes(nllReturns)$pValue = pValue + attributes(nllReturns)$coefficients = coefficients + attributes(nllReturns)$sigma = sigma + attributes(nllReturns)$hessian = myNll$hessian + nllReturns = attributes(nllReturns) + return(nllReturns) +} + +# If we are not using anova, use simple linear regression for nll +if (anova == FALSE){ + # Linear log likelihood function + nllLinear=function(p,x,y){ + B0=p[1] + B1=p[2] + sigma=exp(p[3]) + expected=B0+x*B1 + nll=-sum(dnorm(x=y,mean=expected,sd=sigma,log=TRUE)) + return(nll) + } + # Initial guess for optimization + # First guess tries to approximate the intercept from slope + # Second guess uses end points of the data to approximate the slope + # Guess of the residual standard error is 1 + slope = (y[nExpUnits]-y[1]) / (x[nExpUnits]-x[1]) + firstGuess = y[1]-slope*x[1] + initialGuess_lin=c(firstGuess, + slope, + 1) + # Optimize the nll function + linearFit=optim(par=initialGuess_lin,fn=nllLinear,x=x,y=y,hessian = TRUE) + + # Simple model for nll function + nllSimple_lin<-function(p,x,y){ + B0=p[1] + sigma=exp(p[2]) + expected=B0 + nll=-sum(dnorm(x=y,mean=expected,sd=sigma,log=TRUE)) + return(nll) + } + + # Initial guess for simple model + simpleGuess_lin=c(firstGuess,1) + + # Optimize the simple model + simpleFit=optim(par=simpleGuess_lin,fn=nllSimple_lin,x=x,y=y, hessian = TRUE) + + # Compute the LR statistic + teststat_lin=2*(simpleFit$value-linearFit$value) + + # Compute the degrees of freedom + df_lin=length(linearFit$par)-length(simpleFit$par) + + # P value from nll function + pValue=1-pchisq(teststat_lin,df_lin) + + # Coefficients from nll function + coefficients = linearFit$par[1:2] + + # Sigma from nll function + sigma = exp(linearFit$par[3]) + + # Return the pValue, coefficients, and sigma + nllReturns = NULL + attributes(nllReturns)$pValue = pValue + attributes(nllReturns)$coefficients = coefficients + attributes(nllReturns)$sigma = sigma + attributes(nllReturns)$hessian = linearFit$hessian + nllReturns = attributes(nllReturns) + return(nllReturns) +} + +} + + + + + + + + diff --git a/README.md b/README.md index 2854435..5b3579c 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,19 @@ -# biocomputing_StatsGroupProject +# BIOS 60318 Introduction To Biocomputing Final Project +## Statistical analysis of ANOVA versus linear regression using LSE and MLLE +### Francisco Huizar, Ian Kluper, Nate Hermann -This is a repo for the statiatics group project in BIOS30318/60318 at the -University of Notre Dame. +This is a repository for the FALL2018 BIOS 60318 final project. +Here we compare ANOVA versus linear regression using two different methods: +1. Least Squares Estimate (LSE) +2. Maximum Log-Likelihood Estimate (MLLE) -When completed the students' forked repos should contain: -1. A script generating ANOVA and regression results for the provided data sets. +This repository contains: +1. A script generating ANOVA and regression results for the antibiotics.csv and sugar.csv data sets [BIOS60318_FinalProject.R](https://github.com/fjhuizar/biocomputing_StatsGroupProject/blob/master/RCode/BIOS60318_FinalProject.R) +2. A document displaying figures, and interpretation of results can be found [here](https://github.com/fjhuizar/biocomputing_StatsGroupProject/blob/master/Final%20Writeup/BIOS60318-FinalProject-Writeup.pdf) +3. A script implementing power analysis to compare ANOVA and linear regression [BIOS60318_FinalProject.R](https://github.com/fjhuizar/biocomputing_StatsGroupProject/blob/master/RCode/BIOS60318_FinalProject.R) +4. A document with figures and text summarizing findings of the power analysis can be found [here](https://github.com/fjhuizar/biocomputing_StatsGroupProject/blob/master/Final%20Writeup/BIOS60318-FinalProject-Writeup.pdf) -2. A word document displaying figures and interpretationof results, including specification of the simple and complex models and likelihood ratio test results. - -3. A script implementing the power analysis. - -4. A word document with figures and text summarizing findings of the power analysis. This includes answers to the final question on the assignment. +To run the code yourself: +1. Fork this repository to your own, and clone that repository onto your computer +2. Open the RCodes folder, and run [BIOS60318_FinalProject.R](https://github.com/fjhuizar/biocomputing_StatsGroupProject/blob/master/RCode/BIOS60318_FinalProject.R) diff --git a/Simulation Results/10000_Simulation_Results.xlsx b/Simulation Results/10000_Simulation_Results.xlsx new file mode 100755 index 0000000..5a95f91 Binary files /dev/null and b/Simulation Results/10000_Simulation_Results.xlsx differ diff --git a/Simulation Results/1000_Simulation_Results.xlsx b/Simulation Results/1000_Simulation_Results.xlsx new file mode 100755 index 0000000..7bf149e Binary files /dev/null and b/Simulation Results/1000_Simulation_Results.xlsx differ