wheat <- read.table("C:/temp/wheatear tag.csv", header=T, sep="," ) head(wheat) dim(wheat) wheat$sex <- factor(wheat$sex) wheat$age <- factor(wheat$age) modela <- glm(surv ~ age + sex + order + wing + permass + stalk, family=binomial, data=wheat) summary(modela) #######Run dredge library(MuMIn) options(na.action = "na.fail") dd <- dredge(modela, extra = c("R^2")) top.models <- get.models(dd, subset = delta <4) head(dd, 10) avmodel <- model.avg(top.models) summary(avmodel) table(wheat$sex, wheat$age) ####Logger bird versus control data=rbind(c(14,10),c(45,36)) data chisq.test(data) data1=rbind(c(24,14),c(81,45)) data1 data1 colnames(data1)=c("Available","Resighted") rownames(data1)=c("Logger","Control") barplot(t(data1), xlab="Type", main="Geolocator tag effect 2015", ylab="Frequency", beside=TRUE,las=1) legend(1,60, colnames(data1), lty = c(-1, -1, -1), fill=c("black","lightgrey"),bty = 1 , cex=0.75) text(2.5,40, "Chi-square = 0.9, P = 0.34") text(5.5,40, "55.5%") text(2.5,20, "58.3%") propcont1 <- 45/81 propcont1 proplog1 <- 14/24 proplog1 data6=rbind(c(propcont1, proplog1)) data6 colnames(data6)=c("Control","Logger") barplot(t(data6), xlab="", ylab="Proportion resighted the following year", beside=TRUE, ylim=c(0,0.8)) legend(1,0.8, colnames(data6), lty = c(-1, -1, -1), fill=c("black","lightgrey"),bty = 1 , cex=1,) text(2,0.7, "N = 45/81 control vs 12/24 logger") text(2,0.6, "Chi-square < 0.01, P = 0.99") wheat1 <- read.table("C:/temp/Wheatear Tag effect ringing data 2015.csv", header=T, sep="," ) head(wheat1) dim(wheat1) str(wheat1) options(na.action = "na.omit") model1 <- lm(mass ~ geo, data=wheat1) summary(model1) model2 <- lm(wing ~ geo, data=wheat1) summary(model2) table(wheat1$age) wheat1$age1 <- factor(wheat1$age) wheat1$sex1 <- factor(wheat1$sex) dim(wheat1) wheat2 <- subset(wheat1, mass >-1) wheat3 <- subset(wheat2, age >-1) wheat4 <- subset(wheat3, wing >-1) wheat5 <- subset(wheat4, sex >-1) dim(wheat5) model3 <- glm(survive ~ geo + age1 + sex1 + wing + mass + sex1*mass, family=binomial, data=wheat5) summary(model3) model3b <- glm(survive ~ geo + age1 + sex1 + wing + mass + geo*mass, family=binomial, data=wheat5) summary(model3b) model3a <- glm(survive ~ geo + age1 + sex1 + wing + mass, family=binomial, data=wheat5) summary(model3a) library(MuMIn) options(na.action = "na.fail") dd <- dredge(model3, extra = c("R^2")) top.models <- get.models(dd, subset = delta <4) head(dd, 10) avmodel <- model.avg(top.models) summary(avmodel) ##second model model4 <- glm(survive ~ age1 + mass + wing , family=binomial, data=wheat5) summary(model4) ##################################################################################################### ####Logger bird versus control 2014 data=rbind(c(279,73),c(130,38)) data chisq.test(data) data # print table, remember columns are head up, head down, and rows are male and female colnames(data)=c("Available","Resighted") rownames(data)=c("Control","Logger") barplot(t(data), xlab="Type", main="Geolocator tag effect 2014", ylab="Frequency", beside=TRUE) legend(5,250, colnames(data), lty = c(-1, -1, -1), fill=c("black","lightgrey"),bty = 1 , cex=0.75,) text(4,200, "Chi-square = 0.14, P = 0.70") text(2.5,80, "26.1%") text(5.5,50, "29.2%") propcont2 <- 73/279 propcont2 proplog2 <- 38/130 proplog2 ################Annual return rate proportions data4=rbind(c(propcont1,proplog1),c(propcont2,proplog2)) data4 colnames(data4)=c("Control","Logger") rownames(data4)=c("2013","2014") barplot(t(data4), xlab="", main="Geolocator tag effect 2013 & 2014", ylab="Proportion resighted the following winter", beside=TRUE, ylim=c(0,0.50)) legend(5,0.5, colnames(data4), lty = c(-1, -1, -1), fill=c("black","lightgrey"),bty = 1 , cex=0.75,) text(4,0.4, "Chi-square = 0.14, P = 0.70") ###Overall difference by year data5=rbind(c(73,279),c(38,130)) data5 chisq.test(data5) #####Overall pooled proportion return by year totpropc <- (73+16)/(279+37) totpropl <- (38+10)/(130 + 26) totpropc totpropl data6=rbind(c(totpropc, totpropl)) data6 chisq.test(data6) colnames(data6)=c("Control","Logger") barplot(t(data6), xlab="", main="Geolocator tag effect 2013 & 2014 pooled", ylab="Proportion resighted the following winter", beside=TRUE, ylim=c(0,0.50)) legend(1,0.5, colnames(data6), lty = c(-1, -1, -1), fill=c("black","lightgrey"),bty = 1 , cex=1,) text(2,0.38, "N = 89/316 control vs 48/156 logger") text(2,0.41, "Chi-square = 0.0012, P = 0.97") ########################Harness effect ####code 1 long elas; 2 short elas; 3 no elas; 4 long tie; 5 short tie; 6 no tie #### material 1 elastic; 2 tie #### stalk 3 long; 2 short; 1 none data8=rbind(c(10,20),c(1,16),c(13,18),c(4,12),c(4,9),c(6,17)) data8 chisq.test(data8) colnames(data8)=c("Resighted","Available") rownames(data8)=c("Longelas","Longtie","Shortelas","Shorttie","Noneelas","Nonetie") barplot(t(data8), xlab="Logger light stalk length and attachment material", main="Geolocator harness effect 2014", ylab="Frequency", beside=TRUE) legend(12,20, colnames(data8), lty = c(-1, -1, -1), fill=c("black","lightgrey"),bty = 1 , cex=0.75,) text(4,200, "Chi-square = 0.42, P = 0.51") fisher.test(data8) dat <- read.table("C:/temp/harness.csv", header=T, sep="," ) head(dat) dim(dat) dat$year <- factor(dat$year) dat$material <- factor(dat$material) dat$stalk1 <- factor(dat$stalk) ###Any difference by year? dat1 <- subset(dat, material==1 & stalk==3) dim(dat1) dat1 modela <- glm(return ~ year, family=binomial, data=dat1) summary(modela) data2=rbind(c(10,20),c(10,16)) data2 chisq.test(data2) type1 <- -0.4700 type1u <- -0.4700 +0.4031 type1l <- -0.4700 - 0.4031 type11 <- (exp(type1))/(1 + (exp(type1))) type11u <- (exp(type1u))/(1 + (exp(type1u))) type11l <- (exp(type1l))/(1 + (exp(type1l))) type2 <- -0.4700 + -0.2231 type2u <- -0.4700 + 0.5590 -0.2231 type2l <- -0.4700 - 0.5590 -0.2231 type22 <- (exp(type2))/(1 + (exp(type2))) type22u <- (exp(type2u))/(1 + (exp(type2u))) type22l <- (exp(type2l))/(1 + (exp(type2l))) ### Plot par(mfrow=c(1,1)) y.meansall <- c(type11, type22) barx <- barplot(y.meansall,ylim=c(0,0.52), names.arg=c("2013", "2014"), ylab="Probability of resighting +/- 1SE", xlab="Year") lower <- c(type11u, type22u) upper <- c(type11l, type22l) ##################install Hmisc library(Hmisc) errbar(barx[,1], y.meansall, lower, upper, add=T, xlab="", lty=1, lwd=1) text(0.8,0.51, "N = 26") text(1.8,0.51, "N = 30") ###############Overall model harness type model1 <- glm(return ~ material + stalk1 + material*stalk1, data=dat, family=binomial) summary(model1) drop1(model1,~.,test="F") dim(dat) model2 <- glm(return ~ material + stalk1, data=dat, family=binomial) summary(model2) table(dat$material) model3 <- glm(return ~ material, data=dat, family=binomial) summary(model3) AIC(model1,model2,model3) ####################Medium stalk predictions for material summary(model2) type1 <- -0.3854 + 0.1127 type1u <- -0.3854 + 0.4459 + 0.1127 type1l <- -0.3854 - 0.4459 + 0.1127 type11 <- (exp(type1))/(1 + (exp(type1))) type11u <- (exp(type1u))/(1 + (exp(type1u))) type11l <- (exp(type1l))/(1 + (exp(type1l))) type2 <- -0.3854 + -0.9640 + 0.1127 type2u <- -0.3854 + -0.9640 + 0.4236+ 0.1127 type2l <- -0.3854 + -0.9640 - 0.4236+ 0.1127 type22 <- (exp(type2))/(1 + (exp(type2))) type22u <- (exp(type2u))/(1 + (exp(type2u))) type22l <- (exp(type2l))/(1 + (exp(type2l))) ### Plot par(mfrow=c(1,1)) y.meansall <- c(type11, type22) barx <- barplot(y.meansall,ylim=c(0,0.6), names.arg=c("Elastic", "Tie"), ylab="Probability of resighting +/- 1SE", xlab="Attachment material") lower <- c(type11u, type22u) upper <- c(type11l, type22l) ##################install Hmisc errbar(barx[,1], y.meansall, lower, upper, add=T, xlab="", lty=1, lwd=1) text(0.8,0.58, "N = 100") text(1.8,0.58, "N = 56") ####################Elastic - different stalk lengths summary(model2) type1 <- -0.3854 type1u <- -0.3854 + 0.4459 type1l <- -0.3854 - 0.4459 type11 <- (exp(type1))/(1 + (exp(type1))) type11u <- (exp(type1u))/(1 + (exp(type1u))) type11l <- (exp(type1l))/(1 + (exp(type1l))) type2 <- -0.3854 + 0.1127 type2u <- -0.3854 + 0.1127 + 0.5034 type2l <- -0.3854 + 0.1127 - 0.5034 type22 <- (exp(type2))/(1 + (exp(type2))) type22u <- (exp(type2u))/(1 + (exp(type2u))) type22l <- (exp(type2l))/(1 + (exp(type2l))) type3 <- -0.3854 + -0.3343 type3u <- -0.3854 + -0.3343 + 0.4905 type3l <- -0.3854 + -0.3343 - 0.4905 type33 <- (exp(type3))/(1 + (exp(type3))) type33u <- (exp(type3u))/(1 + (exp(type3u))) type33l <- (exp(type3l))/(1 + (exp(type3l))) y.meansall <- c(type11, type22, type33) barx <- barplot(y.meansall,ylim=c(0,0.6), names.arg=c("None", "Medium","Long"), ylab="Probability of resighting +/- 1SE", xlab="Logger stalk length") lower <- c(type11u, type22u, type33u) upper <- c(type11l, type22l, type33l) ##################install Hmisc errbar(barx[,1], y.meansall, lower, upper, add=T, xlab="", lty=1, lwd=1) text(0.7,0.58, "N = 36") text(1.9,0.58, "N = 47") text(3.1,0.58, "N = 73") table(dat$stalk) ###Stalk as covariate dat$stalkc[dat$stalk == 1] <- 0 dat$stalkc[dat$stalk == 2] <- 5 dat$stalkc[dat$stalk == 3] <- 10 model11 <- glm(return ~ material + stalkc + material*stalkc, data=dat, family=binomial) summary(model11) drop1(model11,~.,test="F") model22 <- glm(return ~ material + stalkc, data=dat, family=binomial) summary(model22) AIC(model22, model2) ######YOUR MODEL summary(model22) ###########specify your variables to work with generic code below dat$y <- dat$return ######Add your y variable dat$x <- dat$stalkc ######Add your x variable dat$v2 <- dat$material ########Add your parameter variables sumtab <- coef(summary(model22)) write.csv(sumtab,"c:/temp/summarytable.csv") test10 <- read.csv("c:/temp/summarytable.csv") test10 library(reshape) ## rename variables test1 <- rename(test10, c(X="X", Estimate="est", Std..Error="se" )) test1 int1 <- test1[1,c("est")] int1se<- test1[1,c("se")] par2 <- test1[2,c("est")] par2se<- test1[2,c("se")] par3 <- test1[3,c("est")] par3se<- test1[3,c("se")] ###Predicted values ##X values for predicted lines min<-min(dat$x) min max<-max(dat$x) max newx<-seq(max,min,length=1000) ########################################################################################## #######################################PLOTTING RELATIONSHIP WITH DISTANCE ####Predicted values for stalk length from the model above (when material=TIE) pmean <- int1 + par2 + (newx* par3) plower <- int1 +int1se + par2 + (newx* (par3+par3se)) pupper <- int1-int1se + par2 + (newx* (par3-par3se)) ####Predicted values for stalk length from the model above (when material=ELASTIC) pmeana <- int1 + (newx* par3) plowera <- int1 +int1se + (newx* (par3+par3se)) puppera <- int1 -int1se + (newx* (par3-par3se)) YAxisLabel <- "Probability of resighting in following year +/- 1SE" XAxisLabel <- "Stalk length (mm) "########## #####Back transform because of log link binomial function ##TIE pmean1 <- (exp(pmean))/(1 + (exp(pmean))) plower1 <- (exp(plower))/(1 + (exp(plower))) pupper1 <- (exp(pupper))/(1 + (exp(pupper))) ###ELASTIC pmean1a <- (exp(pmeana))/(1 + (exp(pmeana))) plower1a <- (exp(plowera))/(1 + (exp(plowera))) pupper1a <- (exp(puppera))/(1 + (exp(puppera))) par(mfrow=c(1,1)) sunflowerplot(dat$x, dat$y, xlim=c(0,10), ylab=YAxisLabel, xlab=XAxisLabel, frame.plot=FALSE, add=F, seg.col="grey") text(1,0.4, "Elastic") text(1,0.2, "Tie") lines(newx, pmean1, col="red", lwd=2) lines(newx, plower1,lty=3, col="red") lines(newx, pupper1, lty=3, col="red") lines(newx, pmean1a, col="blue", lwd=2) lines(newx, plower1a,lty=3, col="blue") lines(newx, pupper1a, lty=3, col="blue")