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.
