search() ls(2) dim(2) dim(agb) head(agb) plot(agb~ht,data=agb) plot(I(log(agb))~I(log(ht)),data=agb) lm(agb~ForestType+ht,data=agb) ls(2) ditch(2) attach('~/Desktop/agb.rdata') ls(2) head(agb) lm(agb~ForestType+ht,data=agb) lm(agb~ForestType+ht,data=agb) agb$logdbh=log(agb$dbh) agb$logagb=log(agb$agb) TreeMass=agb save(TreeMass,file='~/Desktop/TreeMass.rdata') head(TreeMass) table(TreeMass$ForestType) tapply(TreeMass$Rainfall,TreeMass$ForestType,range) tapply(TreeMass$Alt,TreeMass$ForestType,range) weird=TreeMass$ForestType=='Moist-' table(weird) TreeMass$ForestType[weird] weird=TreeMass$ForestType=='Moist-' & !is.na(TreeMass$ForestType ) weird=TreeMass$ForestType=='Moist-' & !is.na(TreeMass$ForestType) TreeMass$ForestType[weird] TreeMass$ForestType[weird]='Moist' save(TreeMass,'~/meetings_workshops/Rmodeling/treemass.rdata') save(TreeMass,'~/meetings_workshops/Rmodeling/treemass.rdata') treemass=TreeMass save(treemass,file='~/meetings_workshops/Rmodeling/treemass.rdata') head(treemass) lm(logagb~logdbh+ForestType,data=treemass) fit=lm(logagb~logdbh+ForestType,data=treemass) coef(fit) abline(coef(fit)[1:2]) abline(coef(fit)[1],sum(coef(fit)[2:3])) abline(sum(coef(fit)[c(1,3)]),coef(fit)[2]) ls(2) plot(logagb~logdbh,data=treemass) abline(sum(coef(fit)[c(1,3)]),coef(fit)[2]) abline(sum(coef(fit)[c(1,3)]),coef(fit)[2],col='red') sum(coef(fit)[c(1,3)]) abline(coef(fit)[1:2],col='red') abline(sum(coef(fit)[c(1,3)]),coef(fit)[2],col='green') abline(sum(coef(fit)[c(1,4)]),coef(fit)[2],col='blue') plot(logagb~logdbh,data=treemass,cex=.5,pch=16) abline(sum(coef(fit)[c(1,4)]),coef(fit)[2],col='blue') abline(sum(coef(fit)[c(1,3)]),coef(fit)[2],col='green') abline(coef(fit)[1:2],col='red') treemass$ctrx=treemass$logdbh-mean(treemass$logdbh) plot(logagb~ctrx,data=treemass,cex=.5,pch=16) fit=lm(logagb~ctrx+ForestType,data=treemass) abline(coef(fit)[1:2],col='red') coef(fit) system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops") system("~/meetings_workshops/Rmodeling/CreateCourse.bat outline Tupper2012 Public/Workshops") treemass$ctrx=treemass$logdbh-mean(treemass$logdbh) plot(logagb~logdbh,data=treemass,cex=.5,pch=16) fit=lm(logagb~logdbh,data=treemass) coef(fit) abline(coef(fit)) abline(coef(fit),col='red') fit=lm(logagb~ctrx+ForestType,data=treemass) abline(fit) fit=lm(logagb~logdbh,data=treemass) plot(logagb~logdbh,data=treemass,cex=.5,pch=16) fit=lm(logagb~logdbh,data=treemass) search() options(width=125) search() options(width=80) search() ls(6) ls() ls(6) search() ls(6) head(treeht) search() ls(2) head(treeheight) head(treemass) search()[1:4] attach("~/meetings_workshops/Rmodeling/treemass.rdata") rm(treemass) ditch(2) attach("~/meetings_workshops/Rmodeling/treemass.rdata") ls(2) head(treemass) options(width=60) head(treemass) options(width=65) head(treemass) ls(2) detach(2) attach("~/meetings_workshops/Rmodeling/treemass.rdata") ls(2) head(treemass) fit=lm(logagb~logdbh,data=treemass) plot(logagb~logdbh,data=treemass,cex=.5,pch=16) plot(treemass$logagb~treemass$logdbh,cex=.5,pch=16) summary(fit) coef(fit) fit=lm(logagb~logdbh,data=treemass) coef(fit) abline(fit) abline(fit,col='red') plot(logagb~logdbh,data=treemass,cex=.5,pch=16) fit=lm(logagb~logdbh,data=treemass) abline(fit) abline(fit,col='red') coef(fit) source('~/meetings_workshops/Rmodeling/teaching.functions.r') args(sumsq.linearmodel) sumsq.linearmodel(beta=coef(fit),x=treemass$logdbh,y=treemass$logagb) sumsq.linearmodel(beta=coef(fit),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(5,0),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(5,2),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(5,2),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(-5,2),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(0,.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(0,.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) sumsq.linearmodel(beta=c(0,1.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) ?optim ?optim optim(par=c(0,0),fn=sumsq.linearmodel) ?optim optim(par=c(0,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb)) optim(par=c(0,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb) sumsq.linearmodel(beta=c(0,1.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE) ?optim optim(par=c(0,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb) coef(fit) optim(par=c(10,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb) optim(par=c(10,-10),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb) ?optim ?optim control=list(fnscale=(-1))q ?optim llike.linearmodel(param=c(-2,2),x=treeheight$logdbh,y=treeheight$logdbh) dnorm(1,mean=0,sd=1) dnorm(0:5,mean=0,sd=1) dnorm(-5:5,mean=0,sd=1) plot(-5:5,dnorm(-5:5,mean=0,sd=2) ) plot(-5:5,dnorm(-5:5,mean=0,sd=2,log=TRUE)) ?dnorm plot(-5:5,dnorm(-5:5,mean=1:11,sd=2,log=TRUE)) dnorm(-5:5,mean=1:11,sd=2,log=TRUE) dnorm(-5:5,mean=1,sd=2,log=TRUE) dnorm(-5:5,mean=1:11,sd=1:11,log=TRUE) llike.linearmodel(param=c(-2,2),x=treeheight$logdbh,y=treeheight$logdbh) llike.linearmodel(param=c(-2,2,1),x=treeheight$logdbh,y=treeheight$logdbh) llike.linearmodel(param=c(-2,2,1),x=treemass$logdbh,y=treemass$logdbh) debug(llike.linearmodel) llike.linearmodel(param=c(-2,2,1),x=treemass$logdbh,y=treemass$logdbh) param summary(x) SD summary(yhat) plot(x,y) Q llike.linearmodel(param=c(-2,2,1),x=treemass$logdbh,y=treemass$logagb) plot(x,y) points(x,yhat,col='red') length(llike) plot(x,llike) x11() points(x,yhat,col='red') plot(x,y) points(x,yhat,col='red') plot(y,yhat,col='red') plot(y,llike,col='red') total optim(par=c(10,-10),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb) optim(par=c(10,-10),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb) Q undebug(llike.linearmodel) optim(par=c(10,-10),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb) optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb) optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1))) coef(fit) llike.linearmodel(param=c(-2.14,2.51,.43),x=treemass$logdbh,y=treemass$logagb) llike.linearmodel(param=c(-2.14,2.51,1),x=treemass$logdbh,y=treemass$logagb) llike.linearmodel(param=c(-2.14,2.51,.1),x=treemass$logdbh,y=treemass$logagb) llike.linearmodel(param=c(-2.14,2.51,.42),x=treemass$logdbh,y=treemass$logagb) llike.linearmodel(param=c(-2.14,2.51,.41),x=treemass$logdbh,y=treemass$logagb) llike.linearmodel(param=c(-2.14,2.51,.44),x=treemass$logdbh,y=treemass$logagb) optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1))) plot(logagb~logdbh,data=treemass,cex=.5,pch=16) fit=lm(logagb~logdbh,data=treemass) sd(resid(fit)) abline(fit,col='red') summary(fit) optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1))) ?optim search() ls(5) head(cecrin) search() opt=optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1)),method='SANN') opt$par opt$par optim(par=c(10,-10),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb) fit1=optim(par=c(0,0,0),fn=llike.linearmodel,x=cerin$dbh1,y=cecrin$gr12) fit1=optim(par=c(0,0,0),fn=llike.linearmodel,x=cerin$dbh1,y=cecrin$gr12) fit1=optim(par=c(0,0,0),fn=llike.linearmodel.full,x=cerin$dbh1,y=cecrin$gr12) fit1=optim(par=c(0,0,0.1),fn=llike.linearmodel.full,x=cerin$dbh1,y=cecrin$gr12) fit1=optim(par=c(0,0,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12) fit1$par fit1=optim(par=c(0,0,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1))) fit1$par fit1$value plot(gr12~dbh1,data=cecrin) abline(fit1$par[1:2]) fit1=optim(par=c(0,-1,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1))) fit1$par abline(fit1$par[1:2]) fit1=optim(par=c(0,-1,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1)),method='SANN') fit1$par abline(fit1$par[1:2]) fit1$par fit1$value fit2=optim(par=fit1$par,fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1))) fit2$par fit2$value coef(lm(gr12~dbh1,data=cecrin)) abline(fit2$par[1:2],'red') abline(fit2$par[1:2],col='red') fit3=optim(par=c(fit2$par,0),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1))) fit3$value abline(fit3$par[1:2],col='blue') fit3$par abline(fit3$par[3:4],col='blue',lty='dashed') abline(sum(fit3$par[1:2],fit3$par[3:4]),col='blue') fit3$par[1:2]+fit3$par[3:4]) fit3$par[1:2]+fit3$par[3:4] fit3$par abline(fit3$par[1:2]+fit3$par[3:4],col='blue') plot(gr12~dbh1,data=cecrin) abline(fit3$par[1:2],col='blue') abline(fit3$par[1:2]+fit3$par[3:4],col='blue') abline(fit3$par[1:2]+2*fit3$par[3:4],col='blue') abline(fit3$par[1:2]-2*fit3$par[3:4],col='blue') fit2$par plot(gr12~dbh1,data=cecrin) abline(fit2$par[1:2],col='blue') abline(c(sum(fit2$par[c(1,3)]),fit2$par[2]),col='blue') abline(fit3$par[1:2]+2*fit3$par[3:4],col='blue') abline(fit3$par[1:2]+fit3$par[3:4],col='blue') system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops") source('~/meetings_workshops/Rmodeling/RfuncBayesWorkshop.r') source('~/meetings_workshops/Rmodeling/RcmdBayesWorkshop.r') Q graphics.off() sample=10 weights=rnorm(n=sample,mean=mu,sd=sigma) # If the mean were 125 and the SD 25, the likelihood is dnorm(x=weights,mean=mu,sd=sigma) # There are 10 points so 10 likelihoods, and it is easier to work with log-likelihood from the start. # We always need the sum of the log-likelihoods, which is the (log) likelihood of the entire set of data. The absolute # probability is not useful. loglike=dnorm(x=weights,mean=mu,sd=sigma,log=TRUE) totalLogLike=sum(loglike) cat(1) cat('1','\n') theta=0.76 N=10 dbinom(x=10,size=N,prob=theta) .76^10 dbinom(x=9,size=N,prob=theta) dbinom(x=0:10,size=N,prob=theta) dbinom(x=11,size=N,prob=theta) like=dbinom(x=11,size=N,prob=theta) sum(like) like=dbinom(x=0:10,size=N,prob=theta) sum(like) S=8 N dbinom(x=S,size=N,prob=.2) dbinom(x=S,size=N,prob=.3) dbinom(x=S,size=N,prob=.4) dbinom(x=S,size=N,prob=.9) dbinom(x=S,size=N,prob=.8) llike.binom=function(par,S,N) return( dbinom(x=S,size=N,prob=par,log=TRUE)) llike.binom ?optimize optimize(f=llike.binom,N=N,S=S) optimize(f=llike.binom,N=N,S=S,interval=c(0,1)) optimize(f=llike.binom,N=N,S=S,interval=c(0,1),maximum=TRUE) optimize(f=llike.binom,N=10,S=8,interval=c(0,1),maximum=TRUE) optimize(f=llike.binom,N=10,S=7,interval=c(0,1),maximum=TRUE) 7/10 theta=seq(0,1,by=0.01) postlike=dbinom(x=S,size=N,prob=theta) dev.set(2) plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob") theta S N postlike=dbinom(x=80,size=100,prob=theta) plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob") x11() postlike=dbinom(x=8,size=10,prob=theta) plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob") optimize(f=llike.binom,N=10,S=8,interval=c(0,1),maximum=TRUE) optimize(f=llike.binom,N=100,S=80,interval=c(0,1),maximum=TRUE) optimize(f=llike.binom,N=1000,S=800,interval=c(0,1),maximum=TRUE) postlike=dbinom(x=8,size=10,prob=theta) postlike=dbinom(x=800,size=1000,prob=theta) plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob") postlike=dbinom(x=8,size=10,prob=theta) plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob") theta=0.76 dbinom(x=S,size=N,prob=theta) theta=0.74 dbinom(x=S,size=N,prob=theta) theta=0.72 dbinom(x=S,size=N,prob=theta) theta=0.78 dbinom(x=S,size=N,prob=theta) theta=0.82 dbinom(x=S,size=N,prob=theta) S N dbinom(x=S,size=N,prob=.8) postlike=dbinom(x=S,size=N,prob=theta) sum(postlike) k=sum(0.01*postlike) post=postlike/k sum(0.01*post) lines(theta,post,type="l",col="red",main="adjusted binomial posterior",ylab="prob") dbinom(8,10,prob=.19) dbinom(8,10,prob=.20) dbinom(8,10,prob=.01) showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=100) showBinomialMCMC(thetastart=.2,S=8,N=10) Q showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=100) graphics.off() showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=1000) thetaMCMC=showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=1000,showstep=4000) summary(thetaMCMC) quantile(thetaMCMC,prob=.95) quantile(thetaMCMC,prob=.05) abline(v=.931,col='red') abline(v=.494,col='red') dev.set(3) abline(v=.931,col='red') abline(v=.494,col='red') sd(thetaMCMC) var(thetaMCMC) median(thetaMCMC) mean(thetaMCMC) length(thetaMCMC) length(whichMCMC<.6) length(which(thetaMCMC<.6)) graphics.off() fit3$value system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops") save.history(file='~/meetings_workshops/Rmodeling/Tupper2012/Tues7May.txt') savehistory(file='~/meetings_workshops/Rmodeling/Tupper2012/Tues7May.txt') args(metrop1step) args(linearfit.Bayes) source('~/meetings_workshops/Rmodeling/modelfitBayes.intro.r') args(linearfit.Bayes) linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=20,showstep=1,burnin=10) args(CI) CI source('Rut') CI linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=20,showstep=1,burnin=10) mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=20,showstep=1,burnin=10) mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=2000,showstep=100,burnin=1000) str(mod) mod$best mod$CI head(mod$full) head(mod$full,25) plot(mod[,1]) plot(mod$full[,1]) hist(mod$full[,1]) hist(mod$full[,1],breaks=25) quantile(mod$full[,1],prob=.95) quantile(mod$full[,1],prob=c(.05,.95)) mod$CI plot(mod$full[,1]) plot(mod$full[,2]) plot(mod$full[,3]) plot(mod$full[1001:2000,3]) quantile(mod$full[1001:2000,1],prob=c(.05,.95)) quantile(mod$full[1:100,1],prob=c(.05,.95)) mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=4000,showstep=100,burnin=1000) mod$CI mod$best colMeans(mod$full) colMeans(mod$full[1001:4000,]) plot(mod$full) plot(data.frame(mod$full)) plot(data.frame(mod$full[1001:4000,])) mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=4000,showstep=3000,burnin=1000) plot(data.frame(mod$full[1001:4000,])) mod$best colMeans(mod$full[1001:4000,]) quantile(mod$full[,1],prob=c(.025,.975)) head(mod$full) plot(gr12~dbh1,data=cecrin) abline(mod$best[1:2]) for(i in 1000:1100) abline(mod$full[i,1:2],col='gray') abline(mod$best[1:2]) for(i in 1000:4000) abline(mod$full[i,1:2],col='gray') abline(mod$best[1:2]) head(mod$full) source('~/meetings_workshops/Rmodeling/modelfitBayes.intro.r') source('~/meetings_workshops/Rmodeling/modelfitBayes.intro.r') mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=4000,showstep=3000,burnin=1000) oneparam args(llike.linearmodel.inter) llike.linearmodel.inter(inter=8.46,param=oneparam,x=x,y=y) Q mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=350,showstep=10,burnin=100) llike.linearmodel.inter(inter=8.46,param=oneparam,x=x,y=y) llike.linearmodel.inter(inter=9.46,param=oneparam,x=x,y=y) llike.linearmodel.inter(inter=7.46,param=oneparam,x=x,y=y) testllike=numeric() testpar1=1:20 for(i in 1:20) testllike[i]=llike.linearmodel.inter(inter=testpar1[i],param=oneparam,x=x,y=y) testllike plot(testpar1,testllike) oneparam Q source('~/meetings_workshops/Rmodeling/modelfitBayes.intro.r') mod=linearfit.Bayes(start=c(1,0,1),data=cecrin,xcol='dbh1',ycol='gr12',steps=350,showstep=10,burnin=100) n oneparam param[249,] param[248,] args(llike.linearmodel.inter) llike.linearmodel.inter(8.34,param=oneparam,x=x,y=y) llike.linearmodel.inter(oneparam[1],param=oneparam,x=x,y=y) llike.linearmodel.inter(9.34,param=oneparam,x=x,y=y) llike.linearmodel.inter(7.34,param=oneparam,x=x,y=y) testint=1:20 testlike=numeric() for(i in 1:20) testlike[i]=llike.linearmodel.inter(testint[i],param=oneparam,x=x,y=y) plot(testint,testlike) oneparam oneparam[2]=.1 oneparam for(i in 1:20) testlike[i]=llike.linearmodel.inter(testint[i],param=oneparam,x=x,y=y) plot(testint,testlike) testint=(-10):20 for(i in 1:20) testlike[i]=llike.linearmodel.inter(testint[i],param=oneparam,x=x,y=y) plot(testint,testlike) for(i in 1:31) testlike[i]=llike.linearmodel.inter(testint[i],param=oneparam,x=x,y=y) plot(testint,testlike) debug(llike.linearmodel.full) c n param SD summary(yhat) param plot(x,y) points(x,yhat) Q growth=growthModel.MultiSpecies() CTFSworkshop() growth=growthModel.MultiSpecies() source('~/meetings_workshops/Rmodeling/survivalGibbsHier.r') growth=growthModel.MultiSpecies() CTFSworkshop(data=TRYE) CTFSworkshop(data=TRUE) growth=growthModel.MultiSpecies() attach("~/meetings_workshops/Rmodeling/BCIlist.rdata") growth=growthModel.MultiSpecies() dim(growth) head(growth) hist(growth$Int) hist(growth$Int,breaks=25) hist(growth$Slope,breaks=25) hist(growth$Slope,breaks=25,xlim=c(-.2,.2)) hist(growth$Slope,breaks=25,xlim=c(-1,1)) hist(growth$Slope,breaks=250,xlim=c(-1,1)) hist(growth$Slope,breaks=250,xlim=c(-2,5)) plot(lgr12~sdbh1,data=BCIlist$tet2pa) for(i in 1:307) abline(growth[i,1:2]) for(i in 1:307) abline(drp(growth[i,1:2])) good=subset(growth,N>=10) good=matrix(subset(growth,N>=10)) good=subset(growth,N>=10) dim(good) for(i in 1:307) abline(drp(good[i,1:2])) head(growth,25) head(growth,20) summary(growth) hist(growth$Slope,breaks=250,xlim=c(-2,5)) hist(growth$Slope,breaks=25) hist(growth$Slope,breaks=250,xlim=c(-2,5)) plot(growth$Inter,growth$Slope) plot(growth$Inter,growth$Slope,xlim=c(-5,10)) plot(growth$Inter,growth$Slope,xlim=c(-5,10),ylim=c(-2,5)) plot(growth$Inter,growth$Slope,xlim=c(-5,10),ylim=c(-2,2)) plot(growth$Inter,growth$Slope) subset(growth,Inter>20) subset(growth,-Slope>15) savehistory(file='~/meetings_workshops/Rmodeling/Tupper2012/Wed8May.txt')