turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

- Home
- /
- Analytics
- /
- Stat Procs
- /
- PROC MCMC error

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

01-16-2014 09:35 PM

Hello All;

I am trying to create a Bayesian Logistic Hierarchical using proc MCMC. I am using the exact dataset and codes provided by SAS in the example in the following link,

the only difference between my model and the example model is that the response variable in the example is a continues variable and has a normal distribution while in my model the response is a dichotomous variable and thus has a binary distribution.

I create a new variable in the data-set, called overweight which is equal to 1 if the Rat's weight is above average and zero otherwise and then try to model the probability of being over weight as a logistic function of age.

data rats;

array days[5] (8 15 22 29 36);

input weight @@;

subject = ceil(_n_/5);

index = mod(_n_-1, 5) + 1;

age = days[index];

drop index days:;

datalines;

151 199 246 283 320 145 199 249 293 354

147 214 263 312 328 155 200 237 272 297

135 188 230 280 323 159 210 252 298 331

141 189 231 275 305 159 201 248 297 338

177 236 285 350 376 134 182 220 260 296

160 208 261 313 352 143 188 220 273 314

154 200 244 289 325 171 221 270 326 358

163 216 242 281 312 160 207 248 288 324

142 187 234 280 316 156 203 243 283 317

157 212 259 307 336 152 203 246 286 321

154 205 253 298 334 139 190 225 267 302

146 191 229 272 302 157 211 250 285 323

132 185 237 286 331 160 207 257 303 345

169 216 261 295 333 157 205 248 289 316

137 180 219 258 291 153 200 244 286 324

;

run;

**data rats;**

**set rats;**

**if weight>242 then overweight=1; else overweight=0;**

**run;**

**proc mcmc data=rats nmc=10000 outpost=postout**

**seed=17 init=random;**

**ods select Parameters REParameters PostSummaries;**

**array theta[2] alpha beta;**

**array theta_c[2];**

**array Sig_c[2,2];**

**array mu0[2];**

**array Sig0[2,2];**

**array S[2,2] (0.02 0 0 20);**

**begincnst;**

**call zeromatrix(mu0);**

**call identity(Sig0);**

**call mult(Sig0, 1000, Sig0);**

**endcnst;**

**parms theta_c Sig_c {121 0 0 0.26} var_y;**

**prior theta_c ~ mvn(mu0, Sig0);**

**prior Sig_c ~ iwish(2, S);**

**prior var_y ~ igamma(0.01, scale=0.01);**

**random theta ~ mvn(theta_c, Sig_c) subject=subject;**

**mu = logistic(alpha + beta * age+var_y );**

**model overweight~ binary(mu);**

**run;**

When I run the model, I get this error message,

**Observation 2 of the response variable has a value of 0, and it does not produce a positive log-likelihood value. Possible causes are: the response variable is outside of the support set of the likelihood function or some distribution parameters**

Would you please help me with this issue?

are missing.

Accepted Solutions

Solution

01-17-2014
02:25 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

01-17-2014 02:25 PM

Well, the big difference is in the distributions, so when you say "the only difference between my model and the example model is that the response variable in the example is a continues variable and has a normal distribution while in my model the response is a dichotomous variable and thus has a binary distribution", you are making a ginormous change. I'll through two things out, neither of which is going to get this to start converging ( I don't think). First, since the Monte Carlo part of the procedure depends on the parameter starting values, make sure that the variance and covariance estimates refer to the logits of the data. Using the values from the normal distribution won't suffice.

Second, why? Why dichoomize the data when you have a very good analysis treating the response variable as continuous? Dichotomization throws away a LOT of information, and immediately reduces the power of the analysis. For example, see the paper presented by Dr. Stephen Senn at the 2011 FDA/industry Statistics Workshop, in the short course, Statistical Issues in Drug Development. From Section 3, slides 29 and 30:

Slide 29:

Dichotomisation

Prospects for a Cure

- I am pessimistic
- Most physicians seem happy with dichotomisation
- Most statisticians seem happy to indulge them
- "That's not my dpartment' syndrome

- We have to bring home the following message:

Slide 30

**DICHOTOMISATION IS VERY SILLY!**

Your problem points out one of the drawbacks Dr. Senn was trying to make.

Steve Denham

All Replies

Solution

01-17-2014
02:25 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

01-17-2014 02:25 PM

Well, the big difference is in the distributions, so when you say "the only difference between my model and the example model is that the response variable in the example is a continues variable and has a normal distribution while in my model the response is a dichotomous variable and thus has a binary distribution", you are making a ginormous change. I'll through two things out, neither of which is going to get this to start converging ( I don't think). First, since the Monte Carlo part of the procedure depends on the parameter starting values, make sure that the variance and covariance estimates refer to the logits of the data. Using the values from the normal distribution won't suffice.

Second, why? Why dichoomize the data when you have a very good analysis treating the response variable as continuous? Dichotomization throws away a LOT of information, and immediately reduces the power of the analysis. For example, see the paper presented by Dr. Stephen Senn at the 2011 FDA/industry Statistics Workshop, in the short course, Statistical Issues in Drug Development. From Section 3, slides 29 and 30:

Slide 29:

Dichotomisation

Prospects for a Cure

- I am pessimistic
- Most physicians seem happy with dichotomisation
- Most statisticians seem happy to indulge them
- "That's not my dpartment' syndrome

- We have to bring home the following message:

Slide 30

**DICHOTOMISATION IS VERY SILLY!**

Your problem points out one of the drawbacks Dr. Senn was trying to make.

Steve Denham