#
# Ship data from McCullagh and Nelder (1989), Section 6.3.
# With SAS GLMMIX, one obtains:
#/*---poisson regression (log-linear model) with random effects,
#     parameterization for TYPE matches McCullagh and Nelder's---*/
#glimmix(data=ship,
#  procopt=order=data,
#  stmts=%str(
#     class type year period;
#     model y = type / solution;
#     random year|period;
#     estimate 'E vs. Others' type -1 -1 -1 4 -1 / divisor=4 cl;
#  ),
#  error=poisson,
#  link=log,
#  offset=service
#  Covariance Parameter
#            Estimates
#     Cov Parm        Estimate
#     year              0.1174
#     period           0.07066
#     year*period            0
#     Residual          1.6702
#
#    Solution for Fixed Effects
#   Effect       type    Estimate       Error      DF    t Value    Pr > |t|
#   Intercept             -5.6799      0.3286       1     -17.28      0.0368
#   type         B        -0.5798      0.2277      23      -2.55      0.0180
#   type         C        -0.6984      0.4248      23      -1.64      0.1138
#   type         D       -0.08703      0.3746      23      -0.23      0.8183
#   type         E         0.3301      0.3046      23       1.08      0.2897
#   type         A              0           .       .        .         .
#
# lmer:
#
# Formula: incidents ~ type + (1 | period/year) + (1 | year) 
#    Data: ships[ships$service > 0, ]             lmer (Ver 9995.2)
#  Family: poisson(log link)                      Family: poisson(log link)
#       AIC     BIC    logLik deviance                 AIC      BIC    logLik deviance
#  72.01782 84.2287 -28.00891 56.01782            72.88705 83.57157 -29.44352 58.88705
# Random effects:                                Random effects:
#  Groups      Name        Variance Std.Dev.      Groups      Name        Variance Std.Dev.
#  period/year (Intercept) 0.04927  0.22197       period/year (Intercept) 0.028332 0.16832
#  year        (Intercept) 0.11173  0.33426       year        (Intercept) 0.064246 0.25347
#  Residual                1.73905  1.31873      
# # of obs: 34, groups: period/year, 6; year, 4
# Estimated scale (compare to 1)  1.217912 
# Fixed effects:
#              Estimate Std. Error  z value  Pr(>|z|)    
# (Intercept) -5.573280   0.217482 -25.6264 < 2.2e-16 ***
# typeB       -0.605271   0.176118  -3.4367 0.0005888 ***
# typeC       -0.716560   0.328634  -2.1804 0.0292262 *  
# typeD       -0.061478   0.289935  -0.2120 0.8320749    
# typeE        0.309430   0.235577   1.3135 0.1890160    
#
# Generalized linear mixed model fit using Laplace
# Formula: incidents ~ type + (1 | period/year) + (1 | year)
#    Data: ships[ships$service > 0, ]
#  Family: poisson(log link)
#      AIC      BIC   logLik deviance
#  72.1928 82.87732 -29.0964  58.1928
# Random effects:
#  Groups      Name        Variance Std.Dev.
#  period/year (Intercept) 0.053025 0.23027
#  year        (Intercept) 0.107726 0.32822
# # of obs: 34, groups: period/year, 6; year, 4
# 
# Estimated scale (compare to 1)  1.217912
# 
# Fixed effects:
#              Estimate Std. Error  z value  Pr(>|z|)
# (Intercept) -5.582762   0.250875 -22.2531 < 2.2e-16 ***
# typeB       -0.585531   0.176731  -3.3131 0.0009226 ***
# typeC       -0.708410   0.328865  -2.1541 0.0312316 *
# typeD       -0.063597   0.290330  -0.2190 0.8266119
# typeE        0.312139   0.235757   1.3240 0.1855066
# 
# And with the alpha version of lmer2
# > lmer(formula = incidents ~ type + (1 | period) + (1 | year),
#  +  data = ships[ships$service>0,], family = poisson, offset = log(service))
#  Generalized linear mixed model fit by the Laplace approximation
#  Formula: incidents ~ type + (1 | period) + (1 | year)
#     Data: ships[ships$service > 0, ]
#     AIC   BIC logLik deviance
#   69.13 79.81 -27.56    55.13
#   Random effects:
#    Groups Name        Variance Std.Dev.
#    period (Intercept) 0.046202 0.21495
#    year   (Intercept) 0.102408 0.32001
#   Number of obs: 34, groups: period, 2; year, 4
#   Fixed effects:
#               Estimate Std. Error z value Pr(>|z|)
#   (Intercept) -5.69649    0.27339 -20.836  < 2e-16 ***
#   typeB       -0.57005    0.17659  -3.228  0.00125 **
#   typeC       -0.69512    0.32883  -2.114  0.03453 *
#   typeD       -0.08378    0.29015  -0.289  0.77277
#   typeE        0.32827    0.23578   1.392  0.16383
#
set ple -1
set che off
set imp -1
set loc service qua
set loc incidents qua
set loc type mar
read pedigree inline
#year           period         service y   type
1965-69 1  m1975-79 f1975-79 x 9.9218  53 B/B 
1965-69 2  m1975-79 f1975-79 x 6.5162  1  C/C 
1965-69 3  m1975-79 f1975-79 x 5.2575  0  D/D 
1965-69 4  m1975-79 f1975-79 x 6.0799  7  E/E 
1965-69 5  m1975-79 f1975-79 x 6.9985  4  A/A 
1965-69 6  m1960-74 f1960-74 x 6.9985  3  A/A 
1965-69 7  m1960-74 f1960-74 x 10.2615 58 B/B 
1965-69 8  m1960-74 f1960-74 x 6.6606  0  C/C 
1965-69 9  m1960-74 f1960-74 x 5.6630  0  D/D 
1965-69 10 m1960-74 f1960-74 x 6.6708  7  E/E 
1970-64 11 m1960-74 f1960-74 x 7.3212  6  A/A 
1970-64 12 m1960-74 f1960-74 x 8.8628  12 B/B 
1970-64 13 m1960-74 f1960-74 x 6.6631  6  C/C 
1970-64 14 m1960-74 f1960-74 x 5.8551  2  D/D 
1970-64 15 m1960-74 f1960-74 x 7.0536  5  E/E 
1970-64 16 m1975-79 f1975-79 x 8.1176  18 A/A 
1970-64 17 m1975-79 f1975-79 x 9.4803  44 B/B 
1970-64 18 m1975-79 f1975-79 x 7.5746  2  C/C 
1970-64 19 m1975-79 f1975-79 x 7.0967  11 D/D 
1970-64 20 m1975-79 f1975-79 x 7.6783  12 E/E 
1975-69 21 m1975-79 f1975-79 x 7.7160  11 A/A 
1975-69 22 m1975-79 f1975-79 x 8.8702  18 B/B 
1975-69 23 m1975-79 f1975-79 x 5.6131  1  C/C 
1975-69 24 m1975-79 f1975-79 x 7.6261  4  D/D 
1975-69 25 m1975-79 f1975-79 x 6.2953  1  E/E 
1960-64 26 m1960-74 f1960-74 x 4.8442  0  A/A 
1960-64 27 m1960-74 f1960-74 x 10.7118 39 B/B 
1960-64 28 m1960-74 f1960-74 x 7.0724  1  C/C 
1960-64 29 m1960-74 f1960-74 x 5.5255  0  D/D 
1960-64 30 m1960-74 f1960-74 x 3.8067  0  E/E 
1960-64 31 m1975-79 f1975-79 x 4.1431  0  A/A 
1960-64 32 m1975-79 f1975-79 x 9.7513  29 B/B 
1960-64 33 m1975-79 f1975-79 x 6.3135  1  C/C 
1960-64 34 m1975-79 f1975-79 x 4.6540  0  D/D 
;;;;
run
macro mktyp
  set loc typ%1 aff
  if (istyp type) then typ%1 = (type=="%1/%1")
;;;;
mktyp { B C D E }
reg incidents = typB typC typD typE offset service poisson
# set chains 4
# set burn 1000
# set iter 2000
fpm incidents nqtl 0 c s cov typB + typC + typD + typE \
    lik poisson offset service
# fpm incidents nqtl 0 c cov typB + typC + typD + typE \
#     lik poisson offset service
# lrt
