# read 2007 file ##### pp <- read.csv("flPhenology.xy.C1.2007.csv") # add nearest neighbor distances for nn 1 to 40 #### library(FNN) kk <- knn.dist(pp[, c("x", "y")], k = 40) kk <- as.data.frame(kk) # str(kk) names(kk) <- paste("nn", 1:40, sep= "") pp <- cbind(pp, kk) # add core and edge designation #### pp$core <- "core" mx <- mean(pp$x) my <- mean(pp$y) pp$dist2Core <- (mx - pp$x)^2 + (my - pp$y)^2 median(pp$dist2Core) pp[pp$dist2Core > median(pp$dist2Core), "core"] <- "edge" pp$core <- as.factor(pp$core) plot(y ~ x, asp = 1, data = pp, pch = 19, col = gray(as.numeric(pp$core)/2.5)) # seed set for models #### pp$ss <- pp$full/(pp$empty + pp$full) pp$tot <- pp$empty + pp$full # delete some records #### dim(pp) # 529 table(pp$full, useNA="always") with(pp[is.na(pp$full), ], points(x, y, pch = ".", col = "red")) pp <- pp[!is.na(pp$avpeak_jd),] # records missing dates pp <- pp[!is.na(pp$full),] # reconds missing seed set dim(pp) table(pp$full, useNA="always") # seed set for glms #### sf4 <- cbind(pp$full,pp$empty) # str(sf4) # head(sf4) # dim(pp) # model selex for 2007 #### # try glm m1 <- glm(sf4 ~ core * avpeak_jd * log(nn12), data = pp, binomial) m2 <- glm(sf4 ~ core * avpeak_jd + log(nn12) * core + log(nn12) * avpeak_jd, data = pp, binomial) anova(m2, m1, test = "Chi") # **, but summary(m1) # dispersion # try quasibinomial m1 <- glm(sf4 ~ core * avpeak_jd * log(nn12), data = pp, quasibinomial) m2 <- glm(sf4 ~ core * avpeak_jd + log(nn12) * core + log(nn12) * avpeak_jd, data = pp, quasibinomial) anova(m2, m1, test = "F") # eliminate 3-way (m2 is best for now) m3 <- glm(sf4 ~ log(nn12) * core + log(nn12) * avpeak_jd, data = pp, quasibinomial) m4 <- glm(sf4 ~ core * avpeak_jd + log(nn12) * avpeak_jd, data = pp, quasibinomial) m5 <- glm(sf4 ~ core * avpeak_jd + log(nn12) * core , data = pp, quasibinomial) anova(m3, m2, test = "F") # no evidence for core X avpeak, p = 0.33 anova(m4, m2, test = "F") # p = 0.48 anova(m5, m2, test = "F") # p = 0.74, eliminate nn X avpeak (m5 is best for now) m6 <- glm(sf4 ~ core : avpeak_jd + log(nn12) + core + avpeak_jd, data = pp, quasibinomial) m7 <- glm(sf4 ~ log(nn12) : core + log(nn12) + core + avpeak_jd, data = pp, quasibinomial) anova(m6, m5, test = "F") # p = 0.46, eliminate nn X core anova(m7, m5, test = "F") # p = 0.35 elim core X avpeak (m6 is best for now) m8 <- glm(sf4 ~ log(nn12) + core + avpeak_jd, data = pp, quasibinomial) anova(m8, m6, test = "F") # 0.39, elim core X avpeak m101 <- glm(sf4 ~ core + avpeak_jd + log(nn12), data = pp, quasibinomial) m102 <- glm(sf4 ~ core + avpeak_jd , data = pp, quasibinomial) m103 <- glm(sf4 ~ core + log(nn12), data = pp, quasibinomial) m104 <- glm(sf4 ~ avpeak_jd + log(nn12), data = pp, quasibinomial) anova(m102, m101, test = "F") # p = 0.132 # do parametric bs on this test anova(m103, m101, test = "F") # highly signif anova(m104, m101, test = "F") # highly signif # best model is m102, little support for m101 m13 <- glm(sf4 ~ avpeak_jd, data = pp, quasibinomial) m14 <- glm(sf4 ~ core , data = pp, quasibinomial) m15 <- glm(sf4 ~ 1 , data = pp, quasibinomial) anova(m13, m102, test = "F") # highly signif anova(m14, m102, test = "F") # highly signif m101 # best model plus space m102 # best model par(mfcol = c(2,2)) plot(m102) # qq plot doesn't look great, but residuals are OK par(mfcol = c(1,1)) # save the data frame and models m101 & m6 for making figures #### pp.2007 <- pp mm.2007a <- m102 mm.2007b <- m101 mm.2007c <- m6 # save(pp.2007, mm.2007a, mm.2007b, mm.2007c, file = "forFig2007new1.RData") # use these results for the table #### anova(m1, m2, m5, m6, m101, m102, m13, test = "F") anova(m102, m14, test = "F")