# Use a Gibbs sampler to fit a linear model of y on one variable x, with either constant SD or a linearly increasing SD. # data is dataframe including x (predictor), y (observations) # Initial parameters are start -- # intercept then slope # the third parameter is SD if constant # third and fourth are intercept and slope of SD if SD is variable linearfit.Bayes=function(start,data,xcol="dbh1",ycol="status2",steps=20,showstep=1,burnin=500) { no.param=length(start) param=matrix(nrow=steps,ncol=no.param) param[1,]=start stepsize=rep(1,no.param) llike=numeric() x=data[,xcol] y=data[,ycol] ## Calculate likelihood of data at starting point for(i in 2:steps) { ### Update the intercept oneparam=param[i-1,] metropResult=metrop1step(func=llike.linearmodel.inter,start.param=param[i-1,1],adjust=1.01,target=0.25, param=oneparam,x=x,y=y,scale.param=stepsize[1]) param[i,1]=metropResult[1] stepsize[1]=metropResult[2] # browser() ### Update the slope oneparam=c(param[i,1],param[i-1,2:no.param]) metropResult=metrop1step(func=llike.linearmodel.slope,start.param=param[i-1,2],adjust=1.01,target=0.25, param=oneparam,x=x,y=y,scale.param=stepsize[2]) param[i,2]=metropResult[1] stepsize[2]=metropResult[2] ### Update the first SD parameter oneparam=c(param[i,1:2],param[i-1,3:no.param]) metropResult=metrop1step(func=llike.linearmodel.SD1,start.param=param[i-1,3],adjust=1.01,target=0.25, param=oneparam,x=x,y=y,scale.param=stepsize[3]) param[i,3]=metropResult[1] stepsize[3]=metropResult[2] ### If there are four parameters, update the second SD parameter if(no.param==4) { oneparam=c(param[i,1:3],param[i-1,4]) metropResult=metrop1step(func=llike.linearmodel.SD2,start.param=param[i-1,4],adjust=1.01,target=0.25, param=oneparam,x=x,y=y,scale.param=stepsize[4]) param[i,4]=metropResult[1] stepsize[4]=metropResult[2] } # browser() ### Calculate and save likelihood after finishing each step llike[i]=metropResult[4] #### cat results if desired if(i%%showstep==0) cat("Step ", i, ": Param: ", round(param[i,],3), "llike: ", round(llike[i],4), "\n") } keep=(-1)*(1:burnin) best=colMeans(param[keep,]) CI=apply(param[keep,],2,CI) return(list(best=best,CI=CI,full=param,llike=llike)) } # Likelihood function for a linear model, where the first argument is the intercept parameter alone. llike.linearmodel.inter=function(inter,param,x,y) { param[1]=inter return(llike.linearmodel.full(param=param,x=x,y=y)) } # Likelihood function for a linear model, where the first argument is the slope parameter alone. llike.linearmodel.slope=function(slope,param,x,y) { param[2]=slope return(llike.linearmodel.full(param=param,x=x,y=y)) } # Likelihood function for a linear model, where the first argument is the first SD parameter alone. llike.linearmodel.SD1=function(SD1,param,x,y) { param[3]=SD1 return(llike.linearmodel.full(param=param,x=x,y=y)) } # Likelihood function for a linear model, where the first argument is the second SD parameter alone. llike.linearmodel.SD2=function(SD2,param,x,y) { param[4]=SD2 return(llike.linearmodel.full(param=param,x=x,y=y)) }