CTFSworkshop(data=TRUE)
grwfull300$dbhctr=(grwfull300$dbh1-mean(grwfull300$dbh1))
mod=lmer(LnGrowth~1+dbhctr+(1+dbhctr|Species),data=grwfull300)
xyplot(LnGrowth~dbhctr|Species,data=grwfull300)
dotplot(ranef(mod),postVar=TRUE)
dotplot(ranef(mod,postVar=TRUE))
xyplot(ht~dbh|species,data=treeht)
xyplot(lht~ldb|species,data=treeht)
xyplot(lht~ldb|species,data=treeht,
panel = function(x,y)
 {
panel.xyplot(x,y)
panel.abline(lm(y~x))
 }
)
head(treemass)
options(width=75)
head(treemass)
colnames(treemass)
mod=lmer(logagb~1+logdbh+(1+logdbh|Latin)+(1+logagb|ForestType),data=treemass)
treemass$ldbhctr=treemass$logdbh-mean(treemass$logdbh)
mod=lmer(logagb~1+ldbhctr+(1+ldbhctr|Latin)+(1+ldbhctr|ForestType),data=treemass)
fixef(mod)
coef(mod)$ForestType
head(coef(mod)$Latin)
colnames(treemass)
table(treemass$Locality)
mod=lmer(logagb~1+ldbhctr+(1+ldbhctr|Latin)+(1+ldbhctr|Locality),data=treemass)
display(mod)
head(coef(mod)$Local)
xyplot(logagb~logdbh|Locality,data=treemass)
xyplot(logagb~logddbh|ForestType,data=treemass)
xyplot(logagb~logdbh|ForestType,data=treemass)
xyplot(logagb~logdbh|ForestType,data=treemass,
   panel = function(x,y)
 {
panel.xyplot(x,y)
panel.abline(lm(y~x))
 }  
)
colnames(grwfull300)
options(width=60)
colnames(grwfull300)
mod=lmer(LnGrowth=1+dbhctr+(1+dbhctr|Species),data=grwfull300)
mod=lmer(LnGrowth~1+dbhctr+(1+dbhctr|Species),data=grwfull300)
fixef(mod)
dotplot(ranef(mod,postVar=TRUE))
dotplot(ranef(mod))
dotplot(ranef(mod,postVar=TRUE))
display(mod,digits=5)
dotplot(ranef(mod,postVar=TRUE),xlim=c(-.005,.005))
xyplot(LnGrowth~dbhctr|Species,data=grwfull300)
xyplot(ht~dbh|species,treeht)
xyplot(lht~ldb|species,treeht)
xyplot(lht~ldb|species,data=treeht, 
   panel = function(x,y)
 {
panel.xyplot(x,y)
panel.abline(lm(y~x))
 }  
   )
mod2=lmer(LnGrowth~1+dbhctr+(1+dbhctr|Species)+(1+dbhctr|census),data=grwfull300)
mod3=lmer(LnGrowth~1+dbhctr+(1|Species)+(1|census),data=grwfull300)
display(mod2,digits=5)
dotplot(ranef(mod2,postVar=TRUE))
dotplot(ranef(mod2,postVar=TRUE),xlim=c(-.005,.005))
dotplot(ranef(mod2,postVar=TRUE),xlim=c(-.001,.001))
coef(mod2)$census
dotplot(ranef(mod2,postVar=TRUE))
lm(LnGrowth~census,data=grwfull300)$coef
coef(mod2)$census
head(coef(mod2)$Species)
mod2=lmer(LnGrowth~1+dbhctr+(1+dbhctr|Species)+(1+dbhctr|census),data=grwfull300)
str(sleepstudy)
xyplot(Reaction~Days|Subject,data=sleepstudy)
plot(Reaction~Days,data=sleepstudy)
modslp=lmer(Reaction~1+Days+(1+Days|Subject),data=sleepstudy)
coef(modslp)$Sub
display(modslp)
colnames(treemass)
fit=lm(ldbhctr~logdbh,data=treemass)
aic(fit)
AIC(fit)
fit=lm(ldbhctr~1+logdbh,data=treemass)
AIC(fit)
fit=lm(ldbhctr~1+ldbhctr,data=treemass)
fit=lm(logagb~1+ldbhctr,data=treemass)
AIC(fit)
fit=lm(logagb~1+logdbh,data=treemass)
AIC(fit)
fit=lm(logagb~1+ldbhctr+I(ldbhctr^2),data=treemass)
fit=lm(logagb~1+ldbhctr,data=treemass)
fitq=lm(logagb~1+ldbhctr+I(ldbhctr^2),data=treemass)
AIC(fit,fitq)
mod=lmer(logagb~1+ldbhctr+I(ldbhctr^2)+(1+ldbhctr+I(ldbhctr^2)|Species),data=treemass)
mod=lmer(logagb~1+ldbhctr+I(ldbhctr^2)+(1+ldbhctr+I(ldbhctr^2)|species),data=treemass)
mod=lmer(logagb~1+ldbhctr+I(ldbhctr^2)+(1+ldbhctr+I(ldbhctr^2)|Latin),data=treemass)
AIC(fit,fitq,mod)
dim(subset(treemass,is.na(Latin)))
dim(subset(treemass,!is.na(Latin)))
mod=lmer(logagb~1+ldbhctr+(1+ldbhctr|Latin),data=treemass)
AIC(mod,fit)
fit=lm(logagb~1+ldbhctr,data=treemass)
mod=lmer(logagb~1+logdbh+(1+logdbh|Latin),data=treemass)
fit=lm(logagb~1+logdbh,data=treemass)
AIC(mod,fit)
fit=lm(logagb~1+ldbhctr,data=treemass)
mod=lmer(logagb~1+logdbh+(1+logdbh|Latin),data=treemass)
AIC(mod,fit)
treemassSp=subset(treemass,!is.na(Latin))
dim(treemass)
dim(treemassSp
)
fit=lm(logagb~1+logdbh,data=treemassSp)
mod=lmer(logagb~1+logdbh+(1+logdbh|Latin),data=treemassSp)
AIC(fit,mod)
plot(logagb~lodbh,data=treemassSp)
plot(logagb~logdbh,data=treemassSp)
modLoc=lmer(logagb~1+logdbh+(1+logdbh|Latin)+(1+logdbh|Locality),data=treemassSp)
anova(mod,modLoc)
fitFT=lm(logagb~1+logdbh+ForestType,data=treemassSp)
coef(fitFT)
fitFT=lm(logagb~1+logdbh+ForestType+logdbh*ForestType,data=treemassSp)
fitFT=lm(logagb~1+logdbh+ForestType,data=treemassSp)
fitFTInt=lm(logagb~1+logdbh+ForestType+logdbh*ForestType,data=treemassSp)
coef(fitFTInt)
cfdry=coef(fitFTInt)[1:2]
cfdr
cfdry
cfmoist=cfdry+coef(fitFTInt)[c(3,5)]
cfmoist
cfwet=cfdry+coef(fitFTInt)[c(4,6)]
plot(logagb~logdbh,data=treemassSp,cex=.3)
abline(cfdry,col='red',lwd=2)
abline(cfmoist,col='green',lwd=2)
abline(cfwet,col='blue',lwd=2)
modFT=lmer(logagb~1+logdbh+(1+logdbh|ForestType),data=treemassSp)
coef(modFT)
cfdry
cfmoist
cfwet
table(treemassSp$ForestType)
AIC(modFT,fitFTInt)
modmej=lmer(logagb~1+logdbh+ForestType+ForestType*logdbh+
(1+logdbh|Latin)+(1+logdbh|Locality),data=treemassSp)
(1+logdbh|Latin)+(1+logdbh|Locality),data=treemassSp,verbose=TRUE)
modmej=lmer(logagb~1+logdbh+ForestType+ForestType*logdbh+
(1+logdbh|Latin)+(1+logdbh|Locality),data=treemassSp,verbose=TRUE)
display(modmej)
hmod6=lmer(lht ~ ldb + I(ldb^2) + (1+ldb+I(ldb^2)|species), data=treeht)
str(treeht)
treeht$species=as.factor(treeht$species)
xyplot(lht~ldb|species, data=treeht, subset=species%in%levels(species)[1:25], ran=coef(hmod6)$species, fix = fixef(hmod6),
   panel = function(x, y, ran, fix)
 {
panel.xyplot(x,y)
x1 = seq(-3,3,0.5)
i=packet.number()
x2 = ran[i,1] * 1 + ran[i,2]*x1 + ran[i,3]*x1^2
x3 = fix[1] * 1 + fix[2]*x1 + fix[3]*x1^2
panel.lines(x2~x1, col='red')
panel.lines(x3~x1, lty=3)
     }  
  )
treemassSp=subset(treemass,!is.na(Latin))
savehistory('~/meetings_workshops/Rmodeling/Cenpat/history12OctAM.txt')
S=8
N=10
theta=0.76
dbinom(x=S,size=N,prob=theta)
dbinom(x=0:10,size=N,prob=theta)
sum(dbinom(x=0:10,size=N,prob=theta))
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.80
dbinom(x=S,size=N,prob=theta)
theta=0.82
dbinom(x=S,size=N,prob=theta)
S/N
theta=seq(0,1,by=0.01)
postlike=dbinom(x=S,size=N,prob=theta)
head(theta)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
dbinom(x=0:10,size=N,prob=theta)dbinom(x=0:10,size=N,prob=theta)
dbinom(x=0:10,size=N,prob=theta)dbinom(x=0:10,size=N,prob=theta)
plot(0:10,dbinom(x=0:10,size=N,prob=theta))
plot(0:10,dbinom(x=0:10,size=N,prob=.76))
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
plot(0:10,dbinom(x=0:10,size=N,prob=.76))
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
sum(0.01*postlike)
post=postlike/k
k=sum(0.01*postlike)
post=postlike/k
plot(theta,post,type='l')
lines(theta,postlike,lty='dotted')
betaprob=dbeta(theta,shape1=S+1,shape2=N-S+1)
points(theta,betaprob)
qbeta(p=c(0.025,0.975),shape1=S+1,shape2=N-S+1)
abline(v=qbeta(p=c(0.025,0.975),shape1=S+1,shape2=N-S+1))
betaprob=dbeta(theta,shape1=S+1,shape2=N-S+1)
S=80
N=100
betaprob=dbeta(theta,shape1=S+1,shape2=N-S+1)
plot(theta,betaprob,type='l')
abline(v=qbeta(p=c(0.025,0.975),shape1=S+1,shape2=N-S+1))
qbeta(p=c(0.025,0.975),shape1=S+1,shape2=N-S+1)
S=800
N=1000
qbeta(p=c(0.025,0.975),shape1=S+1,shape2=N-S+1)
plot(allS,like,pch=16,main="binomial likelihood")
allS=0:N
plot(allS,like,pch=16,main="binomial likelihood")
plot(allS,like,pch=16,main="binomial likelihood")like=dbinom(x=allS,size=N,prob=theta)
like=dbinom(x=allS,size=N,prob=theta)
plot(allS,like,pch=16,main="binomial likelihood")like=dbinom(x=allS,size=N,prob=theta)
plot(allS,like,pch=16,main="binomial likelihood")like=dbinom(x=allS,size=N,prob=theta)
plot(allS,like,pch=16,main="binomial likelihood")
like=dbinom(x=,size=N,prob=theta)
N=10
allS=0:N
like=dbinom(x=allS,size=N,prob=theta)
data.frame(allS,like)
plot(allS,like,pch=16,main="binomial likelihood")
theta=.76
allS=0:N
like=dbinom(x=allS,size=N,prob=theta)
data.frame(allS,like)
x11(height=5.5,width=7,xpos=1200,ypos=1)
# windows(height=5.5,width=7,xpos=1200,ypos=1)
plot(allS,like,pch=16,main="binomial likelihood")
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")
S=8
postlike=dbinom(x=S,size=N,prob=theta)
dev.set(2)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
k=sum(0.01*postlike)
k
1/k
post=postlike/k
abline(theta,post)
lines(theta,post)
showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=100)
x=showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=100)
showBinomialMCMC(thetastart=.2,S=8,N=10)
Q
betaprob=dbeta(theta,shape1=S+1,shape2=N-S+1)
betaprob=dbeta(theta,shape1=S+1,shape2=N-S+1)
theta=seq(0,1,by=0.01)
betaprob=dbeta(theta,shape1=S+1,shape2=N-S+1)
plot(theta,betaprob,type='l')
rbeta(n=5,shape1=S+1,shape2=N-S+1)
N
S
rbeta(n=5,shape1=S+1,shape2=N-S+1)
rbeta(n=5,shape1=S+1,shape2=N-S+1)
randdraw=rbeta(n=10000,shape1=S+1,shape2=N-S+1)
x11()
hist(randdraw,breaks=100)
betaprob[which(theta==.4)]
betaprob[which(theta==.39)]
betaprob[which(theta==.41)]
betaprob[which(theta==.42)]
showBinomialMCMC(thetastart=.2,S=8,N=10)
Q
showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE)
showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=1000)
runif(1)
runif(1)
runif(1)
runif(1)
?runif
runif(1,min=5,max=10)
runif(10,min=5,max=10)
runif(10)
graphics.off()
args(modelfitBayes)
args(modelfit.Bayes)
debug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
Q
debug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
Q
debug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
Q
undebug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
debug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
Q
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
Q
undebug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
debug(modelfit.Bayes)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
length(start)
start
head(param)
plot(x,y)
allparam
metropResult
head(param)
allparam
head(param)
head(param)
head(param)
Q
undebug(modelfit.Bayes)
head(param)
modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model)
fit=modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model,steps=2000,burnin=500,showstep=100)
names(fit)
head(fit$full)
plot(fit$full[,1])
plot(fit$full[,2])
fit$best
fit=modelfit.Bayes(start=c(0,0,1),xcol='logdbh',ycol='logagb',data=treemass, model=linear.model,steps=2000,burnin=1000,showstep=100)
fit$best
lm(logagb~logdbh,data=treemass)$coef
names(fit)
fit$CI
hist(fit$full[,2])
hist(fit$full[1001:2000,2])
plot(fit$full[,1])
fit$best
plot(fit$full[,1],fit$llike)
plot(fit$full[,1],fit$llike,ylim=c(-5000,0))
hist(llike[1001:2000])
hist(fit$llike[1001:2000])
plot(fit$full[,1],fit$llike,ylim=c(-1400,-1450))
plot(fit$full[,1],fit$llike,ylim=c(-1450,-1400))
plot(fit$full[,1],fit$llike,ylim=c(-15000,-1400))
plot(fit$full[,1],fit$llike,ylim=c(-1500,-1400))
plot(fit$full[1001:2000,1],fit$llike[1001:2000],ylim=c(-1500,-1400))
plot(fit$full[1001:2000,1],fit$llike[1001:2000],ylim=c(-1425,-1400))
plot(fit$full[1001:2000,1],fit$llike[1001:2000],ylim=c(-1425,-1410))
plot(fit$llike)
plot(fit$llike[1001:2000])
plot(fit$llike)
plot(fit$llike[1001:2000])
fit$best
plot(logagb~logdbh,data=treemass,cex=.3)
abline(fit$best[1:2],col='red',lwd=2)
for(i in 1900:2000) abline(fit$full[i,1:2],col='gray')
abline(fit$best[1:2],col='red',lwd=2)
for(i in 1000:2000) abline(fit$full[i,1:2],col='gray',lwd=.5)
abline(fit$best[1:2],col='red',lwd=2)
plot(logagb~logdbh,data=treemass,cex=.3,xlim=c(1.5,2))
plot(logagb~logdbh,data=treemass,cex=.3,xlim=c(1.5,2),ylimc-c(2,4))
plot(logagb~logdbh,data=treemass,cex=.3,xlim=c(1.5,2),ylim=c(2,4))
plot(logagb~logdbh,data=treemass,cex=.3,xlim=c(1.5,2),ylim=c(1.5,2.5))
plot(logagb~logdbh,data=treemass,cex=.3,xlim=c(1.5,2),ylim=c(1,3.5))
fit=modelfit.Bayes(start=c(0,0,1),xcol='momcat',ycol='Wt',data=pupsize, model=linear.model,steps=2000,burnin=1000,showstep=100)
fit$best
plot(Wt~momcat,data=pupsize)
abline(fit$best[1:2],col='red',lwd=2)
for(i in 1000:2000) abline(fit$full[i,1:2],col='gray',lwd=.5)
abline(fit$best[1:2],col='red',lwd=2)
for(i in 1000:2000) abline(fit$full[i,1:2],col='blue',lwd=.5)
abline(fit$best[1:2],col='red',lwd=2)
plot(Wt~momcat,data=pupsize,ylim=c(50,200))
abline(fit$best[1:2],col='red',lwd=1)
for(i in 1000:2000) abline(fit$full[i,1:2],col='blue',lwd=.5)
abline(fit$best[1:2],col='red',lwd=1)
fit$best
fit=modelfit.Bayes(start=c(0,0,1),xcol='momcat',ycol='Wt',data=pupsize, model=linear.model,steps=2000,burnin=1000,showstep=100)
plot(Wt~momcat,data=pupsize,ylim=c(50,200))
abline(fit$best[1:2],col='red',lwd=1)
fit$best
tail(fit$full)
for(i in 1000:2000) abline(fit$full[i,1:2],col='blue',lwd=.5)
abline(fit$best[1:2],col='red',lwd=1)
curve(linear.model(param=fit$best[1:2],x=x),add=TRUE,col='yellow')
for(i in 1500:200) curve(linear.model(param=fit$full[i,1:2],x=x),add=TRUE,col='orange')
for(i in 1500:2000) curve(linear.model(param=fit$full[i,1:2],x=x),add=TRUE,col='orange')
curve(linear.model(param=fit$best[1:2],x=x),add=TRUE,col='yellow')
plot(data.frame(fit$full[1001:2000,])
)
plot(Wt~momcat,data=pupsize,ylim=c(50,200))
curve(linear.model(param=fit$best[1:2],x=x),add=TRUE,col='yellow')
curve(linear.model(param=fit$best[1:2],x=x),add=TRUE,col='red')
for(i in 1500:2000) curve(linear.model(param=fit$full[i,1:2],x=x),add=TRUE,col='orange')
plot(data.frame(fit$full[1001:2000,]))
plot(Wt~momcat,data=pupsize,ylim=c(50,200))
plot(data.frame(fit$full[1001:2000,]))
savehistory('~/meetings_workshops/Rmodeling/Cenpat/history12OctPM.txt')