I tried to find differences between delta G values in amino acids between 9 species. Here is a link to the table that I read in. (The header lists the species): https://www.dropbox.com/s/5l1kxytx94hhlsk/SPU_002551.txt?dl=0
In my program, I cut selected the highest and lowest 25% of the data to compare, though I am not sure which Anova to use. The repeated measures Anova assuming Sphericity seems to require 3 columns of data though I only have 2 and cannot add a 3rd position column because everything has been mixed up by selecting parts of the data. I cannot use the simple One-way Anova either because the variance at 25% of the data is not constant.
Is there any other type of Anova I could use to check if the delta G values between species are different while also having Greenhouse-Geisser and Huynh-Feldt Corrections?
Example of nonsensical output from Univariate Type III Repeated-Measures ANOVA Assuming Sphericity:
SS num Df Error SS den Df F Pr(>F)
(Intercept) 9.6914e+10 1 2656941741 88 3209.875 < 2.2e-16 ***
design 5.0190e+06 8 15829627 704 27.901 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
design 0.0022698 4.6815e-87
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
design 0.4192 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
HF eps Pr(>F[HF])
design 0.4378048 6.190704e-18
These p-values are too small for the slight differences I am trying to detect
Code that produces the nonsensical result:
options(contrasts=c("contr.sum","contr.poly"))
#adhere to the sum-to-zero convention for effect weights
#--------------------------------------------------------------------------------------------
#Getting input from anova.py
#--------------------------------------------------------------------------------------------
setwd("/Users/antonysagayaraj/Desktop/Anova")
info <- read.table("infoForRScript.txt", header=FALSE, sep = ',')
type <- toString(info[1,1])
cutLength <- as.numeric(info[1,2])
geneName <- toString(info[1,3])
setwd("/Users/antonysagayaraj/Desktop/Anova/CompiledText")
urchinframe <- read.table(paste(geneName,".txt",sep=""), header=TRUE, sep = ',')
#-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
#MAIN (A) Array
#-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
#-------------------------------------------------------------------------------------------
#Extracting columns as numeric + sorting
#-------------------------------------------------------------------------------------------
sortAscending <- function(SpeciesArray) { #sorts the dG values from the bottom up
return(sort(SpeciesArray, decreasing = FALSE, na.last = NA))
}
sortDescending <- function(SpeciesArray) { #sorts the dG values from the top down
return(sort(SpeciesArray, decreasing = TRUE, na.last = NA))
}
cutAscending <- function(SpeciesArray,cutLength2) { #cuts the highest dG values
temp <- sortAscending(SpeciesArray)
return(temp[1:cutLength2])
}
cutDescending <- function(SpeciesArray,cutLength2) { #cuts the lowerst dG values
temp <- sortDescending(SpeciesArray)
return(temp[1:cutLength2])
}
Afragilis1 <- urchinframe[,c("A..fra")]
#print(Afragilis1)
Hpulcherrimus2 <- urchinframe[,c("H..pul")]
#print(Hpulcherrimus2)
Pdepressus3 <- urchinframe[,c("P..dep")]
#print(Pdepressus3)
Sdroebachiensis4 <- urchinframe[,c("S..dro")]
#print(Sdroebachiensis4)
Sfranciscanus5 <- urchinframe[,c("S..fran")]
#print(Sfranciscanus5)
Sintermedius6 <- urchinframe[,c("S..int")]
#print(Sintermedius6)
Snudus7 <- urchinframe[,c("S..nud")]
#print(Snudus7)
Spallidus8 <- urchinframe [,c("S..pal")]
#print(Spallidus8)
Spurpuratus9 <- urchinframe [,c("S..prp")]
#print(Spurpuratus9)
if (type == "percentage") { #makes sure that the percentage is converted into an amino acid number to be cut
cutLength <- cutLength*length(Afragilis1)
}
#ANOVA FOR ASCENDING----------------------------------------------------------------------------------------Uses the top part of the data
dGlength <- length(cutAscending(Afragilis1,cutLength)) #gets the length of the cut data
ascendingaovurchin = data.frame(
c(rep("Afra",dGlength),rep("Hpul",dGlength),rep("Pdep",dGlength),rep("Sdro",dGlength),rep("Sfra",dGlength),rep("Sint",dGlength),rep("Snud",dGlength),rep("Spal",dGlength),rep("Sprp",dGlength)))
# makes a dataframe for use in the anova
colnames(ascendingaovurchin)= c("Species")
speciesList = list(cutAscending(Afragilis1,cutLength),cutAscending(Hpulcherrimus2,cutLength),cutAscending(Pdepressus3,cutLength),cutAscending(Sdroebachiensis4,cutLength),cutAscending(Sfranciscanus5,cutLength),cutAscending(Sintermedius6,cutLength),cutAscending(Snudus7,cutLength),cutAscending(Spallidus8,cutLength),cutAscending(Spurpuratus9,cutLength)) #list of sorted and cut dG values
for (x in 1:9) { #Goes down the species list and copies all the sorted and cut delta G values into the 2nd column of the data frame
a <- x
for (n in 1:(length(Sfranciscanus5))) {
ascendingaovurchin[n + ((a-1)*dGlength),2] <- speciesList[[a]][n]
}
}
matrix <- with(ascendingaovurchin, cbind(V2[Species=="Afra"], V2[Species=="Hpul"], V2[Species=="Pdep"], V2[Species=="Sdro"], V2[Species=="Sfra"], V2[Species=="Sint"], V2[Species=="Snud"], V2[Species=="Spal"], V2[Species=="Sprp"]))
model <- lm(matrix ~ 1)
design <- factor(c("Afra", "Hpul", "Pdep", "Sdro", "Sfra", "Sint", "Snud", "Spal", "Sprp"))
library(car)
aov <- Anova(model, idata=data.frame(design), idesign=~design, type="III")
summary(aov, multivariate=F)