In this example, I will illustrate the use of the Berkson measurement error model on data from the BRFSS from Texas.
Traditionally, we only assume error in our outcome variables: e.g. yi /alpha+/beta∗xi+ei
But often, both y and x have errors in their measurement. This leads to a joint model for both x and y. Two basic forms of this include the Berkson and classical error models.
The Berkson model represents the situation where the observed predictors are less variable than their actual underlying true covariates, such as if you have a covariate that is recorded in intervals (income quartiles), vs the true continuous predictor (income in dollars). The former is less variable that the latter.
The classical model uses a latent variable specificaiton for the true underlying covariate. This is equivalent to stating that a variable x is an imperfect measure of the true covariate z.
We will use the BRFSS data for the state of Texas for our example, and use BMI as a continous outcome, and obesity status outcome (BMI >= 30) as a dichotomous outcome. To examine the Berkson model, I will use a categorized measure of poverty at the county level. To examine the classical error model, I will use the ACS estimates of the errors in the poverty rate estimates directly.
First we load our data and recode some variables:
library(rjags)## Loading required package: coda
## Loading required package: lattice
## Linked to JAGS 3.4.0
## Loaded modules: basemod,bugslibrary(dplyr)##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)
library(lme4)## Loading required package: Matrix
## Loading required package: Rcppload("~/Google Drive/dem7903_App_Hier/data/brfss_11.Rdata")
nams<-names(brfss_11)
newnames<-gsub("_", "", nams)
names(brfss_11)<-tolower(newnames)
brfss_11$statefip<-sprintf("%02d", brfss_11$state )
brfss_11$cofip<-sprintf("%03d", brfss_11$cnty )
brfss_11$cofips<-paste(brfss_11$statefip, brfss_11$cofip, sep="")
brfss_11$obese<-ifelse(brfss_11$bmi5/100 >=30, 1,0)
brfss_11$black<-recode(brfss_11$racegr2, recodes="2=1; 9=NA; else=0", as.factor.result=F)
brfss_11$white<-recode(brfss_11$racegr2, recodes="1=1; 9=NA; else=0", as.factor.result=F)
brfss_11$other<-recode(brfss_11$racegr2, recodes="3:4=1; 9=NA; else=0", as.factor.result=F)
brfss_11$hispanic<-recode(brfss_11$racegr2, recodes="5=1; 9=NA; else=0", as.factor.result=F)
#education level
brfss_11$lths<-recode(brfss_11$educa, recodes="1:3=1;9=NA; else=0", as.factor.result=F)
brfss_11$coll<-recode(brfss_11$educa, recodes="5:6=1;9=NA; else=0", as.factor.result=F)
brfss_11$agez<-scale(brfss_11$age, center=T, scale=T)
brfss_11$bmiz<-scale(brfss_11$bmi5/100, center=T, scale=T)
brfss_11$lowinc<-recode(brfss_11$incomg, recodes = "1:3=1; 4:5=0; else=NA")
#Read in ACS poverty rates at the county level
#poverty is measured with error, and we know the error in this case
acsecon<-read.csv("~/Google Drive/dem7903_App_Hier/data/aff_download/ACS_10_5YR_DP03_with_ann.csv")
acsecon$povrate<-acsecon[, "HC03_VC156"]
acsecon$poverr<-acsecon[, "HC04_VC156"]
#the next way I do it is the classical way, where you have a predictor that is given in ranges, not actual values.
#So I make grouped poverty rates using quantiles
quantile(acsecon$povrate, p=c(0, .25, .5, .75, 1))## 0% 25% 50% 75% 100%
## 0.0 7.4 10.4 14.1 44.9acsecon$pov_cut<-recode(acsecon$povrate, recodes = "0:7.4=1; 7.41:10.4=2; 10.41:14.1=3; 14.11:44.9=4")
acsecon$cofips<-substr(acsecon$GEO.id, 10,14)
acsecon<-acsecon[, c("cofips", "povrate", "poverr", "pov_cut")]
head(acsecon)## cofips povrate poverr pov_cut
## 1 01001 7.5 1.4 2
## 2 01003 9.1 1.0 2
## 3 01005 19.9 2.9 4
## 4 01007 9.4 2.3 2
## 5 01009 10.0 1.5 2
## 6 01011 22.6 5.6 4#and merge the data back to the brfss data
merged<-merge(x=brfss_11, y=acsecon, by.x="cofips", by.y="cofips", all.x=T)## Warning in merge.data.frame(x = brfss_11, y = acsecon, by.x = "cofips", :
## column names 'cholchk', 'mrace' are duplicated in the resultNext, I subset the data to select the observations from Texas.
tx<-subset(merged, subset = state=="48"& is.na(bmiz)==F& is.na(black)==F& is.na(lths)==F& is.na(lowinc)==F)
nwncos<-table(tx$cofips)
nwncos #Number of people within counties##
## 48029 48133 48157 48201 48303 48329 48423 48439 48453
## 863 498 772 1185 603 421 452 483 859tx$conum<-rep(1:length(unique(tx$cofips)), nwncos[nwncos!=0])
length(unique(tx$conum)) #Number of counties## [1] 9The first model will ignore the errors in the predictor.
model1<-"
model{
for (i in 1:n){
bmi[i]~dnorm(mu[i], tau)
mu[i]<-b[1]*black[i]+b[2]*hisp[i]+b[3]*other[i]+b[4]*lths[i]+b[5]*coll[i]+u[cofips[i]]
}
for (j in 1:ncos)
{
u[j]~dnorm(muu[j], tau_u)
muu[j]<-u0+gam*(pov[j]-mean(pov[]))
}
#priors
u0~dnorm(0, .01)
for(j in 1:5) { b[j]~dnorm(0, .01)}
gam~dnorm(0, .01)
tau<-pow(sd, -2)
sd~dunif(0,100)
tau_u<-pow(sd_u, -2)
sd_u~dunif(0,100)
}
"
#get initial values from a linear model
b<-as.numeric(coef(lm(bmiz~black+hispanic+other+lths+coll+scale(povrate, center=T), tx)))
b## [1] 0.02931728 0.50613065 0.17991108 -0.15825431 0.10769865 -0.06863772
## [7] 0.03491620#make the data
dat<-list(bmi=as.numeric(tx$bmiz), obese=tx$obese, black=tx$black, hisp=tx$hispanic, other=tx$other, lths=tx$lths, coll=tx$coll, age=as.numeric(tx$agez), lowinc=tx$lowinc, pov=as.numeric(tapply(tx$povrate, tx$cofips, mean)), n=length(tx$bmiz),cofips=tx$conum, ncos=length(unique(tx$cofips)))
init.rng1<-list(".RNG.seed" = 1234, ".RNG.name" = "base::Mersenne-Twister", u0=b[1], b=b[2:6], gam=b[7])
init.rng2<-list(".RNG.seed" = 5678, ".RNG.name" = "base::Mersenne-Twister", u0=b[1], b=b[2:6], gam=b[7])
#start the model
mod1<-jags.model(file=textConnection(model1), data=dat,inits =list(init.rng1, init.rng2) , n.chains=2)## Warning in jags.model(file = textConnection(model1), data = dat, inits =
## list(init.rng1, : Unused variable "obese" in data## Warning in jags.model(file = textConnection(model1), data = dat, inits =
## list(init.rng1, : Unused variable "age" in data## Warning in jags.model(file = textConnection(model1), data = dat, inits =
## list(init.rng1, : Unused variable "lowinc" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 43133
##
## Initializing model#next, we update the model, this is the "burn in" period
update(mod1, 20000)
#sample samples from each chain, thinning every 25th iteration
samps<-coda.samples(mod1, variable.names=c("u0", "gam", "b", "sd", "sd_u"), n.iter=2000, n.thin=25)
effectiveSize(samps)## b[1] b[2] b[3] b[4] b[5] gam sd
## 1968.8843 1233.3564 3138.9558 1033.6617 225.0964 975.5451 2567.3602
## sd_u u0
## 175.0144 202.3266gelman.diag(samps, multivariate = F)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## b[1] 1.00 1.00
## b[2] 1.00 1.00
## b[3] 1.00 1.02
## b[4] 1.00 1.00
## b[5] 1.00 1.01
## gam 1.00 1.00
## sd 1.00 1.00
## sd_u 1.01 1.03
## u0 1.00 1.01#Numerical summary of each parameter, here I also include the 90% credible interval:
summary(samps, quantiles = c(.025, .975))##
## Iterations = 21001:23000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## b[1] 0.50987 0.045039 0.0007121 0.0010216
## b[2] 0.17939 0.038017 0.0006011 0.0011476
## b[3] -0.15389 0.055694 0.0008806 0.0009948
## b[4] 0.11227 0.051677 0.0008171 0.0016181
## b[5] -0.06265 0.032617 0.0005157 0.0021903
## gam 0.01329 0.007990 0.0001263 0.0002623
## sd 1.00169 0.009091 0.0001437 0.0001794
## sd_u 0.04018 0.027738 0.0004386 0.0021597
## u0 0.02518 0.034242 0.0005414 0.0024216
##
## 2. Quantiles for each variable:
##
## 2.5% 97.5%
## b[1] 0.4230474 0.5977690
## b[2] 0.1062045 0.2558031
## b[3] -0.2640193 -0.0454144
## b[4] 0.0139422 0.2165754
## b[5] -0.1275151 -0.0004197
## gam -0.0030158 0.0295003
## sd 0.9842130 1.0198156
## sd_u 0.0004914 0.1027521
## u0 -0.0383918 0.0944713dic.samples(mod1, n.iter = 1000, type = "pD")## Mean deviance: 17434
## penalty 11.9
## Penalized deviance: 17445Here, I have error in the poverty rate because the rate has been measured in an interval fashion, instead of a continuous fashion. One thing you may encounter in this type of situation is where you have a calibration estimate for the intervaled data, such as a regression of the true value on the intervaled value, so you know the parameters from such an equation. I estimate these calibration values below using a linear model of the intervaled estimate on the true estimate.
model2<-"
model{
for (i in 1:n){
bmi[i]~dnorm(mu[i], tau)
mu[i]<-b[1]*black[i]+b[2]*hisp[i]+b[3]*other[i]+b[4]*lths[i]+b[5]*coll[i]+u[cofips[i]]
}
for (j in 1:ncos)
{
u[j]~dnorm(muu[j], tau_u)
muu[j]<-u0+gam*(z[j])
z[j]~dnorm(muz[j], tauz)
muz[j]<-alphaz+bz*povcut[j]
}
#priors
for(j in 1:5) { b[j]~dnorm(0, .01)}
gam~dnorm(0,.01)
u0~dnorm(0, .0001)
tau<-pow(sd, -2)
sd~dunif(0,100)
tau_u<-pow(sd_u, -2)
sd_u~dunif(0,100)
}
"
#get calibration values for error model.
summary(lm(povrate~pov_cut, tx))##
## Call:
## lm(formula = povrate ~ pov_cut, data = tx)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5502 -0.3620 0.1557 0.8439 1.0498
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.36798 0.04104 82.06 <2e-16 ***
## pov_cut 3.09407 0.01467 210.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9128 on 6134 degrees of freedom
## Multiple R-squared: 0.8788, Adjusted R-squared: 0.8788
## F-statistic: 4.447e+04 on 1 and 6134 DF, p-value: < 2.2e-16#Make the data
dat<-list(bmi=as.numeric(tx$bmiz), obese=tx$obese, black=tx$black, hisp=tx$hispanic, other=tx$other, lths=tx$lths, coll=tx$coll, age=as.numeric(tx$agez), lowinc=tx$lowinc, povcut=tapply(tx$pov_cut, tx$cofips, median ), n=length(tx$bmiz),cofips=tx$conum, ncos=length(unique(tx$cofips)), alphaz=3.367, bz=3.094, tauz=1/(.9128^2))
#get initial values
b<-as.numeric(fixef(lmer(bmiz~black+hispanic+other+lths+coll+pov_cut+(1|cofips), tx)))
b## [1] -0.04088886 0.50858872 0.17811710 -0.15601629 0.11105274 -0.06681006
## [7] 0.02686786VarCorr(lmer(bmiz~black+hispanic+other+lths+coll+pov_cut+(1|cofips), tx))## Groups Name Std.Dev.
## cofips (Intercept) 0.046201
## Residual 1.001549#initialize
init.rng1<-list( b=b[2:6], u0=b[1], gam=b[7], sd=1.001549, sd_u=.046201)
init.rng2<-list( b=b[2:6], u0=b[1], gam=b[7], sd=1.001549, sd_u=.046201)
#Start the model
mod2<-jags.model(file=textConnection(model2), data=dat,inits =list(init.rng1, init.rng2) , n.chains=2)## Warning in jags.model(file = textConnection(model2), data = dat, inits =
## list(init.rng1, : Unused variable "obese" in data## Warning in jags.model(file = textConnection(model2), data = dat, inits =
## list(init.rng1, : Unused variable "age" in data## Warning in jags.model(file = textConnection(model2), data = dat, inits =
## list(init.rng1, : Unused variable "lowinc" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 43143
##
## Initializing model#next, we update the model, this is the "burn in" period
update(mod2, 30000)
#sample samples from each chain, thinning every 25th iteration
samps2<-coda.samples(mod2, variable.names=c("u0", "b","gam", "sd", "sd_u", "z", "sd", "sd_u"), n.iter=2000, n.thin=25)## Warning in jags.samples(model, variable.names, n.iter, thin, type = "trace", : Failed to set trace monitor for sd
## Monitor already exists and cannot be duplicated## Warning in jags.samples(model, variable.names, n.iter, thin, type = "trace", : Failed to set trace monitor for sd_u
## Monitor already exists and cannot be duplicatedeffectiveSize(samps2)## b[1] b[2] b[3] b[4] b[5] gam
## 2138.60942 1748.75166 2911.11641 803.02335 249.60109 68.92078
## sd sd_u u0 z[1] z[2] z[3]
## 2405.04009 261.94953 68.54170 2681.46889 2990.04280 2309.37881
## z[4] z[5] z[6] z[7] z[8] z[9]
## 3124.52492 2670.47419 2536.43120 1820.33461 2411.54651 1517.70422gelman.diag(samps2, multivariate = F)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## b[1] 1.00 1.00
## b[2] 1.00 1.02
## b[3] 1.00 1.00
## b[4] 1.00 1.00
## b[5] 1.00 1.01
## gam 1.04 1.06
## sd 1.00 1.00
## sd_u 1.00 1.00
## u0 1.03 1.06
## z[1] 1.00 1.00
## z[2] 1.00 1.00
## z[3] 1.00 1.01
## z[4] 1.00 1.01
## z[5] 1.00 1.01
## z[6] 1.00 1.00
## z[7] 1.00 1.01
## z[8] 1.00 1.00
## z[9] 1.00 1.01#Numerical summary of each parameter, here I also include the 90% credible interval:
summary(samps2, quantiles = c(.025, .975))##
## Iterations = 31001:33000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## b[1] 0.51106 0.044334 0.0007010 0.0009696
## b[2] 0.17915 0.037841 0.0005983 0.0009042
## b[3] -0.15364 0.054971 0.0008692 0.0010195
## b[4] 0.11343 0.053182 0.0008409 0.0018869
## b[5] -0.06245 0.033981 0.0005373 0.0021496
## gam 0.01189 0.010112 0.0001599 0.0013474
## sd 1.00162 0.008988 0.0001421 0.0001833
## sd_u 0.05399 0.032977 0.0005214 0.0020772
## u0 -0.11208 0.122702 0.0019401 0.0157495
## z[1] 12.79380 0.902000 0.0142619 0.0188566
## z[2] 15.68188 0.884428 0.0139840 0.0162108
## z[3] 6.44646 0.902590 0.0142712 0.0191741
## z[4] 12.58502 0.885850 0.0140065 0.0158862
## z[5] 12.72573 0.912685 0.0144308 0.0184077
## z[6] 9.68164 0.921951 0.0145773 0.0188387
## z[7] 12.50257 0.922895 0.0145923 0.0225793
## z[8] 9.75881 0.914454 0.0144588 0.0187690
## z[9] 12.34640 0.946765 0.0149697 0.0252258
##
## 2. Quantiles for each variable:
##
## 2.5% 97.5%
## b[1] 0.424325 0.59770
## b[2] 0.108111 0.25275
## b[3] -0.262130 -0.04670
## b[4] 0.006049 0.21667
## b[5] -0.130353 0.00225
## gam -0.008807 0.03153
## sd 0.984054 1.01956
## sd_u 0.006818 0.13507
## u0 -0.356820 0.14689
## z[1] 11.062171 14.55954
## z[2] 13.979474 17.37629
## z[3] 4.691435 8.21881
## z[4] 10.832859 14.38868
## z[5] 10.959678 14.46606
## z[6] 7.846071 11.49392
## z[7] 10.724212 14.33990
## z[8] 7.968003 11.56446
## z[9] 10.485059 14.25753dic.samples(mod2, n.iter = 1000, type = "pD")## Mean deviance: 17434
## penalty 12.48
## Penalized deviance: 17446Here the poverty rate is treated as a latent variable with known error (from the ACS margin of error). Alternatively, you could put a prior on the latent variable variance.
model3<-"
model{
for (i in 1:n){
bmi[i]~dnorm(mu[i], tau)
mu[i]<-b[1]*black[i]+b[2]*hisp[i]+b[3]*other[i]+b[4]*lths[i]+b[5]*coll[i]+u[cofips[i]]
}
for (j in 1:ncos)
{
u[j]~dnorm(muu[j], tau_u)
muu[j]<-u0+gam*(povtrue[j]-mean(povtrue[]))
povrate[j]~dnorm(povtrue[j],taupov[j] )
povtrue[j]~dnorm(0, taup)T(0,100)
taupov[j]<-pow(poverr[j],-2)
}
#priors
u0~dnorm(0, .01)
for(j in 1:5) { b[j]~dnorm(0, .01)}
gam~dnorm(0, .01)
tau<-pow(sd, -2)
sd~dunif(0,100)
tau_u<-pow(sd_u, -2)
sd_u~dunif(0,100)
taup~dgamma(.01, .01)
}
"
dat<-list(bmi=as.numeric(tx$bmiz), obese=tx$obese, black=tx$black, hisp=tx$hispanic, other=tx$other, lths=tx$lths, coll=tx$coll, age=as.numeric(tx$agez), lowinc=tx$lowinc, povrate=as.numeric(tapply(tx$povrate, tx$cofips, mean)), poverr=tapply(tx$poverr/1.645, tx$cofips, mean), n=length(tx$bmiz),cofips=tx$conum, ncos=length(unique(tx$cofips)))
#get initial values
b<-as.numeric(fixef(lmer(bmiz~black+hispanic+other+lths+coll+scale(povrate, center=T)+(1|cofips), tx)))
b## [1] 0.03040590 0.50906581 0.17779706 -0.15478146 0.11036624 -0.06592581
## [7] 0.03446866VarCorr(lmer(bmiz~black+hispanic+other+lths+coll+scale(povrate, center=T)+(1|cofips), tx))## Groups Name Std.Dev.
## cofips (Intercept) 0.035714
## Residual 1.001545init.rng1<-list( u0=b[1], b=b[2:6], gam=b[7])
init.rng2<-list( u0=b[1], b=b[2:6], gam=b[7])
mod3<-jags.model(file=textConnection(model3), data=dat,inits =list(init.rng1, init.rng2) , n.chains=2)## Warning in jags.model(file = textConnection(model3), data = dat, inits =
## list(init.rng1, : Unused variable "obese" in data## Warning in jags.model(file = textConnection(model3), data = dat, inits =
## list(init.rng1, : Unused variable "age" in data## Warning in jags.model(file = textConnection(model3), data = dat, inits =
## list(init.rng1, : Unused variable "lowinc" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 43158
##
## Initializing model#next, we update the model, this is the "burn in" period
update(mod3, 20000)
#sample 1000 samples from each chain, thinning every 25th iteration
samps3<-coda.samples(mod3, variable.names=c("u0", "b","gam", "sd", "sd_u", "povtrue", "povrate"), n.iter=2000, n.thin=25)
effectiveSize(samps3)## b[1] b[2] b[3] b[4] b[5] gam
## 2437.5915 1568.5121 3013.7831 1065.8596 341.2738 1425.7628
## povrate[1] povrate[2] povrate[3] povrate[4] povrate[5] povrate[6]
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## povrate[7] povrate[8] povrate[9] povtrue[1] povtrue[2] povtrue[3]
## 0.0000 0.0000 0.0000 3872.4346 2049.5267 3807.1865
## povtrue[4] povtrue[5] povtrue[6] povtrue[7] povtrue[8] povtrue[9]
## 4000.0000 3556.0354 2703.2562 2918.1094 4000.0000 3761.9943
## sd sd_u u0
## 2547.5197 380.4240 342.8834gelman.diag(samps3, multivariate = F)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## b[1] 1.00 1.00
## b[2] 1.00 1.00
## b[3] 1.00 1.00
## b[4] 1.00 1.00
## b[5] 1.00 1.01
## gam 1.00 1.00
## povrate[1] NaN NaN
## povrate[2] NaN NaN
## povrate[3] NaN NaN
## povrate[4] NaN NaN
## povrate[5] NaN NaN
## povrate[6] NaN NaN
## povrate[7] NaN NaN
## povrate[8] NaN NaN
## povrate[9] NaN NaN
## povtrue[1] 1.00 1.00
## povtrue[2] 1.00 1.00
## povtrue[3] 1.00 1.00
## povtrue[4] 1.00 1.00
## povtrue[5] 1.00 1.00
## povtrue[6] 1.00 1.00
## povtrue[7] 1.00 1.01
## povtrue[8] 1.00 1.00
## povtrue[9] 1.00 1.00
## sd 1.00 1.00
## sd_u 1.01 1.06
## u0 1.00 1.00#Numerical summary of each parameter:
summary(samps3, quantiles = c(.025, .975))##
## Iterations = 21001:23000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## b[1] 0.50856 0.045210 0.0007148 0.0009170
## b[2] 0.17741 0.037301 0.0005898 0.0009550
## b[3] -0.15491 0.055716 0.0008810 0.0010226
## b[4] 0.10834 0.050594 0.0008000 0.0015581
## b[5] -0.06702 0.031530 0.0004985 0.0017093
## gam 0.01372 0.008967 0.0001418 0.0002375
## povrate[1] 13.20000 0.000000 0.0000000 0.0000000
## povrate[2] 15.90000 0.000000 0.0000000 0.0000000
## povrate[3] 6.10000 0.000000 0.0000000 0.0000000
## povrate[4] 13.70000 0.000000 0.0000000 0.0000000
## povrate[5] 12.30000 0.000000 0.0000000 0.0000000
## povrate[6] 10.10000 0.000000 0.0000000 0.0000000
## povrate[7] 11.30000 0.000000 0.0000000 0.0000000
## povrate[8] 10.40000 0.000000 0.0000000 0.0000000
## povrate[9] 11.10000 0.000000 0.0000000 0.0000000
## povtrue[1] 13.20419 0.183786 0.0029059 0.0029572
## povtrue[2] 15.20978 1.828042 0.0289039 0.0406905
## povtrue[3] 6.10252 0.300538 0.0047519 0.0048759
## povtrue[4] 13.68942 0.183155 0.0028959 0.0028962
## povtrue[5] 12.33436 0.532059 0.0084126 0.0089287
## povtrue[6] 10.20148 0.907504 0.0143489 0.0175023
## povtrue[7] 11.21962 0.588596 0.0093065 0.0109008
## povtrue[8] 10.40920 0.178102 0.0028160 0.0028163
## povtrue[9] 11.05568 0.300233 0.0047471 0.0049099
## sd 1.00186 0.009315 0.0001473 0.0001846
## sd_u 0.04662 0.028250 0.0004467 0.0014430
## u0 0.03107 0.034948 0.0005526 0.0018998
##
## 2. Quantiles for each variable:
##
## 2.5% 97.5%
## b[1] 0.416773 0.594974
## b[2] 0.102011 0.249291
## b[3] -0.263098 -0.044232
## b[4] 0.009332 0.207569
## b[5] -0.128957 -0.006447
## gam -0.003974 0.031573
## povrate[1] 13.200000 13.200000
## povrate[2] 15.900000 15.900000
## povrate[3] 6.100000 6.100000
## povrate[4] 13.700000 13.700000
## povrate[5] 12.300000 12.300000
## povrate[6] 10.100000 10.100000
## povrate[7] 11.300000 11.300000
## povrate[8] 10.400000 10.400000
## povrate[9] 11.100000 11.100000
## povtrue[1] 12.838419 13.566852
## povtrue[2] 11.574008 18.802500
## povtrue[3] 5.511780 6.692156
## povtrue[4] 13.331411 14.040259
## povtrue[5] 11.274370 13.354635
## povtrue[6] 8.401609 11.953778
## povtrue[7] 10.083676 12.359108
## povtrue[8] 10.061355 10.756515
## povtrue[9] 10.486926 11.637176
## sd 0.984061 1.020340
## sd_u 0.008032 0.116552
## u0 -0.035262 0.100764dat$povrate## [1] 13.2 15.9 6.1 13.7 12.3 10.1 11.3 10.4 11.1dat$poverr## 48029 48133 48157 48201 48303 48329 48423
## 0.1823708 2.0060790 0.3039514 0.1823708 0.5471125 0.9118541 0.6079027
## 48439 48453
## 0.1823708 0.3039514dic.samples(mod3, n.iter = 1000, type = "pD")## Mean deviance: 17444
## penalty 20.8
## Penalized deviance: 17464In this case, the first model, without any measurement error fits best, followed by the berkson model, then the classical model, using the DIC.