############################################ # STAT 740 R examples # Tim Hanson # U. South Carolina # October 15, 2017 ############################################ ############################################ # Independence M-H for Challenger data ############################################ y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0) t=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,73,75,75,76,76,78,79,81) X=cbind(rep(1,23),t) # design matrix # log-posterior function m=matrix(c(168.75,-2.402,-2.402,0.03453),2,2) m.inv=solve(m) lp=function(theta){ o=sum(X%*%theta*y-log(1+exp(X%*%theta)))-0.5*t(theta)%*%m.inv%*%theta o[1,1] } f=glm(y~t,family="binomial") bi=200; gi=10200 b=matrix(0,gi+bi,2) # keep MCMC iterates in b colnames(b)=c("beta0","beta1") v0=1.1*vcov(f); v0sqrt=t(chol(v0)); v0inv=solve(v0) b[1,]=f$coef; m=1; accept=0 llp=lp(b[1,]) # last log-posterior repeat{ m=m+1 bn=as.vector(t(f$coef+v0sqrt%*%rnorm(2))) nlp=lp(bn) lrho=nlp-llp-0.5*(b[m-1,]-f$coef)%*%v0inv%*%(b[m-1,]-f$coef)+ 0.5*(bn-f$coef)%*%v0inv%*%(bn-f$coef) if(log(runif(1))gi){break} } cat("Acceptance rate=",accept/gi,"\n") library(coda) b=mcmc(data=b,start=bi+1,end=gi,thin=1) summary(b) plot(b) autocorr.plot(b) # mixing? cumuplot(b) ############################################ # Random walk M-H for Challenger data ############################################ y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0) t=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,73,75,75,76,76,78,79,81) X=cbind(rep(1,23),t) # design matrix # log-posterior function m=matrix(c(168.75,-2.402,-2.402,0.03453),2,2) m.inv=solve(m) lp=function(theta){ o=sum(X%*%theta*y-log(1+exp(X%*%theta)))-0.5*t(theta)%*%m.inv%*%theta o[1,1] } f=glm(y~t,family="binomial") bi=200; gi=10200 b=matrix(0,gi+bi,2) # keep MCMC iterates in b colnames(b)=c("beta0","beta1") v0=0.7*vcov(f); v0sqrt=t(chol(v0)); v0inv=solve(v0) b[1,]=f$coef; m=1; accept=0 llp=lp(b[1,]) # last log-posterior repeat{ m=m+1 bn=as.vector(t(b[m-1,]+v0sqrt%*%rnorm(2))) nlp=lp(bn) lrho=nlp-llp if(log(runif(1))gi){break} } cat("Acceptance rate=",accept/gi,"\n") library(coda) b=mcmc(data=b,start=bi+1,end=gi,thin=1) summary(b) plot(b) autocorr.plot(b) # mixing? cumuplot(b) ############################################ # Gibbs sampling: normal data by hand # note: no tuning parameter w/ Gibbs sampler ############################################ # generate some data n=200; mu=5; sigma=2 x=rnorm(n,mu,sigma) # prior a=0.0001; b=0.0001 p=0.0001; mu.m=0 # Gibbs sampler bi=200; gi=5200 b=matrix(0,gi+bi,2) # keep MCMC iterates in b colnames(b)=c("mu","tau") b[1,]=c(mean(x),1/var(x)); m=1 xbar=mean(x); s2=sum((x-xbar)^2)/n; repeat{ m=m+1 v=1/(b[m-1,2]*n+p) b[m,1]=rnorm(1,v*(b[m-1,2]*n*xbar+p*mu.m),sqrt(v)) b[m,2]=rgamma(1,a+0.5*n,rate=b+0.5*n*(s2+(b[m,1]-xbar)^2)) if(floor(m/100)==(m/100)){cat("it=",m," mu=",b[m,1]," tau=",b[m,2],"\n")} if(m>gi){break} } library(coda) b=mcmc(data=b,start=bi+1,end=gi,thin=1) summary(b) hist(1/sqrt(b[,2]),freq=F,xlab="sigma",ylab="posterior") autocorr.plot(b) # mixing? cumuplot(b) ############################################ # Gibbs sampling: normal data via JAGS # need to install JAGS first ############################################ library(R2jags) # starting values for JAGS start=list(list("mu"=xbar,"tau"=1/s2)) model=function(){ # model for JAGS for(i in 1:n){ x[i]~dnorm(mu,tau) } mu~dnorm(0,0.0001) tau~dgamma(0.0001,0.0001) sigma<-1/sqrt(tau) } dat=list("x","n") # data for JAGS par=c("mu","sigma") f=jags(data=dat,inits=start,parameters.to.save=par, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=5) print(f) traceplot(f) plot(f) autocorr.plot(f) # mixing? ############################################ # Finite mixture of normals hand-coded # uses latent z's ############################################ library(MASS) # has galaxies data set library(MCMCpack) # has rdirichlet x=galaxies/1000 # x=data, J=# components, gi=number of MCMC iterates # mu ~ N(m,1/p), tau ~ Gamma(a,b), pi ~ Dirich(w,...,w) # bi=burnin # no "error trapping" -- sometimes doesn't work... fm=function(x,J,gi,bi,m,p,a,b,w){ n=length(x); nc=1:J; z=1:n delta=(max(x)-min(x))/(J+1) pi=rep(1/J,J); mu=min(x)+delta*1:J; sd=rep(delta/2,J) o=list(pi=matrix(0,gi,J),mu=matrix(0,gi,J),sd=matrix(0,gi,J)) gibbs=0 repeat{ gibbs=gibbs+1 for(i in 1:n){z[i]=sample(J,1,prob=pi*dnorm(x[i],mu,sd))} for(j in 1:J){ nc[j]=sum(z==j) v=1/(p+nc[j]/sd[j]^2) mu[j]=rnorm(1,v*(p*m+nc[j]*mean(x[z==j])/sd[j]^2),sqrt(v)) sd[j]=1/sqrt(rgamma(1,a+0.5*nc[j],rate=b+0.5*sum((x[z==j]-mu[j])^2))) } pi=rdirichlet(1,nc+w)[1,] if(ceiling(gibbs/1000)==gibbs/1000){cat("it=",gibbs,"\n")} if(gibbs>gi+bi){break} if(gibbs>bi){ o$mu[gibbs-bi,]=mu o$sd[gibbs-bi,]=sd o$pi[gibbs-bi,]=pi } } o } f=fm(x,3,10000,2000,mean(x),0.001,0.001,0.001,1) J=3 # number of components hist(x,freq=F,breaks=10) r=max(x)-min(x) g=seq(min(x)-r/4,max(x)+r/4,length=1000) de=Vectorize(function(x){mean(rowSums(dnorm(x,f$mu,f$sd)*f$pi))}) lines(g,de(g),col="blue") ############################################ # Finite mixture of normals in JAGS # uses latent z's ############################################ library(R2jags) library(MASS) # has galaxies data set x=galaxies/1000 m=mean(x); n=length(x) J=3 # starting values for JAGS start=list(list("mu"=rep(mean(x),J),"tau"=rep(1/var(x),J), "pi"=rep(1/J,J),"z"=sample(J,n,replace=T))) model=function(){ # model for JAGS for(i in 1:n){ x[i]~dnorm(mu[z[i]],tau[z[i]]) z[i]~dcat(pi[1:J]) } for(j in 1:J){ mu[j]~dnorm(m,0.001) tau[j]~dgamma(0.001,0.001) alpha[j]<-1 } pi[1:J] ~ ddirch(alpha[1:J]) } dat=list("x","n","J","m") # data for JAGS par=c("mu","tau","pi") f=jags(data=dat,inits=start,parameters.to.save=par, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=5) print(f) traceplot(f) plot(f) # can add more if not converged! # discards the the entire first run f=update(f,n.iter=10000) f=as.mcmc(f)[[1]] # use first chain (only chain!) f=list(mu=f[,2:(J+1)],pi=f[,(J+2):(2*J+1)],sd=1/sqrt(f[,(2*J+2):(3*J+1)])) hist(x,freq=F,breaks=10) r=max(x)-min(x) g=seq(min(x)-r/4,max(x)+r/4,length=1000) de=Vectorize(function(x){mean(rowSums(dnorm(x,f$mu,1/sqrt(f$sd))*f$pi))}) lines(g,de(g),col="blue") ############################################ # Finite mixture of normals in JAGS # uses "zeros trick" # no latent z's ############################################ library(R2jags) library(MASS) # has galaxies data set x=galaxies/1000 m=mean(x); n=length(x); J=3 y=rep(0,n) # data to use in "trick" # starting values for JAGS start=list(list("mu"=rep(mean(x),J),"tau"=rep(1/var(x),J),"pi"=rep(1/J,J))) model=function(){ # model for JAGS for(i in 1:n){ for(j in 1:J){ to.add[i,j]<-dnorm(x[i],mu[j],1/sqrt(tau[j]))*pi[j] } L[i]<- -log(sum(to.add[i,1:J])) y[i]~dpois(L[i]) } for(j in 1:J){ mu[j]~dnorm(m,0.001) tau[j]~dgamma(0.001,0.001) alpha[j]<-1 } pi[1:J] ~ ddirch(alpha[1:J]) } dat=list("x","n","J","m","y") # data for JAGS par=c("mu","tau","pi") f=jags(data=dat,inits=start,parameters.to.save=par, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=5) print(f) traceplot(f) plot(f) f=as.mcmc(f)[[1]] # use first chain (only chain!) f=list(mu=f[,2:(J+1)],pi=f[,(J+2):(2*J+1)],sd=1/sqrt(f[,(2*J+2):(3*J+1)])) hist(x,freq=F,breaks=10) r=max(x)-min(x) g=seq(min(x)-r/4,max(x)+r/4,length=1000) de=Vectorize(function(x){mean(rowSums(dnorm(x,f$mu,1/sqrt(f$sd))*f$pi))}) lines(g,de(g),col="blue") ############################################ # Challenger data via JAGS ############################################ y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0) t=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,73,75,75,76,76,78,79,81) f=glm(y~t,family="binomial") library(R2jags) n=length(y) # starting values for JAGS start=list(list("beta"=f$coef)) model=function(){ # model for JAGS for(i in 1:n){ y[i]~dbern(pi[i]) pi[i]<-ilogit(beta[1]+beta[2]*t[i]) } for(i in 1:2){ beta[i]~dnorm(0,0.0001) } OR<-exp(beta[2]) } chal=list("y","t","n") # data for JAGS param=c("beta","OR") f=jags(data=chal,inits=start,parameters.to.save=param, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=1) print(f) traceplot(f) b=as.mcmc(f)[[1]] library(coda) summary(b) plot(b) autocorr.plot(b) # mixing? Better than random-walk M-H! cumuplot(b) ############################################ # Ache hunting via JAGS ############################################ library(R2jags) d=read.table("http://people.stat.sc.edu/hansont/stat740/ache.txt",header=T) attach(d) # makes 'age', 'days', and 'monkeys' available to R n=length(age) id=1:47 m=monkeys; a=age; t=days f=glmer(monkeys~age+(1|id),family="poisson",offset=log(days),nAGQ=100) # starting values for JAGS start=list(list("beta"=c(0,0),"tau"=1,"u"=rep(0,n))) start=list(list("beta"=f@beta,"tau"=1/f@theta^2,"u"=f@u)) model=function(){ # model for JAGS for(i in 1:n){ m[i]~dpois(mu[i]) mu[i]<-t[i]*exp(beta[1]+beta[2]*a[i]+u[i]) u[i]~dnorm(0,tau) } for(i in 1:2){ beta[i]~dnorm(0,0.0001) } tau~dgamma(0.0001,0.0001) sigma<-1/sqrt(tau) OR<-exp(beta[2]) } kills=list("m","a","t","n") # data for JAGS param=c("beta","sigma","OR") f=jags(data=kills,inits=start,parameters.to.save=param, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=5) print(f) traceplot(f) plot(f) ############################################ # predicted number of kills from new hunter # aged 50 years over 5-day trek ############################################ model=function(){ # model for JAGS for(i in 1:n){ m[i]~dpois(mu[i]) mu[i]<-t[i]*exp(beta[1]+beta[2]*a[i]+u[i]) u[i]~dnorm(0,tau) } for(i in 1:2){ beta[i]~dnorm(0,0.0001) } tau~dgamma(0.0001,0.0001) sigma<-1/sqrt(tau) u.new~dnorm(0,tau) mu.new<-5*exp(beta[1]+beta[2]*50+u.new) y.new~dpois(mu.new) } param=c("beta","sigma","mu.new","y.new") f=jags(data=kills,inits=start,parameters.to.save=param, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=5) print(f) traceplot(f) fm=as.mcmc(f) summary(fm) plot(fm) library(lattice) densityplot(fm) xyplot(fm) geweke.diag(fm) raftery.diag(fm) geweke.plot(fm) ############################################ # Adaptive M-H for Challenger data ############################################ y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0) t=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,73,75,75,76,76,78,79,81) X=cbind(rep(1,23),t) # design matrix # log-posterior function m=matrix(c(168.75,-2.402,-2.402,0.03453),2,2) m.inv=solve(m) lp=function(theta){ o=sum(X%*%theta*y-log(1+exp(X%*%theta)))-0.5*t(theta)%*%m.inv%*%theta o[1,1] } f=glm(y~t,family="binomial") bi=200; gi=5200 b=matrix(0,gi+bi,2) # keep MCMC iterates in b colnames(b)=c("beta0","beta1") b[1,]=f$coef; mm=b[1,]; m=1; accept=0 v0=vcov(f); vm=matrix(0,2,2) llp=lp(b[1,]) # last log-posterior repeat{ m=m+1 if(m>bi){ cov=vm } else{ cov=v0 } bn=as.vector(t(b[m-1,]+t(chol(cov))%*%rnorm(2))) nlp=lp(bn) if(log(runif(1))<(nlp-llp)){ b[m,]=bn llp=nlp accept=accept+1 } else{ b[m,]=b[m-1,] } mo=mm # old mean for recursive variance update mm=((m-1)*mm+b[m,])/m # new mean vm=((m-1)*(vm+mo%*%t(mo))+b[m,]%*%t(b[m,]))/m-mm%*%t(mm) # new var if(floor(m/1000)==(m/1000)){cat("it=",m," beta=",b[m,],"\n")} if(m>gi){break} } cat("Acceptance rate=",accept/gi,"\n") library(coda) b=mcmc(data=b,start=bi+1,end=gi,thin=1) summary(b) plot(b) ############################################ # Adaptive M-H for Challenger data # using package "adaptMCMC" # EASY!!! Just need log-posterior, # starting values, and crude scale matrix ############################################ library(adaptMCMC) af=MCMC(lp,n=5000,init=f$coef,scale=vcov(f), adapt=TRUE,acc.rate=0.25) plot(af$samples) af=convert.to.coda(af) summary(af) plot(af) autocorr.plot(af) cumuplot(af) ############################################ # Gibbs sampling: censored normal data ############################################ library(MASS) # has VA data library(survival) # has survfit & survreg # censored data histogram chist=function(x,d,J){ a=min(x-(max(x)-min(x))/10000); b=max(x); J=10 delta=(b-a)/J t=seq(a,b,length=(J+1)) f=summary(survfit(Surv(x,d)~1),times=t) y=c(0,rep(f$surv[1:J]-f$surv[2:(J+1)],1,each=2),0)/delta x=rep(t,1,each=2) plot(x,y,type="l",col="blue",xlab="censored data", ylab="density",ylim=c(0,max(y*1.1))) } VAs=subset(VA,prior==0) attach(VAs) f=survreg(Surv(log(stime),status)~1,dist="gaussian") # the easy way! chist(stime,status,10) # histogram x=seq(0,600,length=1000) lines(x,dlnorm(x,f$coef,f$scale),col="red") # plot of log-normal fit summary(f) # MLEs of mu and sigma; compare to values you get via E-M x=log(stime) n=length(x) # prior a=0.0001; b=0.0001 p=0.0001; mu.m=0 cens.index=which(status==0) cens=x[cens.index] num.cens=length(cens) # Gibbs sampler bi=200; gi=5200 b=matrix(0,gi+bi,2) # keep MCMC iterates in b colnames(b)=c("mu","tau") b[1,]=c(mean(x),1/var(x)); m=1 repeat{ m=m+1 # update right-censored log-survival times for(j in 1:num.cens){ # update right-censored values u=runif(1,pnorm(cens[j],b[m-1,1],1/sqrt(b[m-1,2]))) x[cens.index[j]]=qnorm(u,b[m-1,1],1/sqrt(b[m-1,2])) } xbar=mean(x); s2=sum((x-xbar)^2)/n; v=1/(b[m-1,2]*n+p) b[m,1]=rnorm(1,v*(b[m-1,2]*n*xbar+p*mu.m),sqrt(v)) b[m,2]=rgamma(1,a+0.5*n,rate=b+0.5*n*(s2+(b[m,1]-xbar)^2)) if(floor(m/100)==(m/100)){cat("it=",m," mu=",b[m,1]," tau=",b[m,2],"\n")} if(m>gi){break} } library(coda) b=mcmc(data=b,start=bi+1,end=gi,thin=1) summary(b) plot(b) ############################################ # JAGS: censored normal data ############################################ library(MASS) # has VA data VAs=subset(VA,prior==0) attach(VAs) # x should have NA for censored values # cens should have Inf for not-censored x=log(stime); n=length(x) cens=rep(Inf,n) cens[which(status==0)]=x[which(status==0)] x[which(status==0)]=NA ind=1-status # 1 if right-censored x.start=rep(NA,n) x.start[which(status==0)]=cens[which(status==0)]+0.01 xbar=mean(x); s2=var(x) # starting values for JAGS start=list(list("mu"=xbar,"tau"=1/s2,"x"=x.start)) model=function(){ # model for JAGS for(i in 1:n){ x[i]~dnorm(mu,tau) ind[i]~dinterval(x[i],cens[i]) } # only part of model that's different from above is 'dinterval' mu~dnorm(0,0.0001) tau~dgamma(0.0001,0.0001) sigma<-1/sqrt(tau) } dat=list("x","n","cens","ind") # data for JAGS par=c("mu","sigma") f=jags(data=dat,inits=start,parameters.to.save=par, n.chains=1,n.iter=21000,n.burnin=1000,model.file=model, n.thin=5) print(f) traceplot(f) plot(f) autocorr.plot(f) # mixing?