Showing posts with label statistics. Show all posts
Showing posts with label statistics. Show all posts

Wednesday, November 26, 2014

11 24 2014 Non Parametric Test for Differences in Size

y1boxplot.R
require(ggplot2)
## Loading required package: ggplot2
require(plyr)
## Loading required package: plyr
require(splitstackshape)
## Loading required package: splitstackshape
## Loading required package: data.table
require(nparcomp)
## Loading required package: nparcomp
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: splines
## Loading required package: TH.data
y1size=read.csv('Y1size.csv')
#creates dataframe and reads in the CSV file for sizes
View(y1size)
#check data
y1size$Date<-as.Date(y1size$Date, "%m/%d/%Y")
#make R understand dates
y1meansize<-ddply(y1size,.(Date,Site,Pop),summarise, mean_size=mean(Length.mm,na.rm=T))
#create table of ave size for outplant and year one for each pop at each site
#print it out
print(y1meansize)
##          Date       Site Pop mean_size
## 1  2013-08-16    Fidalgo  2H     10.67
## 2  2013-08-16    Fidalgo  2N     11.60
## 3  2013-08-16    Fidalgo  2S     11.25
## 4  2013-08-16 Manchester  4H     10.53
## 5  2013-08-16 Manchester  4N     13.40
## 6  2013-08-16 Manchester  4S     11.30
## 7  2013-08-16 Oyster Bay  1H     10.49
## 8  2013-08-16 Oyster Bay  1N     10.90
## 9  2013-08-16 Oyster Bay  1S     12.15
## 10 2014-09-19 Oyster Bay  1H     27.96
## 11 2014-09-19 Oyster Bay  1N     34.65
## 12 2014-09-19 Oyster Bay  1S     27.98
## 13 2014-10-17    Fidalgo  2H     24.40
## 14 2014-10-17    Fidalgo  2N     29.10
## 15 2014-10-17    Fidalgo  2S     28.91
## 16 2014-10-24 Manchester  4H     21.49
## 17 2014-10-24 Manchester  4N     24.37
## 18 2014-10-24 Manchester  4S     23.99
#now we need to create subsets for each site for out plant and end of year 1
outmany1<-ddply(y1size,.(Length.mm,Pop,Tray,Sample,Area),subset,Date=="2013-08-16"&Site=="Manchester")
outfidy1<-ddply(y1size,.(Length.mm,Pop,Tray,Sample,Area),subset,Date=="2013-08-16"&Site=="Fidalgo")
outoysy1<-ddply(y1size,.(Length.mm,Pop,Tray,Sample,Area),subset,Date=="2013-08-16"&Site=="Oyster Bay")
endmany1<-ddply(y1size,.(Length.mm,Pop,Tray,Sample,Area),subset,Date=="2014-10-24"&Site=="Manchester")
endfidy1<-ddply(y1size,.(Length.mm,Pop,Tray,Sample,Area),subset,Date=="2014-10-17"&Site=="Fidalgo")
endoysy1<-ddply(y1size,.(Length.mm,Pop,Tray,Sample,Area),subset,Date=="2014-09-19"&Site=="Oyster Bay")

ggplot()+
  geom_boxplot(data=outmany1,aes(x=Pop,y=Length.mm,fill=Pop))+
  scale_colour_manual(values=c("blue","purple","orange"))+
  scale_fill_manual(values=c("blue","purple","orange"))

plot of chunk unnamed-chunk-1

ggplot()+
  geom_boxplot(data=endmany1,aes(x=Pop,y=Length.mm,fill=Pop))+
  scale_colour_manual(values=c("blue","purple","orange"),guide=F)+
  scale_fill_manual(values=c("blue","purple","orange"), guide=F)+
  ylim(c(0,50))+
  labs(x="Population",y="Average Length (mm)")+
  scale_x_discrete(labels=c("Dabob","Fidalgo","Oyster Bay"))+
  annotate("text", x=c("4N","4H","4S"),y=50, label=c("A","B","A"),size=10)+
  theme_bw()+
  theme(axis.text.x=element_text(size=20),
        axis.title.x=element_text(size=25),
        axis.title.y=element_text(size=25, vjust=2),
        axis.text.y=element_text(size=20))

plot of chunk unnamed-chunk-1

ggplot()+
  geom_boxplot(data=outfidy1,aes(x=Pop,y=Length.mm,fill=Pop))+
  scale_colour_manual(values=c("blue","purple","orange"))+
  scale_fill_manual(values=c("blue","purple","orange"))

plot of chunk unnamed-chunk-1

ggplot()+
  geom_boxplot(data=endfidy1,aes(x=Pop,y=Length.mm,fill=Pop))+
  scale_colour_manual(values=c("blue","purple","orange"),guide=F)+
  scale_fill_manual(values=c("blue","purple","orange"),guide=F)+
  ylim(c(0,50))+
  labs(x="Population",y="Average Length (mm)")+
  scale_x_discrete(labels=c("Dabob","Fidalgo","Oyster Bay"))+
  annotate("text", x=c("2N","2H","2S"),y=50, label=c("A","B","A"),size=10)+
  theme_bw()+
  theme(axis.text.x=element_text(size=20),
        axis.title.x=element_text(size=25),
        axis.title.y=element_text(size=25, vjust=2),
        axis.text.y=element_text(size=20))
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

plot of chunk unnamed-chunk-1

ggplot()+
  geom_boxplot(data=outoysy1,aes(x=Pop,y=Length.mm,fill=Pop))+
  scale_colour_manual(values=c("blue","purple","orange"))+
  scale_fill_manual(values=c("blue","purple","orange"))

plot of chunk unnamed-chunk-1

ggplot()+
  geom_boxplot(data=endoysy1,aes(x=Pop,y=Length.mm,fill=Pop))+
  scale_colour_manual(values=c("blue","purple","orange"),guide=F)+
  scale_fill_manual(values=c("blue","purple","orange"),guide=F)+
  ylim(c(0,50))+
  labs(x="Population",y="Average Length (mm)")+
  scale_x_discrete(labels=c("Dabob","Fidalgo","Oyster Bay"))+
  annotate("text", x=c("1N","1H","1S"),y=50, label=c("B","A","A"),size=10)+
  theme_bw()+
  theme(axis.text.x=element_text(size=20),
        axis.title.x=element_text(size=25),
        axis.title.y=element_text(size=25, vjust=2),
        axis.text.y=element_text(size=20))
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

plot of chunk unnamed-chunk-1

normality<-ddply(y1size,.(Date,Site,Pop),summarize,n=length(Length.mm),sw=shapiro.test(as.numeric(Length.mm))[2])
y1size$Pop2<-y1size$Pop
y1size$Pop2<-revalue(y1size$Pop2,c("1H"="H","2H"="H","4H"="H","1N"="N","2N"="N","4N"="N","1S"="S","2S"="S","4S"="S"))
#Here we subset the data set to only include data from the end of year 1
endy1<-ddply(y1size,.(Length.mm,Site,Pop,Tray,Sample,Area,Pop2),subset,Date>="2014-09-19")
normality<-ddply(endy1,.(Date,Site,Pop),summarize,n=length(Length.mm),sw=shapiro.test(as.numeric(Length.mm))[2])
endy1$log<-log2(endy1$Length.mm)
endy1$asin<-asin(sign(endy1$Length.mm)*sqrt(abs(endy1$Length.mm)))
## Warning: NaNs produced
normality<-ddply(endy1,.(Date,Site,Pop),summarize,n=length(log),sw=shapiro.test(as.numeric(log))[2])

sizekw<-kruskal.test(endy1$Length.mm~endy1$Site,endy1)
print(sizekw)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  endy1$Length.mm by endy1$Site
## Kruskal-Wallis chi-squared = 426.2, df = 2, p-value < 2.2e-16
sizekwpop<-kruskal.test(endy1$Length.mm~endy1$Pop2,endy1)
print(sizekwpop)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  endy1$Length.mm by endy1$Pop2
## Kruskal-Wallis chi-squared = 230, df = 2, p-value < 2.2e-16
require(PMCMR)
## Loading required package: PMCMR
sizenemenyi1<-posthoc.kruskal.nemenyi.test(x=endy1$Length.mm,g=endy1$Site, method="Tukey")
## Warning: Ties are present, p-values are not corrected.
sizenemenyi1
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  endy1$Length.mm and endy1$Site 
## 
##            Fidalgo Manchester
## Manchester < 2e-16 -         
## Oyster Bay 1.1e-12 < 2e-16   
## 
## P value adjustment method: none
sizenemenyi2<-posthoc.kruskal.nemenyi.test(x=endy1$Length.mm,g=endy1$Pop2, method="Tukey")
## Warning: Ties are present, p-values are not corrected.
sizenemenyi2
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  endy1$Length.mm and endy1$Pop2 
## 
##   H       N      
## N < 2e-16 -      
## S 3.4e-14 4.9e-08
## 
## P value adjustment method: none
sizenemenyi3<-posthoc.kruskal.nemenyi.test(x=endy1$Length.mm,g=endy1$Site:endy1$Pop2, method="Tukey")
## Warning: Ties are present, p-values are not corrected.
sizenemenyi3
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  endy1$Length.mm and endy1$Site:endy1$Pop2 
## 
##              Fidalgo:H Fidalgo:N Fidalgo:S Manchester:H Manchester:N
## Fidalgo:N    < 2e-16   -         -         -            -           
## Fidalgo:S    7.8e-14   0.9995    -         -            -           
## Manchester:H 2.7e-07   < 2e-16   < 2e-16   -            -           
## Manchester:N 1.0000    9.2e-14   1.1e-13   7.3e-06      -           
## Manchester:S 0.9786    < 2e-16   8.5e-14   0.0004       0.9880      
## Oyster Bay:H 3.1e-10   0.2781    0.6352    1.2e-14      2.4e-08     
## Oyster Bay:N < 2e-16   1.5e-11   7.2e-13   < 2e-16      < 2e-16     
## Oyster Bay:S 1.8e-09   0.6813    0.9255    9.1e-14      5.7e-08     
##              Manchester:S Oyster Bay:H Oyster Bay:N
## Fidalgo:N    -            -            -           
## Fidalgo:S    -            -            -           
## Manchester:H -            -            -           
## Manchester:N -            -            -           
## Manchester:S -            -            -           
## Oyster Bay:H 9.5e-12      -            -           
## Oyster Bay:N < 2e-16      1.0e-13      -           
## Oyster Bay:S 6.0e-11      1.0000       6.9e-12     
## 
## P value adjustment method: none