Statstics

#Means Tests, ANOVA, and Regression in R library(dplyr) #Two weeks ago weexplored the results from 3 common statistical tests: Students T-tests, ANOVAs, and Linear Regression #Today we will revisit that exercise but add the actual statistical tests that we have learned in lecture #25% of your grade for assignment #5 will directly come from output that we generate in class #Another 25% of your grade for assignment $5 will involve manipulating these data, so... #Be sure to take detailed comments so you remember how to interpret what we have done! #We will be using the iris data set. #all measurements are in centimeters data(iris) summary(iris) head(iris) #------------PART ONE: ONE-SAMPLE T TESTS----------------- #Let's start with a one-sample means test #The overall mean petal length is 3.758 centimeters. #Is the mean petal length for any species statistically different than 3.758cm? #Let's have a look at species setosa test1=t.test(iris$Petal.Length[iris$Species=="setosa"],mu = 3.758) test1 #what does the output tell us? #can we re-create the test statistic by hand? t = (mean(iris$Petal.Length[iris$Species=="setosa"])-3.758)/(sd(iris$Petal.Length[iris$Species=="setosa"])/sqrt(length(iris$Petal.Length[iris$Species=="setosa"]))) t #and we can access the lookup table in R, too! p = pt(t,length(iris$Species[iris$Species=="setosa"]))*2 #the 2 is for a two-tailed test #does our code by hand match the t.test package? #some clunky code above...how can we make it easier to work with? #check the other species... test2=t.test(iris$Petal.Length[iris$Species=="versicolor"],mu = 3.758) test2 test3=t.test(iris$Petal.Length[iris$Species=="virginica"],mu = 3.758) test3 #How many of the species had mean petal lengths that were statistically different from the overall petal length? #Sample code for one-tailed one-sample test test2b=t.test(iris$Petal.Length[iris$Species=="versicolor"],alt="greater",mu = 3.758) test2b #compare to output from test2 test2 #where do you see differences in the output? #----------PART TWO: TWO-SAMPLE T TESTS----------------- #we will need to tell R how to calculate standard error se <- function(x) sd(x)/sqrt(length(x)) #Plot 1.Students T-test #Use of T-test: Compare the means between two groups unique(iris$Species) # we have three species of iris #Question: Does Petal Length differ between setosa and virginica? #Null hypothesis: Petal length does not differ between the two species #Step 1. Calculate Mean and Standard Error for Data PetalLength_ds = group_by(iris, Species) %>% #cluster by species summarize(mean=mean(Petal.Length), #mean petal length se=se(Petal.Length), n=length(Petal.Length)) PetalLength_ttest= PetalLength_ds[-2,] #removing versicolor to plot t test #Step 2. Plot Mean ttest_plot=barplot(PetalLength_ttest$mean,names.arg=c("setosa","virginica"), col="green",ylim=c(0,6),ylab="Average Petal Length",xlab="Species") #basic plot arrows(ttest_plot, PetalLength_ttest$mean-PetalLength_ttest$se, ttest_plot, PetalLength_ttest$mean+PetalLength_ttest$se, angle = 90,code = 3) #So, we see on the graphs that there are large differences between the group means #Our visual clues suggest that these two groups have statistically different mean petal lengths #Do they? Let's run a two-sample t-test to find out... #and let's make our life a little easier... sample1=iris$Petal.Length[iris$Species=="setosa"] sample2=iris$Petal.Length[iris$Species=="virginica"] sample3=iris$Petal.Length[iris$Species=="versicolor"] test4=t.test(sample1,sample2) test4 #by default, the syntax is for a two-tailed test #let's try a one-tailed test...is setosa significantly smaller than virginica? test5=t.test(sample1,sample2,alt="less") test5 test4$p.value test5$p.value #how do the p-values compare between the one-tailed and two-tailed tests? #let's compare the two species with closer means... test6=t.test(sample2,sample3,alt="greater") test6 #what happens if we only look at a fraction of the samples? test7=t.test(sample2[1:7],sample3[1:7],alt="greater") test7 #---------PART THREE: ANOVA-------------- #Plot 2.ANVOA #Use of ANOVA: Compare the means between MORE THAN two groups #Question: Does Petal Length differ between all the iris species? #Null hypothesis: Petal length does not differ between the species #Step 1. Calculate Mean and Standard Error for Data PetalLength_ds #Step 2. Plot Mean ANOVA_plot=barplot(PetalLength_ds$mean,names.arg=c("setosa","versicolor","virginica"), col="green",ylim=c(0,6),ylab="Average Petal Length",xlab="Species") #basic plot arrows(ANOVA_plot, PetalLength_ds$mean-PetalLength_ds$se, ANOVA_plot, PetalLength_ds$mean+PetalLength_ds$se, angle = 90,code = 3) #Plot 3. Looking at the standard Error #Question: What happens when we reduce the sample size to 10 observations per species, what about 5 observations per species? #Null hypothesis: Petal length does not differ between the species #PetalLength_Error = # group_by(iris, Species) %>% #cluster by species # summarize(mean=mean(Petal.Length), #mean petal length # se=sd(Petal.Length)/sqrt(5), # n=length(Petal.Length)) #Error_plot=barplot(PetalLength_Error$mean,names.arg=c("setosa","versicolor","virginica"), col="green",ylim=c(0,6),ylab="Average Petal Length",xlab="Species") #basic plot #arrows(Error_plot, PetalLength_Error$mean-PetalLength_Error$se, Error_plot, PetalLength_Error$mean+PetalLength_Error$se, angle = 90,code = 3) #This plot also points to a statistically significant difference between the groups...is there? test8=aov(iris$Petal.Length~iris$Species) summary(test8) sum_test=unlist(summary(test8)) sum_test["Pr(>F)1"] #what does the aov output tell us? #can we calculate the F value by hand? #numerator = the variance of the sample means times the number of observations per sample numer = var(c(mean(sample1),mean(sample2),mean(sample3))) * length(sample1) #denominator = the mean of the sample variances (for balanced design) denom = mean(c(var(sample1),var(sample2),var(sample3))) ourF = numer/demon ourF #let's look up the p-value for our test statistic #numerator df = number of groups-2 #denominator df = number of groups * n-1 ourP = pf(ourF,2,3*(length(sample1)-1),lower.tail=FALSE) ourP #are the differences between the other three variables (sepal length, sepal width, pedal width) also statistically different? #before you leave, read in the class survey data and try to run one of the tests we have done today. #be sure you can interpret the output! #possible example: compare the mean spiciness preference of Taco Bell vs. Chipotle fans... #note: our class data do not necessarily meet all of the assumptions necessary for all of these tests #try to recognize when the data are especially skewed or otherwise problematic, which could impact the validity of your results #We will pick up with regression on Monday April 16. Good luck on the assignment! #Plot 4. Regression #The relationship between two continous variables (strength and direction) #Question: Are sepal and petal width related? If so, are they positively or negatively related? #Null hypothesis: Petal Width and Petal Length are not related plot(iris$Petal.Width, iris$Petal.Length, pch = 16, cex = 1.3, col = "black", xlab = "Petal Width", ylab = "Petal Length") abline(lm(iris$Petal.Length~iris$Petal.Width),lwd=3,col="red")