# statistical analysis of Springcasting data
# T Crimmins (theresa@usanpn.org) 1-8-21

library(car)

# load files
leafdata <- read.csv(file = "C:/NPN/Manuscripts/Working/Springcasting_evaluation/SpringcastingAnalyses/2018-2020_lilac_leaf_merge.csv")

bloomdata <- read.csv(file = "C:/NPN/Manuscripts/Working/Springcasting_evaluation/SpringcastingAnalyses/2018-2020_lilac_bloom_merge.csv")

# calculate diff between DOY leaf-out observed & leaf-out reported (predicted - observed)

leafdata$PredObs <- (leafdata$SI.x_Leaf_Value - leafdata$first_yes_doy)

bloomdata$PredObs <- (bloomdata$SI.x_Bloom_Value - bloomdata$first_yes_doy)

# filter out records w/NAs for PredObs
leafdata2 <- leafdata[!is.na(leafdata$PredObs), ]

bloomdata2 <- leafdata[!is.na(bloomdata$PredObs), ]

# plot histogram
hist(leafdata2$PredObs)

hist(bloomdata2$PredObs)

#drop outliers
leafdata2 <- leafdata2[-c(244),]

bloomdata2 <- bloomdata[-c(439),]

# ANOVAs on accuracy
summary(aov(leafdata2$PredObs ~ leafdata2$Code))

leaf.av <- (aov(leafdata2$PredObs ~ leafdata2$Code))

summary(aov(bloomdata2$PredObs ~ bloomdata2$Code))

bloom.av <- (aov(bloomdata$PredObs ~ bloomdata$Code))


# box & whisker plot
boxplot(PredObs~Code,data=leafdata2)

boxplot(PredObs~Code,data=bloomdata2)



# multiple comparisons

tukey.test <- TukeyHSD(bloom.av)
tukey.test

tukey.test2 <- TukeyHSD(bloom.av)
tukey.test2

#Levine's test - use to test precision
leveneTest(PredObs ~ Code, data = leafdata2)

leveneTest(PredObs ~ Code, data = bloomdata2)