Double stdev plot

plotting

require(ggplot2)
require(stringr)

halo = read.table('Halo-Size.txt', header= T)
dry = read.table('Dry-Mass.txt', header =T)



# calculation of standard deviation
haloError = qnorm(0.975) * apply(halo, 2, sd)/sqrt(apply(halo, 2, length))
dryError = qnorm(0.975) * apply(dry, 2, sd)/sqrt(apply(dry, 2, length))

# means
haloMean = apply(halo, 2, mean)
dryMean = apply(dry, 2, mean)


# build frame to plot
frame = data.frame(strain = factor(str_extract(names(haloMean),"(?<=X).*"), 
                                   levels = str_extract(names(haloMean),"(?<=X).*")),
                   halo = haloMean,
                   dry = dryMean,
                   haloPos = haloMean + haloError,
                   haloNeg = haloMean - haloError,
                   dryPos = dryMean + dryError,
                   dryNeg = dryMean - dryError)

# plot
p = (ggplot(frame, aes(x = halo, y = dry,color = strain)) 
     + geom_point()
     + geom_errorbar(aes(ymin=dryNeg, ymax= dryPos)) 
     + geom_errorbarh(aes(xmin =haloNeg, xmax = haloPos)))


plot(p)

significance

library(reshape2)

halo = read.table('Halo-Size.txt', header= T)
dry = read.table('Dry-Mass.txt', header =T)

# anova ---------------

halo = melt(halo)
names(halo) = c('colony', 'value')
dry = melt(dry)
names(dry) = c('colony', 'value')

# you want to read the Pr. If the values were matched you could have done a 
# single anova to get signicance for both but alas you cannot
summary(aov(value ~ colony , halo))
summary(aov(value ~ colony , dry))

# this only tells you if the variable (groups) has different distributions than
# each other. for both dimensions, this is true

# comparison between groups -----------------
# you should pick your comparisons beforehand if you can since the 
# significance will drop with the number of comparisons you have. If you are 
# trying to compare one group to all the rest, that is the way to go


# creating all comparisons

groups = unique(halo$colony)
comb = combn(groups,2)

# get p values for each comparison
# using for isn't the nice way to do it in R but it's easier to read
ps = c()
for (i in 1:ncol(comb)){
    # this part depends if you expect your distribution to be normal or not.
    # if it is normal use a t.test if not wilcox.test is a better option. but 
    # your sample size isn't enough to give significance with wilcox
    group1 = halo$value[halo$colony == comb[1,i]]
    group2 = halo$value[halo$colony == comb[2,i]]
    p = t.test(group1,group2)$p.value
    ps = c(ps,p)
}

ps= p.adjust(ps, 'fdr')

# prints significant ones
comb[,ps<0.05]

# for mass -----------
groups = unique(dry$colony)
comb = combn(groups,2)

# get p values for each comparison
# using for isn't the nice way to do it in R but it's easier to read
ps = c()
for (i in 1:ncol(comb)){
    # this part depends if you expect your distribution to be normal or not.
    # if it is normal use a t.test if not wilcox.test is a better option. but 
    # your sample size isn't enough to give significance with wilcox
    group1 = dry$value[dry$colony == comb[1,i]]
    group2 = dry$value[dry$colony == comb[2,i]]
    p = t.test(group1,group2)$p.value
    ps = c(ps,p)
}

ps= p.adjust(ps, 'fdr')

# prints significant ones
comb[,ps<0.05]