To be clear everything "runs" with no error message... the only hint of a problem is at the end of the code: the plot will not fill out/ it is empty.
if anyone has any idea why something like this might happen, i would greatly appreciate it... so i can handle it quickly. thanks in advance. On Mar 28, 2013, at 7:55 PM, Nicole Ford wrote: > i am having problem running my own data. yesterday it was working just fine. > today it is not. this is the code i was using as an example to follow. > this code ALSO worked just fine yesterday, and is no longer working at all. > i suspect it is a problem with either my computer or the software, at this > point. if THIS won't even run.... something is wrong. > > i can assure you this isn't HW.... i know dave, but i am no longer at UW-M > and i have never learned HLMs and i am learning this on my own for my own > research. > > his code is here, along with data. it is short, quick, etc. > > http://www.quantoid.net/936/Lecture7.R > > ### R code from vignette source 'Lecture7.Rnw' > > ################################################### > ### code chunk number 1: opts > ################################################### > options(useFancyQuotes=F) > > > ################################################### > ### code chunk number 2: data1 > ################################################### > library(foreign) > therms <- na.omit(read.dta("http://quantoid.net/936/2008_difftherm.dta";)) > unstate <- unique(therms[,1]) > therms$numstate <- match(therms$state, unstate) > library(runjags) > dat <- dump.format(list( > N = nrow(therms), J=length(unstate), > y = therms$difftherm, > numstate = therms$numstate > )) > > > ################################################### > ### code chunk number 3: exchange > ################################################### > exchange.mod <- "model{ > for(i in 1:N){ > y[i] ~ dnorm(mu, tau) > } > mu ~ dnorm(0,.001) > tau ~ dgamma(.1,.1) > }" > exchange.out <- run.jags(exchange.mod, > data=dat, burnin=10000, sample=50000, > thin=5, monitor=c("mu", "tau"), > monitor.deviance=T, monitor.pd=T, > silent.jags=T) > > > > ################################################### > ### code chunk number 4: exchange > ################################################### > FE.mod <- "model{ > for(i in 1:N){ > y[i] ~ dnorm(mu[numstate[i]], tau[numstate[i]]) > } > for(j in 1:J){ > mu[j] ~ dnorm(0,.001) > tau[j] ~ dgamma(.1,.1) > } > }" > FE.out <- run.jags(FE.mod, > data=dat, burnin=10000, sample=50000, > thin=5, monitor=c("mu", "tau"), > monitor.deviance=T, monitor.pd=T, > silent.jags=T) > > > ################################################### > ### code chunk number 5: exchange > ################################################### > hier.mod <- "model{ > for(i in 1:N){ > y[i] ~ dnorm(mu[numstate[i]], tau[numstate[i]]) > } > for(j in 1:J){ > mu[j] ~ dnorm(theta,nu) > tau[j] ~ dgamma(a,b) > } > theta ~ dnorm(0,.01) > nu ~ dgamma(.1,.1) > a ~ dunif(0,1000) > b ~ dunif(0,1000) > }" > hier.out <- run.jags(hier.mod, > data=dat, burnin=10000, sample=100000, > thin=10, monitor=c("mu", "tau", "theta", "nu", "a", "b"), > monitor.deviance=T, monitor.pd=T, > silent.jags=T) > > > ################################################### > ### code chunk number 6: sums > ################################################### > hier.chains <- combine.mcmc(hier.out$mcmc) > FE.chains <- combine.mcmc(FE.out$mcmc) > exchange.chains <- combine.mcmc(exchange.out$mcmc) > > mu.bar <- apply(FE.chains[, grep("mu\\[", colnames(FE.chains))], 2, mean) > mu.bar2 <- apply(hier.chains[, grep("mu\\[", colnames(hier.chains))], 2, mean) > ns <- aggregate(therms$numstate, list(therms$stateabb), length) > plot(mu.bar, mu.bar2, cex=sqrt(ns[,2])/3, > xlab = "FE mu[j]", > ylab = "Hierarchical mu[j]") > abline(a=0, b=1) > > > ################################################### > ### code chunk number 7: dotchart > ################################################### > fe.mu <- FE.chains[,grep("mu\\[", colnames(FE.chains))] > fe.ci <- t(apply(fe.mu, 2, quantile, c(.5,.025,.975))) > rownames(fe.ci) <- unstate > fe.ci <- fe.ci[order(fe.ci[,1]), ] > dotchart(fe.ci[order(fe.ci[,1]),1], lcolor="white", pch=16, > xlim=range(c(fe.ci))) > segments(fe.ci[,2], 1:34, fe.ci[,3], 1:34) > mu.ci <- quantile(exchange.chains[,1], c(.5,.025,.975)) > polygon(x=mu.ci[c(2,3,3,2)], > y = c(-1,-1,36,36), > col=rgb(128,128,128,100, maxColorValue=255), > border=NA) > abline(v=mu.ci[1], lty=2, lwd=2) > axis(4, at=1:34, labels=ns[match(rownames(fe.ci), ns[,1]),2], > cex.axis=.75, las=2) > > > ################################################### > ### code chunk number 8: femeans > ################################################### > library(sm) > sm.density(mu.bar, model="normal") > > > ############################ > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
