#Octopus analysis #Import Data setwd("~/Documents/R files Octopus") Mating <- read.csv("~/Documents/R files Octopus/Mating.csv") ##Set factors Mating$maleorder=factor(Mating$maleorder) Mating$size=factor(Mating$size) Mating$m2msize=factor(Mating$m2msize) Mating$m2fsize=factor(Mating$m2fsize) Mating$f2fsize=factor(Mating$f2fsize) #Remove All Wild Males MatingNoNA=subset(Mating,maleorder!=0) MatingNoNA$maleorder=factor(MatingNoNA$maleorder,levels=c(1,2,3)) ##Linear model testing mod=lm(eggsP95~maleorder+matingtime+archnpump+sizeg+remSP,data=MatingNoNA) summary(mod) step(mod) plot(mod) plot(eggsP95NEW~maleorder, data=Mating) plot(eggsP95NEW~sizeg, data=Mating) plot(eggsP95NEW~matingtime, data=Mating) mod3=lm(eggsP95NEW~maleorder*sizeg,data=Mating) summary(mod3) #Percent eggs v order test mod5<-lm(eggsP95NEW~maleorder, data=Mating) summary(mod5) #Leaps leaps=regsubsets(eggsP95~maleorder+matingtime+archnpump+sizeg+remSP,data=MatingNoNA) summary(leaps) plot(leaps, scale="adjr2") plot(leaps, scale="bic") ls=summary(leaps) plot(ls$bic) leaps1=regsubsets(eggsP95~maleorder*matingtime*archnpump*sizeg,data=Mating) summary(leaps1) plot(leaps1, scale="bic") #Exploration # More Graphs plot(eggsP95NEW~maleorder, data=Mating, main="Percentage of Eggs Sired Vs. Male Order", xlab="Male Order", ylab="Percentage", ylim=c(0.0,1.0)) plot(eggsP95NEW~maleorder, data=Mating, main="Percentage of Eggs Sired Vs. Male Order", xlab="Male Order", ylab="Percentage") text(c(0.3, 0.9, -10, -10), "*", cex=2) mod7<-lm(eggsP95NEW~maleorder, data=Mating) summary(mod7) #Size Vs Prop 95 Eggs, with regression line plot(Mating$sizeg, Mating$eggsP95, main="Percentage of Eggs Sired Vs. Male Size", xlab="Male Size in Grams", ylab="Percentage") abline((lm(Mating$eggsP95~Mating$sizeg)), col="red") cor.test(Mating$sizeg, Mating$eggs95) #New Size vs Prop 95NEW, With ylim set plot(Mating$sizeg, Mating$eggsP95NEW, main="Percentage of Eggs Sired Vs. Male Size", xlab="Male Size in Grams", ylab="Percentage", ylim=c(0.0,1.0)) abline((lm(Mating$eggsP95NEW~Mating$sizeg)), col="red") cor.test(Mating$sizeg, Mating$eggs95NEW) #REMOVE NONA plot(MatingNoNA$sizeg, MatingNoNA$eggsP95NEW, main="Percentage of Eggs Sired Vs. Male Size", xlab="Male Size in Grams", ylab="Percentage") abline((lm(MatingNoNA$eggsP95NEW~MatingNoNA$sizeg)), col="red") cor.test(MatingNoNA$eggsP95NEW, MatingNoNA$sizeg) #Mating time vs. prop 95 eggs plot(MatingNoNA$matingtime, MatingNoNA$eggsP95NEW, main="Percentage of Eggs Sired Vs. Mating Time", xlab="Mating Time in Seconds", ylab="Percentage") abline((lm(MatingNoNA$eggsP95NEW~MatingNoNA$matingtime)), col="red") cor.test(MatingNoNA$eggsP95NEW, MatingNoNA$matingtime) #Eggs P vs. archnpump plot(MatingNoNA$archnpump, MatingNoNA$eggsP95NEW, main="Percentage of Eggs Sired Vs. Mating Time", xlab="Number of Arch and Pumps", ylab="Percentage") abline((lm(MatingNoNA$eggsP95NEW~MatingNoNA$archnpump)), col="red") cor.test(MatingNoNA$eggsP95NEW, MatingNoNA$archnpump) #Eggs Vs Archnpump with ylim plot(Mating$archnpump, Mating$eggsP95NEW, main="Percentage of Eggs Sired Vs. Number of Arch and Pumps", xlab="Number of Arch and Pumps", ylab="Percentage", ylim=c(0.0,1.0)) abline((lm(Mating$eggsP95NEW~Mating$archnpump)), col="red") cor.test(Mating$eggsP95NEW, Mating$archnpump) #Subset Regression library(leaps) regsubsets.out<-regsubsets(eggsP95NEW~maleorder+matingtime+archnpump+sizeg+remSP,data=MatingNoNA, nbest = 1, # 1 best model for each number of predictors nvmax = NULL, # NULL for no limit on number of variables force.in = NULL, force.out = NULL, method = "exhaustive") regsubsets.out summary.out <- summary(regsubsets.out) as.data.frame(summary.out$outmat) ## Adjusted R2 plot(regsubsets.out, scale = "adjr2", main = "Adjusted R^2") # Find best formula library(car) layout(matrix(1:2, ncol = 2)) ## Adjusted R2 res.legend <- subsets(regsubsets.out, statistic="adjr2", legend = FALSE, min.size = 2, main = "Adjusted R^2") ## Mallow Cp res.legend <- subsets(regsubsets.out, statistic="cp", legend = FALSE, min.size = 2, main = "Mallow Cp") abline(a = 1, b = 1, lty = 2) res.legend which.max(summary.out$adjr2) best.model <- lm(eggsP95NEW~maleorder+matingtime+archnpump+sizeg+remSP,data=MatingNoNA) summary(best.model) #Without REMSP library(leaps) regsubsets.out<-regsubsets(eggsP95NEW~maleorder+matingtime+archnpump+sizeg,data=MatingNoNA, nbest = 1, # 1 best model for each number of predictors nvmax = NULL, # NULL for no limit on number of variables force.in = NULL, force.out = NULL, method = "exhaustive") regsubsets.out summary.out <- summary(regsubsets.out) as.data.frame(summary.out$outmat) ## Adjusted R2 plot(regsubsets.out, scale = "adjr2", main = "Adjusted R^2") # Find best formula library(car) layout(matrix(1:2, ncol = 2)) ## Adjusted R2 res.legend <- subsets(regsubsets.out, statistic="adjr2", legend = FALSE, min.size = 2, main = "Adjusted R^2") ## Mallow Cp res.legend <- subsets(regsubsets.out, statistic="cp", legend = FALSE, min.size = 2, main = "Mallow Cp") abline(a = 1, b = 1, lty = 2) res.legend which.max(summary.out$adjr2) best.model <- lm(eggsP95NEW~maleorder+matingtime+archnpump+sizeg,data=MatingNoNA) summary(best.model) #Try glmulti glmulti.lm.out <- glmulti(eggsP95NEW~maleorder+matingtime+archnpump+sizeg,data=MatingNoNA, level = 1, # No interaction considered method = "h", # Exhaustive approach crit = "aic", # AIC as criteria confsetsize = 5, # Keep 5 best models plotty = F, report = F, # No plot or interim reports fitfunction = "lm") # lm function ## Show 5 best models (Use @ instead of $ for an S4 object) glmulti.lm.out@formulas ## Show result for the best model summary(glmulti.lm.out@objects[[1]]) #MuMIn model.MO = glm.nb(eggs95NEW~maleorder, data=MatingNoNA) model.MT= glm.nb(eggs95NEW~matingtime, data=MatingNoNA) model.AP=glm.nb(eggs95NEW~archnpump, data=MatingNoNA) model.SG= glm.nb(eggs95NEW~sizeg, data=MatingNoNA) library(MuMIn) all.models=list(model.MO,model.MT,model.AP,model.SG) aic.table=model.sel(all.models) aic.table #Include all variables model.MO = glm.nb(eggs95NEW~maleorder, data=MatingNoNA) model.MT= glm.nb(eggs95NEW~matingtime, data=MatingNoNA) model.AP=glm.nb(eggs95NEW~archnpump, data=MatingNoNA) model.SG= glm.nb(eggs95NEW~sizeg, data=MatingNoNA) model.MOMTAP = glm.nb(eggs95NEW~maleorder+matingtime+archnpump, data=MatingNoNA) model.MOMTSG = glm.nb(eggs95NEW~maleorder+matingtime+sizeg, data=MatingNoNA) model.MOAPSG = glm.nb(eggs95NEW~maleorder+archnpump+sizeg, data=MatingNoNA) model.MTAPSG = glm.nb(eggs95NEW~sizeg+matingtime+archnpump, data=MatingNoNA) model.MOMT = glm.nb(eggs95NEW~maleorder+matingtime, data=MatingNoNA) model.MOAP = glm.nb(eggs95NEW~maleorder+archnpump, data=MatingNoNA) model.MOSG = glm.nb(eggs95NEW~maleorder+sizeg, data=MatingNoNA) model.MTAP= glm.nb(eggs95NEW~matingtime+archnpump, data=MatingNoNA) model.MTSG= glm.nb(eggs95NEW~matingtime+sizeg, data=MatingNoNA) model.APSG=glm.nb(eggs95NEW~archnpump+sizeg, data=MatingNoNA) model.MOMTAPSG= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+maleorder, data=MatingNoNA) all.models=list(model.MO,model.MT,model.AP,model.SG,model.MOMTAP,model.MOMTSG,model.MOAPSG,model.MTAPSG,model.MOMT,model.MOAP,model.MOSG,model.MTAP,model.APSG,model.MOMTAPSG) aic.table=model.sel(all.models) aic.table importance(all.models) #Graph most important model model.MOSG = glm.nb(eggs95NEW~maleorder+sizeg, data=MatingNoNA) plot(model.MOSG) #With Remsp model.MO = glm.nb(eggs95NEW~maleorder, data=MatingNoNA) model.MT= glm.nb(eggs95NEW~matingtime, data=MatingNoNA) model.AP=glm.nb(eggs95NEW~archnpump, data=MatingNoNA) model.SG= glm.nb(eggs95NEW~sizeg, data=MatingNoNA) model.RS=glm.nb(eggs95NEW~remSP, data=MatingNoNA) model.MOMTAP = glm.nb(eggs95NEW~maleorder+matingtime+archnpump, data=MatingNoNA) model.MOMTSG = glm.nb(eggs95NEW~maleorder+matingtime+sizeg, data=MatingNoNA) model.MOAPSG = glm.nb(eggs95NEW~maleorder+archnpump+sizeg, data=MatingNoNA) model.MOMTRS = glm.nb(eggs95NEW~maleorder+matingtime+remSP, data=MatingNoNA) model.MOAPRS = glm.nb(eggs95NEW~maleorder+archnpump+remSP, data=MatingNoNA) model.MOSGRS = glm.nb(eggs95NEW~maleorder+sizeg+remSP, data=MatingNoNA) model.MTAPSG = glm.nb(eggs95NEW~sizeg+matingtime+archnpump, data=MatingNoNA) model.MTAPRS = glm.nb(eggs95NEW~remSP+matingtime+archnpump, data=MatingNoNA) model.MTSGRS = glm.nb(eggs95NEW~sizeg+matingtime+remSP, data=MatingNoNA) model.APSGRS = glm.nb(eggs95NEW~sizeg+remSP+archnpump, data=MatingNoNA) model.MOMT = glm.nb(eggs95NEW~maleorder+matingtime, data=MatingNoNA) model.MOAP = glm.nb(eggs95NEW~maleorder+archnpump, data=MatingNoNA) model.MOSG = glm.nb(eggs95NEW~maleorder+sizeg, data=MatingNoNA) model.MORS = glm.nb(eggs95NEW~maleorder+remSP, data=MatingNoNA) model.MTAP= glm.nb(eggs95NEW~matingtime+archnpump, data=MatingNoNA) model.MTSG= glm.nb(eggs95NEW~matingtime+sizeg, data=MatingNoNA) model.MTRS= glm.nb(eggs95NEW~matingtime+remSP, data=MatingNoNA) model.APSG=glm.nb(eggs95NEW~archnpump+sizeg, data=MatingNoNA) model.APRS=glm.nb(eggs95NEW~archnpump+remSP, data=MatingNoNA) model.SGRS= glm.nb(eggs95NEW~sizeg+remSP, data=MatingNoNA) model.MOMTAPSG= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+maleorder, data=MatingNoNA) model.MTAPRSSG= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+remSP, data=MatingNoNA) model.MOMTAPRS= glm.nb(eggs95NEW~remSP+archnpump+matingtime+maleorder, data=MatingNoNA) model.MOMTRSSG= glm.nb(eggs95NEW~sizeg+remSP+matingtime+maleorder, data=MatingNoNA) model.MORSAPSG= glm.nb(eggs95NEW~sizeg+archnpump+remSP+maleorder, data=MatingNoNA) model.MOMTAPSGRS= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+maleorder+remSP, data=MatingNoNA) all.models=list(model.MO,model.MT,model.AP,model.SG,model.RS,model.MOMTAP,model.MOMTSG,model.MOAPSG,model.MOMTRS,model.MOAPRS,model.MOSGRS,model.MTAPSG,model.MTAPRS,model.MTSGRS,model.APSGRS,model.MORS,model.MOMT,model.MOAP,model.MOSG,model.MTAP,model.MTRS,model.APRS,model.SGRS,model.APSG,model.MOMTAPSG,model.MOMTRSSG,model.MORSAPSG,model.MOMTAPRS,model.MTAPRSSG,model.MOMTAPSGRS) aic.table=model.sel(all.models) aic.table importance(all.models) #ANOVA big.model = glm(eggs95NEW ~ maleorder + sizeg + archnpump + matingtime, data = MatingNoNA) plot(allEffects(big.model)) Anova(big.model) #ANOVA with Stepwise step.model = step(big.model) Anova(step.model) #Best AIC Model library(MuMIn) big.model = glm(eggs95NEW ~ maleorder + sizeg + archnpump + matingtime, data = MatingNoNA, na.action=na.pass) sole.dredge = dredge(big.model, extra = "R^2") #Plot with RemSP big.model = glm(eggs95NEW ~ maleorder + sizeg + archnpump + matingtime + remSP, data = MatingNoNA) plot(allEffects(big.model)) Anova(big.model) plot(allEffects(big.model), ylab= "Percentage", ylim=c(-0.2, 0.8)) #Plot with m2msize big.model = glm(eggsP95NEW ~ maleorder + sizeg + archnpump + matingtime + remSP + m2fsize, data = MatingNoNA) plot(allEffects(big.model), ylim=c(-0.4, 0.8), ylab="Percentage") Anova(big.model) model.MO = glm.nb(eggs95NEW~maleorder, data=MatingNoNA) model.MT= glm.nb(eggs95NEW~matingtime, data=MatingNoNA) model.AP=glm.nb(eggs95NEW~archnpump, data=MatingNoNA) model.SG= glm.nb(eggs95NEW~sizeg, data=MatingNoNA) model.RS=glm.nb(eggs95NEW~remSP, data=MatingNoNA) model.M2F=glm.nb(eggs95NEW~m2f, data=MatingNoNA) model.MOMTAP = glm.nb(eggs95NEW~maleorder+matingtime+archnpump, data=MatingNoNA) model.MOMTSG = glm.nb(eggs95NEW~maleorder+matingtime+sizeg, data=MatingNoNA) model.MOAPSG = glm.nb(eggs95NEW~maleorder+archnpump+sizeg, data=MatingNoNA) model.MOMTRS = glm.nb(eggs95NEW~maleorder+matingtime+remSP, data=MatingNoNA) model.MOAPRS = glm.nb(eggs95NEW~maleorder+archnpump+remSP, data=MatingNoNA) model.MTAPSG = glm.nb(eggs95NEW~sizeg+matingtime+archnpump, data=MatingNoNA) model.MTSGRS = glm.nb(eggs95NEW~sizeg+matingtime+remSP, data=MatingNoNA) model.MTAPRS = glm.nb(eggs95NEW~sizeg+remSP+archnpump, data=MatingNoNA) model.MOMT = glm.nb(eggs95NEW~maleorder+matingtime, data=MatingNoNA) model.MOAP = glm.nb(eggs95NEW~maleorder+archnpump, data=MatingNoNA) model.MOSG = glm.nb(eggs95NEW~maleorder+sizeg, data=MatingNoNA) model.MORS = glm.nb(eggs95NEW~maleorder+remSP, data=MatingNoNA) model.MOM2F=glm.nb(eggs95NEW~maleorder+M2F, data=MatingNoNA) model.MTAP= glm.nb(eggs95NEW~matingtime+archnpump, data=MatingNoNA) model.MTSG= glm.nb(eggs95NEW~matingtime+sizeg, data=MatingNoNA) model.MTRS= glm.nb(eggs95NEW~matingtime+remSP, data=MatingNoNA) model.MTM2F=glm.nb(eggs95NEW~matingtime+M2F, data=MatingNoNA) model.APSG=glm.nb(eggs95NEW~archnpump+sizeg, data=MatingNoNA) model.APRS=glm.nb(eggs95NEW~archnpump+remSP, data=MatingNoNA) model.APM2F=glm.nb(eggs95NEW~archnpump+M2F, data=MatingNoNA) model.SGRS= glm.nb(eggs95NEW~sizeg+remSP, data=MatingNoNA) model.SGM2F =glm.nb(eggs95NEW~sizeg+M2F, data=MatingNoNA) model.RSM2F=glm.nb(eggs95NEW~M2F+remSP, data=MatingNoNA) model.MOMTAPSG= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+maleorder, data=MatingNoNA) model.MTAPRSSG= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+remSP, data=MatingNoNA) model.MOMTAPRS= glm.nb(eggs95NEW~remSP+archnpump+matingtime+maleorder, data=MatingNoNA) model.MOMTRSSG= glm.nb(eggs95NEW~sizeg+remSP+matingtime+maleorder, data=MatingNoNA) model.MORSAPSG= glm.nb(eggs95NEW~sizeg+archnpump+remSP+maleorder, data=MatingNoNA) model.MOMTAPSGRS= glm.nb(eggs95NEW~sizeg+archnpump+matingtime+maleorder+remSP, data=MatingNoNA)