Hello again, WinBugs users!
I am trying to model a data generation process in which no observations on
out-of-sample i's are available, and the sample inclusion is
non-deterministic (more poor-quality social science data). Bloom and
Killingsworth (1985) propose a stochastic truncation model in which an
unobserved random threshold variable Ci and the dependent variable Yi are
distributed bivariate normal with parameters mu Y mu Y sigma2 Y sigma2 C
and sigma YC. This is not implementable as written in WinBugs, of course,
because of the limitation on missing data in a multivariate node (and all
Ci's are missing). DeGroot (1986) points out that it is possible to
calculate a univariate truncated distribution that takes into account the
relationship between Y and C:
f(yi|ci>0,...) = normal(y[i]|mu,sigma2) * Pr(C[i]>0|y[i]) / Pr(C[i] >0)
The tricky part is Pr(C[i]>0|y[i]) which is equal to:
1-CDFnormal(0|D[i],P[i]) where
D[i] = muC[i] + (sigmayc/sigma2y)*(y[i] - muY[i]) and
P[i] = sigma2c - (sigma2yc/sigma2y) [Note: sigma2c is set equal to 1 to
identify/scale the model]
I have been trying, without success, to implement this function in WinBugs
using the "ones trick." Currently, Bugs is reporting an undefined variable,
which I am unable to locate. My code, complete with 10 fake observations on
the covariates, is included below. Any help or hints would be greatly
appreciated,as always!
Thanks,
___________________
Jack Buckley
SUNY Stony Brook
Department of Political Science
[log in to unmask]
(631) 632-4353
---------------------
WinBugs code follows
---------------------
Model
Cbar <- mean(C)
Ybar <- mean(Y)
for(i in 1:N){
product[i] <- (Y[i]-Ybar)*(C[i]-Cbar)
d[i] <- muC[i] + ((sigmayc/pow(sigmay,2))*(Y[i]-muY[i]))
muY[i] <-alpha1 + beta1 * x1[i] + beta2 * x2[i] + beta3 * x3[i]
muC[i] <-alpha2 + gamma1 * x1[i] + gamma2 * x2[i]
C[i] ~ dnorm(0,1)
comp1[i] ~ dnorm(muY[i], tau1)
temp1[i] ~ dnorm(d[i], p)
temp2[i] ~ dnorm(muC[i], 1)
comp2[i] <- pow((1-temp1[i]),-1)
comp3[i] <- pow((1-temp2[i]),-1)
ones[i] <- 1
ones[i] ~ dbern( prob[i] )
prob[i] <- comp1[i]*comp2[i]/comp3[i] / K
}
p <- 1 - (pow(sigmayc,2)/pow(sigmay,2))
alpha ~ dnorm(0, 1.0E-6)
beta1 ~ dnorm(0, 1.0E-6)
beta2 ~ dnorm(0, 1.0E-6)
beta3 ~ dnorm(0, 1.0E-6)
alpha2 ~ dnorm(0, 1.0E-6)
gamma1 ~ dnorm(0, 1.0E-6)
gamma2 ~ dnorm(0, 1.0E-6)
tau1 ~ dgamma(.001, .001)
sigmay <- 1/sqrt(tau1)
productsum <- sum(product)
sigmayc <-mean(productsum)
K <-100
list(Y=c(1,2,3,4,5,6,7,8,9,10),
x1=c(4,6,6,8,9,10,10,14,12,17),
x2=c(0,0,1,0,0,1,1,1,0,1),
x3=c(100,80,78,60,58,43,33,23,30,10),
C=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
comp1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
temp1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
temp2=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
N=10)
list(alpha1=0, alpha2=0, beta1=0, beta2=0, beta3=0, gamma1=0, gamma2=0,
tau1 =1)
-------------------------------------------------------------------
To mail the BUGS list, mail to [log in to unmask]
You can search old messages at www.jiscmail.ac.uk/lists/bugs.html
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
This list is for discussion of modelling issues and the BUGS
software. For help with crashes and error messages, first mail
[log in to unmask]
