#
# Data from CA McGilchrist, CW Aisbett (1991), Regression with frailty in
#  survival analysis.  _Biometrics_ *47*, 461-66.
#
# Recurrence times to catheter infection in 
# patients using portable
# dialysis equipment. Catheters may be removed for reasons other
# than infection, in which case the observation is censored.  Each
# patient has exactly 2 observations.
#
# BUGs example of Weibull frailty     Sib-pair
# regression:                         chains=4 burn=1000
# burnin=1000 iterations=10000        iterations=50000
# 
# par    mean   SD     MC-SE          
# alpha  -4.529 0.9036 0.0624         4.4571  0.8687  0.0638
# GN     0.1265 0.5679 0.0186        -0.0397  0.5219  0.0138
# AN     0.5995 0.5781 0.0220        -0.5586  0.5146  0.0118
# PKD    -1.198 0.8483 0.0315         1.2897  0.7660  0.0225
# Female -1.945 0.5019 0.0280         1.9924  0.4931  0.0254
# r      1.205  0.1711 0.0125         1.1903
# sigma  0.6367 0.3802 0.0316         0.6201  0.2489  0.0173
#
#
# Results from survreg(Surv(time, status) ~ age + sex + disease +
#                      frailty.gaussian(id), data = kidney)
#
#                      coef     se(coef) se2     Chisq DF   p      
# (Intercept)           2.11528 0.7612   0.36962  7.72  1.0 5.5e-03
# age                  -0.00301 0.0127   0.00595  0.06  1.0 8.1e-01
# sex                   1.56195 0.3518   0.17657 19.71  1.0 9.0e-06
# diseaseGN            -0.15705 0.4563   0.22798  0.12  1.0 7.3e-01
# diseaseAN            -0.55862 0.4565   0.22832  1.50  1.0 2.2e-01
# diseasePKD            0.63778 0.6346   0.32905  1.01  1.0 3.1e-01
# frailty.gaussian(id)                           58.62 22.5 4.6e-05
# 
# Scale= 0.55
# 
# Iterations: 6 outer, 28 Newton-Raphson
#      Variance of random effect= 0.628 
# Degrees of freedom for terms=  0.2  0.2  0.3  0.8 22.5  0.9 
# Likelihood ratio test=81.7  on 22.8 df, p=1.53e-08  n= 76 
# 
#
# And survreg(Surv(time, status)~ age + sex + disease + 
#             frailty.gaussian(id, method="aic", caic=T), kidney))
# 
#                Value Std. Error      z        p
# (Intercept)  2.06991     0.6667  3.105 1.90e-03
# age         -0.00277     0.0105 -0.264 7.92e-01
# sex          1.59980     0.3108  5.147 2.65e-07
# diseaseGN   -0.08882     0.3840 -0.231 8.17e-01
# diseaseAN   -0.53129     0.3788 -1.403 1.61e-01
# diseasePKD   0.96965     0.5749  1.686 9.17e-02
# Log(scale)  -0.31929     0.1263 -2.527 1.15e-02
# 
# Scale= 0.727 
# 
# Weibull distribution
# Loglik(model)= -314   Loglik(intercept only)= -340.9
#         Chisq= 53.81 on 12.7 degrees of freedom, p= 5.2e-07 
# Number of Newton-Raphson Iterations:  8 32 
# 
set loc time    qua . recurrence times to catheter infection
set loc status  qua .
set loc age     qua . in years
set loc disease qua . 1=Other 2=GN 2=AN 3=PKD
set loc frail   qua . frailty estimate from McGilchrist 1991
set loc GN      aff . 
set loc AN      aff .
set loc PKD     aff .
read pedigree inline 
# id obs   sex time status age disease frail GN AN PKD 
1  1   x x 1   8    y      28  1       2.3   n  n  n   
1  2   x x 1   16   y      28  1       2.3   n  n  n   
2  3   x x 2   23   y      48  2       1.9   y  n  n   
2  4   x x 2   13   n      48  2       1.9   y  n  n   
3  5   x x 1   22   y      32  1       1.2   n  n  n   
3  6   x x 1   28   y      32  1       1.2   n  n  n   
4  7   x x 2   447  y      31  1       0.5   n  n  n   
4  8   x x 2   318  y      32  1       0.5   n  n  n   
5  9   x x 1   30   y      10  1       1.5   n  n  n   
5  10  x x 1   12   y      10  1       1.5   n  n  n   
6  11  x x 2   24   y      16  1       1.1   n  n  n   
6  12  x x 2   245  y      17  1       1.1   n  n  n   
7  13  x x 1   7    y      51  2       3     y  n  n   
7  14  x x 1   9    y      51  2       3     y  n  n   
8  15  x x 2   511  y      55  2       0.5   y  n  n   
8  16  x x 2   30   y      56  2       0.5   y  n  n   
9  17  x x 2   53   y      69  3       0.7   n  y  n   
9  18  x x 2   196  y      69  3       0.7   n  y  n   
10 19  x x 1   15   y      51  2       0.4   y  n  n   
10 20  x x 1   154  y      52  2       0.4   y  n  n   
11 21  x x 2   7    y      44  3       0.6   n  y  n   
11 22  x x 2   333  y      44  3       0.6   n  y  n   
12 23  x x 2   141  y      34  1       1.2   n  n  n   
12 24  x x 2   8    n      34  1       1.2   n  n  n   
13 25  x x 2   96   y      35  3       1.4   n  y  n   
13 26  x x 2   38   y      35  3       1.4   n  y  n   
14 27  x x 2   149  n      42  3       0.4   n  y  n   
14 28  x x 2   70   n      42  3       0.4   n  y  n   
15 29  x x 2   536  y      17  1       0.4   n  n  n   
15 30  x x 2   25   n      17  1       0.4   n  n  n   
16 31  x x 1   17   y      60  3       1.1   n  y  n   
16 32  x x 1   4    n      60  3       1.1   n  y  n   
17 33  x x 2   185  y      60  1       0.8   n  n  n   
17 34  x x 2   177  y      60  1       0.8   n  n  n   
18 35  x x 2   292  y      43  1       0.8   n  n  n   
18 36  x x 2   114  y      44  1       0.8   n  n  n   
19 37  x x 2   22   n      53  2       0.5   y  n  n   
19 38  x x 2   159  n      53  2       0.5   y  n  n   
20 39  x x 2   15   y      44  1       1.3   n  n  n   
20 40  x x 2   108  n      44  1       1.3   n  n  n   
21 41  x x 1   152  y      46  4       0.2   n  n  y   
21 42  x x 1   562  y      47  4       0.2   n  n  y   
22 43  x x 2   402  y      30  1       0.6   n  n  n   
22 44  x x 2   24   n      30  1       0.6   n  n  n   
23 45  x x 2   13   y      62  3       1.7   n  y  n   
23 46  x x 2   66   y      63  3       1.7   n  y  n   
24 47  x x 2   39   y      42  3       1     n  y  n   
24 48  x x 2   46   n      43  3       1     n  y  n   
25 49  x x 1   12   y      43  3       0.7   n  y  n   
25 50  x x 1   40   y      43  3       0.7   n  y  n   
26 51  x x 2   113  n      57  3       0.5   n  y  n   
26 52  x x 2   201  y      58  3       0.5   n  y  n   
27 53  x x 2   132  y      10  2       1.1   y  n  n   
27 54  x x 2   156  y      10  2       1.1   y  n  n   
28 55  x x 2   34   y      52  3       1.8   n  y  n   
28 56  x x 2   30   y      52  3       1.8   n  y  n   
29 57  x x 1   2    y      53  2       1.5   y  n  n   
29 58  x x 1   25   y      53  2       1.5   y  n  n   
30 59  x x 2   130  y      54  2       1.5   y  n  n   
30 60  x x 2   26   y      54  2       1.5   y  n  n   
31 61  x x 2   27   y      56  3       1.7   n  y  n   
31 62  x x 2   58   y      56  3       1.7   n  y  n   
32 63  x x 2   5    n      50  3       1.3   n  y  n   
32 64  x x 2   43   y      51  3       1.3   n  y  n   
33 65  x x 2   152  y      57  4       2.9   n  n  y   
33 66  x x 2   30   y      57  4       2.9   n  n  y   
34 67  x x 2   190  y      44  2       0.7   y  n  n   
34 68  x x 2   5    n      45  2       0.7   y  n  n   
35 69  x x 2   119  y      22  1       2.2   n  n  n   
35 70  x x 2   8    y      22  1       2.2   n  n  n   
36 71  x x 2   54   n      42  1       0.7   n  n  n   
36 72  x x 2   16   n      42  1       0.7   n  n  n   
37 73  x x 2   6    n      52  4       2.1   n  n  y   
37 74  x x 2   78   y      52  4       2.1   n  n  y   
38 75  x x 1   63   y      60  4       1.2   n  n  y   
38 76  x x 1   8    n      60  4       1.2   n  n  y   
;;;;
run
set loc fem aff
fem=female
reg time = age fem GN AN PKD weibull status
set iter 5000
set chains 4
fpm time lik weibull cens status nqtl 0 c cov age + fem + GN + AN + PKD
