R Under development (unstable) (2024-09-10 r87115 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> options(digits = 4)
> 
> ## Package and data
> library("betareg")
> data("GasolineYield", package = "betareg")
> 
> ## Same results as in Table 1 in Kosmidis and Firth (2010, EJS,
> ## http://dx.doi.org/10.1214/10-EJS579)
> gyML <- betareg(yield ~ batch + temp, data = GasolineYield, link.phi = "identity")
> gyBC <- betareg(yield ~ batch + temp, data = GasolineYield, link.phi = "identity", type = "BC")
> gyBR <- betareg(yield ~ batch + temp, data = GasolineYield, link.phi = "identity", type = "BR")
> 
> ## Coefficients and standard errors
> se <- function(obj, ...) sqrt(diag(vcov(obj, ...)))
> sapply(list(gyML, gyBC, gyBR), coef)
                 [,1]      [,2]      [,3]
(Intercept)  -6.15957  -6.14837  -6.14171
batch1        1.72773   1.72484   1.72325
batch2        1.32260   1.32009   1.31860
batch3        1.57231   1.56928   1.56734
batch4        1.05971   1.05788   1.05677
batch5        1.13375   1.13165   1.13024
batch6        1.04016   1.03829   1.03714
batch7        0.54369   0.54309   0.54242
batch8        0.49590   0.49518   0.49446
batch9        0.38579   0.38502   0.38459
temp          0.01097   0.01094   0.01093
(phi)       440.27839 261.20610 261.03777
> sapply(list(gyML, gyBC, gyBR), se)
                 [,1]      [,2]      [,3]
(Intercept) 1.823e-01 2.360e-01 2.359e-01
batch1      1.012e-01 1.311e-01 1.311e-01
batch2      1.179e-01 1.526e-01 1.526e-01
batch3      1.161e-01 1.503e-01 1.503e-01
batch4      1.024e-01 1.325e-01 1.325e-01
batch5      1.035e-01 1.340e-01 1.340e-01
batch6      1.060e-01 1.373e-01 1.373e-01
batch7      1.091e-01 1.412e-01 1.412e-01
batch8      1.089e-01 1.410e-01 1.410e-01
batch9      1.186e-01 1.535e-01 1.535e-01
temp        4.126e-04 5.339e-04 5.338e-04
(phi)       1.100e+02 6.526e+01 6.522e+01
> 
> ## Same results as in Table 3 in Kosmidis and Firth (2010, EJS,
> ## http://dx.doi.org/10.1214/10-EJS579). BR and BC estimates in the
> ## latter study were calculated in a different way than betareg
> ## computes them which provides some indication on the correctness of
> ## implementation
> gyMLlog <- betareg(yield ~ batch + temp, data = GasolineYield, link.phi = "log")
> gyBClog <- betareg(yield ~ batch + temp, data = GasolineYield, link.phi = "log", type = "BC")
> gyBRlog <- betareg(yield ~ batch + temp, data = GasolineYield, link.phi = "log", type = "BR")
> sapply(list(gyMLlog, gyBClog, gyBRlog), coef)
                      [,1]     [,2]     [,3]
(Intercept)       -6.15957 -6.14837 -6.14259
batch1             1.72773  1.72484  1.72347
batch2             1.32260  1.32009  1.31880
batch3             1.57231  1.56928  1.56758
batch4             1.05971  1.05788  1.05691
batch5             1.13375  1.13165  1.13041
batch6             1.04016  1.03829  1.03729
batch7             0.54369  0.54309  0.54248
batch8             0.49590  0.49518  0.49453
batch9             0.38579  0.38502  0.38465
temp               0.01097  0.01094  0.01093
(phi)_(Intercept)  6.08741  5.71191  5.61608
> sapply(list(gyMLlog, gyBClog, gyBRlog), se)
                       [,1]      [,2]      [,3]
(Intercept)       0.1823247 0.2194387 0.2299769
batch1            0.1012294 0.1218919 0.1277689
batch2            0.1179020 0.1419260 0.1487521
batch3            0.1161045 0.1397813 0.1465144
batch4            0.1023598 0.1232343 0.1291683
batch5            0.1035232 0.1246546 0.1306668
batch6            0.1060365 0.1276729 0.1338262
batch7            0.1091275 0.1313272 0.1376330
batch8            0.1089257 0.1311226 0.1374341
batch9            0.1185933 0.1427822 0.1496578
temp              0.0004126 0.0004966 0.0005204
(phi)_(Intercept) 0.2499001 0.2498567 0.2498428
> 
> ## Fit with temp as a dispersion covariate
> gy2ML <- betareg(yield ~ batch + temp | temp, data = GasolineYield)
> gy2BC <- betareg(yield ~ batch + temp | temp, data = GasolineYield, type = "BC")
> gy2BR <- betareg(yield ~ batch + temp | temp, data = GasolineYield, type = "BR")
> sapply(list(gy2ML, gy2BC, gy2BR), coef)
                      [,1]     [,2]      [,3]
(Intercept)       -5.92324 -5.91682 -6.085354
batch1             1.60199  1.60063  1.680352
batch2             1.29727  1.29591  1.312044
batch3             1.56534  1.56362  1.570254
batch4             1.03007  1.02919  1.048684
batch5             1.15416  1.15318  1.137410
batch6             1.01944  1.01857  1.031739
batch7             0.62226  0.62171  0.571073
batch8             0.56458  0.56416  0.518957
batch9             0.35944  0.35907  0.375556
temp               0.01036  0.01035  0.010782
(phi)_(Intercept)  1.36409  1.98198  4.366840
(phi)_temp         0.01457  0.01148  0.003763
> sapply(list(gy2ML, gy2BC, gy2BR), se)
                       [,1]      [,2]      [,3]
(Intercept)       0.1835262 0.2216756 0.2340856
batch1            0.0638561 0.0860389 0.1165475
batch2            0.0991001 0.1239755 0.1439965
batch3            0.0997392 0.1238178 0.1420904
batch4            0.0632882 0.0855579 0.1171486
batch5            0.0656427 0.0877328 0.1185384
batch6            0.0663510 0.0892694 0.1214629
batch7            0.0656325 0.0892209 0.1242404
batch8            0.0601846 0.0832203 0.1211889
batch9            0.0671406 0.0924409 0.1324988
temp              0.0004362 0.0005247 0.0005415
(phi)_(Intercept) 1.2257812 1.2266705 1.2323159
(phi)_temp        0.0036183 0.0036200 0.0036336
> 
> 
> proc.time()
   user  system elapsed 
   0.92    0.14    1.04