########################################################################

# Logit models.

######################################################################

model

{

for(i in 1:N)

{

beta[i, 1:3] ~ dmnorm(mu.beta[], R[,])

for(j in 1:T[i])

{

y[i,j] <- ty[j,i]

n[i,j] <- tn[j,i]

x[i,j] <- tx[j,i]

y[i,j] ~ dbin(p[i,j], n[i,j]) # likelihood

logit(p[i,j]) <- beta[i, 1] + beta[i,2]*x[i,j]

+ beta[i, 3]*x[i,j]*x[i,j]

}

}

mu.beta[1:3] ~ dmnorm(mean[1:3], prec[1:3, 1:3])

R[1:3, 1:3] ~ dwish(Omega[1:3,1:3], 3)

# prediction of peak ages

for(i in 1:N)

{

peak.age[i] <- -0.5*beta[i, 2]/beta[i,3]

}

}

# Initial values

list(mu.beta = c( -7.69, 0.350, -0.0058),

beta = structure(

.Data = c( -7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058,

-7.69, 0.350, -0.0058),

.Dim =c(10,3)

),

R = structure(

.Data = c(0.1, 0, 0,

0, 0.1, 0,

0, 0, 0.1),

.Dim =c(3,3)

)

)

# Data 1

list(T = c(23, 13, 22, 18, 22, 22, 22, 22, 18, 17), N = 10,

mean = c(0,0,0),

Omega = structure(

.Data = c(0.1, 0, 0,

0, 0.1, 0,

0, 0, 0.1),

.Dim =c(3,3)

),

prec = structure(

.Data = c(1.0E-6, 0, 0,

0, 1.0E-6, 0,

0, 0, 1.0E-6),

.Dim =c(3,3)

)

)

# Data 2

# ty = transpose(y) -> y = transpose(ty)

# Aaron, Greenberg, Killebrew, Mantle, Mays, McCovey, Ott, Ruth, Schmidt, Sosa

ty[,1] ty[,2] ty[,3] ty[,4] ty[,5] ty[,6] ty[,7] ty[,8] ty[,9] ty[,10]

13 0 0 13 20 13 0 0 1 4

27 12 4 23 4 13 1 4 18 15

26 26 5 21 41 18 18 3 36 10

44 36 2 27 51 20 42 2 38 8

30 1 0 37 36 44 25 11 38 33

39 40 42 52 35 18 29 29 38 25

40 58 31 34 29 39 38 54 21 36

34 33 46 42 34 36 23 59 45 40

45 41 48 31 29 31 35 35 48 36

44 2 45 40 40 36 31 41 31 66

24 13 49 54 49 45 33 46 35 63

32 44 25 30 38 39 31 25 40 50

44 25 39 15 47 18 36 47 36 64

39 0 44 35 52 14 27 60 33 49

29 0 17 19 37 29 19 54 37 40

44 0 49 23 22 22 27 46 35 35

38 0 41 22 23 23 30 49 12 14

47 0 28 18 13 7 18 46 6 0

34 0 26 0 28 28 26 41 0 0

40 0 5 0 18 12 21 34 0 0

20 0 13 0 8 15 1 22 0 0

12 0 14 0 6 1 0 6 0 0

10 0 0 0 0 0 0 0 0 0

END

# Data 3

# tx = transpose(x) -> x = transpose(tx)

tx[,1] tx[,2] tx[,3] tx[,4] tx[,5] tx[,6] tx[,7] tx[,8] tx[,9] tx[,10]

20 19 18 19 20 21 17 19 22 20

21 22 19 20 21 22 18 20 23 21

22 23 20 21 23 23 19 21 24 22

23 24 21 22 24 24 20 22 25 23

24 25 22 23 25 25 21 23 26 24

25 26 23 24 26 26 22 24 27 25

26 27 24 25 27 27 23 25 28 26

27 28 25 26 28 28 24 26 29 27

28 29 26 27 29 29 25 27 30 28

29 30 27 28 30 30 26 28 31 29

30 34 28 29 31 31 27 29 32 30

31 35 29 30 32 32 28 30 33 31

32 36 30 31 33 33 29 31 34 32

33 0 31 32 34 34 30 32 35 33

34 0 32 33 35 35 31 33 36 34

35 0 33 34 36 36 32 34 37 35

36 0 34 35 37 37 33 35 38 36

37 0 35 36 38 38 34 36 39 0

38 0 36 0 39 39 35 37 0 0

39 0 37 0 40 40 36 38 0 0

40 0 38 0 41 41 37 39 0 0

41 0 39 0 42 42 38 40 0 0

42 0 0 0 0 0 0 0 0 0

END

# Data 4

# tn = transpose(n) -> n = transpose(tn)

tn[,1] tn[,2] tn[,3] tn[,4] tn[,5] tn[,6] tn[,7] tn[,8] tn[,9] tn[,10]

429 1 10 267 404 157 51 6 19 136

541 371 49 438 110 207 154 69 231 382

555 500 60 371 508 268 399 113 430 218

557 528 23 436 520 194 507 105 382 199

552 40 19 420 513 445 486 259 435 463

575 493 430 434 523 291 453 374 422 334

527 464 336 399 544 422 527 377 410 430

539 405 432 399 517 402 532 459 426 364

519 498 410 415 525 346 539 326 429 468

537 55 410 402 495 452 535 429 283 472

524 230 442 402 536 425 493 448 383 454

489 435 332 299 513 420 476 291 386 436

507 329 471 140 506 272 480 419 412 424

503 0 436 363 487 218 346 451 432 412

544 0 225 285 471 305 486 449 468 374

500 0 471 257 394 268 457 439 442 345

453 0 443 327 417 333 488 457 348 296

437 0 404 338 332 183 332 483 131 0

394 0 342 0 388 372 352 395 0 0

341 0 189 0 294 294 410 369 0 0

311 0 272 0 196 283 53 302 0 0

414 0 242 0 162 90 4 48 0 0

233 0 0 0 0 0 0 0 0 0

END