binplot<- function (size,surv, ncuts = 70,...) { os <- order(size) os.surv <- (surv)[os] os.size <- (size)[os] psz <- tapply(os.size, as.numeric(cut(os.size, ncuts)), mean, na.rm = TRUE) ps <- tapply(os.surv, as.numeric(cut(os.size, ncuts)), mean, na.rm = TRUE) plot(as.numeric(psz), as.numeric(ps),...) } #convert between population and individual data expand.dft <- function(x, na.strings = "NA", as.is = FALSE, dec = ".") { # Take each row in the source data frame table and replicate it # using the Freq value DF <- sapply(1:nrow(x), function(i) x[rep(i, each = x\$Freq[i]), ], simplify = FALSE) # Take the above list and rbind it to create a single DF # Also subset the result to eliminate the Freq column DF <- subset(do.call("rbind", DF), select = -Freq) # Now apply type.convert to the character coerced factor columns # to facilitate data type selection for each column for (i in 1:ncol(DF)) { DF[[i]] <- type.convert(as.character(DF[[i]]), na.strings = na.strings, as.is = as.is, dec = dec) } DF } #NOTE THAT FOR GLMS YOU CAN USE PREDICT INSTEAD OF X lz<-function(x) {log(x+1)} RSS1<-function(m,x) {sum((lz((m))-lz(x))^2)} TSS1<-function(x) {sum((lz(x)-lz(mean(x)))^2)} RSS<-function(m,x) {sum((((m))-(x))^2)} TSS<-function(x) {sum(((x)-(mean(x)))^2)} R2<-function(m,x,logo=F) {if(logo==F) {return(1-RSS(m,x)/TSS(x))} else {return(1-RSS1(m,x)/TSS1(x))} } #summary[[2]] gisum<-function(model) {return(model[[2]]\$summary[1:20,])} gipar<-function(model) {return(model[2]\$BUGSoutput\$sims.list)} Rhat<-function(model) { return(model[[2]]\$summary[,8]) } densplot<-function(model,variable,colo="light gray") { simlist<-gipar(model) simvar<-which(names(simlist)==variable) plot(density(simlist[[simvar]]),main=variable) polygon(density(simlist[[simvar]]),col=colo) } pplotCOEF<-function(sims,lab1="Parameters",cexo=0.5,lab2="Values",listorno="T") { #this is designed to plot the 95% confidence intervals if(listorno=="T") {simlist<-sims simlist\$deviance<-NULL } else{ simlist<-gipar(sims) dumb<-which(names(simlist)=="deviance") simlist<-simlist[-dumb] } variablenames<-names(simlist) xo<-c(1:length(names(simlist))) upo<-max(unlist(lapply(simlist,ups))) downo<-min(unlist(lapply(simlist,lows))) plot(0,xo[1],pch="",ylim=c(downo,upo),xlim=c(1,max(xo)),xlab=lab1,ylab=lab2,xaxt="n") axis(side=1,at=c(1:max(xo)),labels=variablenames) abline(h=0,lty=2,col="gray") for(i in 1:length(xo)) { valz<-simlist[[i]] m<-median(valz) u<-ups(valz) l<-lows(valz) points(m~xo[i],pch=19,cex=cexo,col="black") segments(x0=xo[i],y0=l,y1=u,lty=1,col="black") } } pplot2<-function(input,pvalz,lab1="Predictor",cexo=0.5,lab2="Response",ylimo=c(min(l),max(u)), xlimo=c(min(input),max(input)),maino="",meormed="med") { #this is designed to plot the 95% confidence intervals if(meormed=="me") { m<-apply(pvalz,2,mean) } else{m<-apply(pvalz,2,median)} u<-apply(pvalz,2,ups) l<-apply(pvalz,2,lows) have<-data.frame(l,m,u,input) colnames(have)<-c("lower","mean","upper","input") have<-have[order(have\$input),] m<-have\$m l<-have\$l u<-have\$u ip<-have\$input plot(m~ip,ylim=ylimo,xlim=xlimo,pch=19,cex=cexo,xlab=lab1,ylab=lab2,col="black",main=maino,cex.lab=1.8,cex.axis=1.8) #lines(m~ip,lty=2) for(i in 1:length(m)) { segments(x0=ip[i],y0=l[i],y1=u[i],lty=2,col="gray") #points(m~ip,pch=19,cex=cexo,col="black") } } ups<-function(listie) { want<-length(listie)*0.975 take<-sort(listie)[want] return(take) } lows<-function(listie) { want<-length(listie)*0.025 take<-sort(listie)[want] return(take) } ppointFX<-function(input,predicted) { #this is designed for plotting the 95% confidence intervals around a functional relationship pvalz<-predicted m<-apply(pvalz,2,median) u<-apply(pvalz,2,ups) l<-apply(pvalz,2,lows) have<-data.frame(l,m,u,input) colnames(have)<-c("lower","mean","upper","input") have<-have[order(have\$input),] m<-have\$m l<-have\$l u<-have\$u ip<-have\$input points(m~ip,ylim=c(min(l),max(u)),pch=19,col="gray") lines(m~ip,lty=2) #points(u~ip,pch=19,col="gray") points is to confusing lines(u~ip,lty=2,col="red") #the lines might be good if #points(l~ip,pch=19,col="gray") lines(l~ip,lty=2,col="red") return(have) } pplotFX<-function(input,predicted,lab1,lab2) { #this is designed for plotting the 95% confidence intervals around a functional relationship pvalz<-predicted m<-apply(pvalz,2,median) u<-apply(pvalz,2,ups) l<-apply(pvalz,2,lows) have<-data.frame(l,m,u,input) colnames(have)<-c("lower","mean","upper","input") have<-have[order(have\$input),] m<-have\$m l<-have\$l u<-have\$u ip<-have\$input plot(m~ip,ylim=c(min(l),max(u)),pch=19,col="gray",xlab=lab1,ylab=lab2) lines(m~ip,lty=2) #points(u~ip,pch=19,col="gray") points is to confusing lines(u~ip,lty=2,col="red") #the lines might be good if #points(l~ip,pch=19,col="gray") lines(l~ip,lty=2,col="red") return(have) } pplot<-function(input,predicted,lab1,lab2) { #this is designed to plot the 95% confidence intervals pvalz<-predicted m<-apply(pvalz,2,median) u<-apply(pvalz,2,ups) l<-apply(pvalz,2,lows) have<-data.frame(l,m,u,input) colnames(have)<-c("lower","mean","upper","input") have<-have[order(have\$input),] m<-have\$m l<-have\$l u<-have\$u ip<-have\$input plot(m~ip,ylim=c(min(l),max(u)),pch=19,col="gray",xlab=lab1,ylab=lab2) #lines(m~ip,lty=2) for(i in 1:length(m)) { segments(x0=ip[i],y0=l[i],y1=u[i],lty=2,col="gray70") } #points(u~ip,pch=19,col="gray") points is to confusing #lines(u~ip,lty=2,col="red") #the lines might be good if #points(l~ip,pch=19,col="gray") #lines(l~ip,lty=2,col="red") return(have) } ppoint<-function(input,predicted,lab1,lab2) { #this is designed to plot the 95% confidence intervals pvalz<-predicted m<-apply(pvalz,2,median) u<-apply(pvalz,2,ups) l<-apply(pvalz,2,lows) have<-data.frame(l,m,u,input) colnames(have)<-c("lower","mean","upper","input") have<-have[order(have\$input),] m<-have\$m l<-have\$l u<-have\$u ip<-have\$input points(m~ip,ylim=c(min(l),max(u)),pch=19,col="gray") #lines(m~ip,lty=2) for(i in 1:length(m)) { segments(x0=ip[i],y0=l[i],y1=u[i],lty=2,col="gray70") } #points(u~ip,pch=19,col="gray") points is to confusing #lines(u~ip,lty=2,col="red") #the lines might be good if #points(l~ip,pch=19,col="gray") #lines(l~ip,lty=2,col="red") return(have) } #I'm pretty sure this is wrong-o, wrong-o! binMSPE<-function(y,sims.p) { p<-apply(sims.p,2,mean) mspe<-sum((p-y)^2)+sum(p*(1-p)) m<-sum((p-y)^2) great<-data.frame(mspe,m) names(great)<-c("MSPE","mean only") return(great) } #(2/(length(y)*(length(y)+1))) #an alternative. I'm not sure which of these to use poisMSPE<-function(y,sims.p,its=length(sims.p[1,])) { runny<-rep(0,times=its) for(i in 1:its) { runny[i]<-(1/(length(y)))*sum((log(sims.p[,i]+0.5)-log(y+0.5))^2) } return(runny) } poisMSPE1it<-function(y,sims.p) { runny<-(1/(length(y)))*sum((log(sims.p+0.5)-log(y+0.5))^2) return(runny) } normalMSPE<-function(y,sims.p,sigma) { p<-apply(sims.p,2,mean) mspe<-sum((p-y)^2)+sum(sigma) m<-sum((p)^2) great<-data.frame(mspe,m) names(great)<-c("MSPE","mean only") return(great) } #OR IF SIGMA IS VARIABLE: normalMSPEvar<-function(y,sims.p,sigma) { p<-apply(sims.p,2,mean) sigma.m<-apply(sigma,2,mean) mspe<-sum((p-y)^2)+sum(sigma.m) m<-sum((p)^2) great<-data.frame(mspe,m) names(great)<-c("MSPE","mean only") return(great) } #what's next: #1. experiment with value for total seeds #2. put in covariates for final model! yay! #3.. dealing with missing values ###NOTES: #the lmer model for site with a binomial model is 0 for all random effects...seems problematic credint<-function(x) { gret<-sort(x) up<-round(length(x)*0.95) down<-round(length(x)*0.05) final<-c(gret[down],median(gret),mean(gret),gret[up]) names(final)<-c("lower 5%","median","mean","upper 95%") return(final) } stdize<-function(x) {(x-mean(x,na.rm=T))/(2*sd(x,na.rm=T))} transdis<-function(x,y) {sqrt((521995-x)^2+(1727781-y)^2) }