#script for data analysis.  follows "exploratory data analysis" and
#uses sample MO data with Num.Surveys<8 removed
#based on attached data summaryMO.csv and MOtest.txt

attach(MOcomb8)  #most recent iteration of MO sample data

#find the most meaningful scale for each variable
cor(richness, PROP_DEV_300, use="pairwise.complete.obs")
cor(richness, PROP_DEV_600, use="pairwise.complete.obs")
cor(richness, PROP_DEV_1000, use="pairwise.complete.obs")
cor(richness, PROP_DEV_5000, use="pairwise.complete.obs")
cor(richness, PROP_DEV_10000, use="pairwise.complete.obs")

#test whether these correlations are significant
cor.test(richness, PROP_DEV_300, use="pairwise.complete.obs")
cor.test(richness, PROP_DEV_600, use="pairwise.complete.obs")
cor.test(richness, PROP_DEV_1000, use="pairwise.complete.obs")
cor.test(richness, PROP_DEV_5000, use="pairwise.complete.obs")
cor.test(richness, PROP_DEV_10000, use="pairwise.complete.obs")

#forest.  
cor(richness, PROP_FOR_300, use="pairwise.complete.obs")
cor(richness, PROP_FOR_600, use="pairwise.complete.obs")
cor(richness, PROP_FOR_1000, use="pairwise.complete.obs")
cor(richness, PROP_FOR_5000, use="pairwise.complete.obs")
cor(richness, PROP_FOR_10000, use="pairwise.complete.obs")

#and row crops
cor(richness, PROP_AGR_300, use="pairwise.complete.obs")
cor(richness, PROP_AGR_600, use="pairwise.complete.obs")
cor(richness, PROP_AGR_1000, use="pairwise.complete.obs")
cor(richness, PROP_AGR_5000, use="pairwise.complete.obs")
cor(richness, PROP_AGR_10000, use="pairwise.complete.obs")

#wetland area
cor(richness, WET_AREA_300, use="pairwise.complete.obs")
cor(richness, WET_AREA_600, use="pairwise.complete.obs")
cor(richness, WET_AREA_1000, use="pairwise.complete.obs")
cor(richness, WET_AREA_5000, use="pairwise.complete.obs")
cor(richness, WET_AREA_10000, use="pairwise.complete.obs")

#total road length
cor(richness, T_ROAD_LEN_300, use="pairwise.complete.obs")
cor(richness, T_ROAD_LEN_600, use="pairwise.complete.obs")
cor(richness, T_ROAD_LEN_1000, use="pairwise.complete.obs")
cor(richness, T_ROAD_LEN_5000, use="pairwise.complete.obs")
cor(richness, T_ROAD_LEN_10000, use="pairwise.complete.obs")

#could similarly use AICs from the GLMs for any of these
dev300<-glm(richness ~ PROP_DEV_300, family="poisson", data=MOcomb8)
summary(dev300)

dev600<-glm(richness ~ PROP_DEV_600, family="poisson", data=MOcomb8)
summary(dev600_
        
dev1000<-glm(richness ~ PROP_DEV_1000, family="poisson", data=MOcomb8)
summary(dev1000)

dev5<-glm(richness ~ PROP_DEV_5000, family="poisson", data=MOcomb8)
summary(dev5)
        
dev10<-glm(richness ~ PROP_DEV_10000, family="poisson", data=MOcomb8)
summary(dev10)

#see what a model with all variables looks like
fullGLM<-glm(richness ~ PROP_DEV_300 + PROP_FOR_5000 + WET_AREA_300 + T_ROAD_LEN_300, family="poisson", data=MOcomb8)
summary(fullGLM)

#individual species results should be more interesting (less varied responses to land use)
#try to choose species at 0.3 to 0.7 of all sites
        

#comparing factors across scales for individual species        

        anfwGLM300<-glm(ANFW ~ T_ROAD_LEN_300, family="binomial", data=MOcomb8)
        summary(anfwGLM300)
        
        anfwGLM600<-glm(ANFW ~ T_ROAD_LEN_600, family="binomial", data=MOcomb8)
        summary(anfwGLM600)
        
        anfwGLM1000<-glm(ANFW ~ T_ROAD_LEN_1000, family="binomial", data=MOcomb8)
        summary(anfwGLM1000)
        
        anfwGLM5<-glm(ANFW ~ T_ROAD_LEN_5000, family="binomial", data=MOcomb8)
        summary(anfwGLM5)
        
        anfwGLM10<-glm(ANFW ~ T_ROAD_LEN_10000, family="binomial", data=MOcomb8)
        summary(anfwGLM10)
        
        #sample model for a single species

        licaGLM<-glm(LICA ~ T_ROAD_LEN_300 + PROP_DEV_300 + PROP_FOR_5000, family="binomial", data=MOcomb8)
        summary(licaGLM)


#create new variable that combines primary and secondar road length

MOcomb8$PS_ROAD_LEN_300<-MOcomb8$P_ROAD_LEN_300 + MOcomb8$S_ROAD_LEN_300
MOcomb8$PS_ROAD_LEN_600<-MOcomb8$P_ROAD_LEN_600 + MOcomb8$S_ROAD_LEN_600
MOcomb8$PS_ROAD_LEN_1000<-MOcomb8$P_ROAD_LEN_1000 + MOcomb8$S_ROAD_LEN_1000
MOcomb8$PS_ROAD_LEN_5000<-MOcomb8$P_ROAD_LEN_5000 + MOcomb8$S_ROAD_LEN_5000
MOcomb8$PS_ROAD_LEN_10000<-MOcomb8$P_ROAD_LEN_10000 + MOcomb8$S_ROAD_LEN_10000

attach(MOcomb8)  # need to reattach since the data was modified
#now look at correlation coefficients for PS versus O roads        
cor(richness, PS_ROAD_LEN_300, use="pairwise.complete.obs")
cor(richness, PS_ROAD_LEN_600, use="pairwise.complete.obs")
cor(richness, PS_ROAD_LEN_1000, use="pairwise.complete.obs")
cor(richness, PS_ROAD_LEN_5000, use="pairwise.complete.obs")
cor(richness, PS_ROAD_LEN_10000, use="pairwise.complete.obs")

cor(richness, O_ROAD_LEN_300, use="pairwise.complete.obs")
cor(richness, O_ROAD_LEN_600, use="pairwise.complete.obs")
cor(richness, O_ROAD_LEN_1000, use="pairwise.complete.obs")
cor(richness, O_ROAD_LEN_5000, use="pairwise.complete.obs")
cor(richness, O_ROAD_LEN_10000, use="pairwise.complete.obs")

#compare using AICs from GLMs for the same models
PSROAD<-glm(richness ~ PS_ROAD_LEN_300, family="poisson", data=MOcomb8)
summary(PSROAD)
        
OROAD<-glm(richness ~ O_ROAD_LEN_300, family="poisson", data=MOcomb8)
summary(OROAD)

#compare wetland connectivity to wetland area measures
cor(richness, WET_AREA_300, use="pairwise.complete.obs")
cor(richness, WET_AREA_600, use="pairwise.complete.obs")
cor(richness, WET_AREA_1000, use="pairwise.complete.obs")
cor(richness, WET_CONNECT_10000, use="pairwise.complete.obs")