Improved quasi-likelihood estimation

# Improved quasi-likelihood estimation
### <div class="line-block">Ioannis Kosmidis<br><svg style="height:0.8em;top:.04em;position:relative;fill:#969696;" viewBox="0 0 496 512"><path d="M336.5 160C322 70.7 287.8 8 248 8s-74 62.7-88.5 152h177zM152 256c0 22.2 1.2 43.5 3.3 64h185.3c2.1-20.5 3.3-41.8 3.3-64s-1.2-43.5-3.3-64H155.3c-2.1 20.5-3.3 41.8-3.3 64zm324.7-96c-28.6-67.9-86.5-120.4-158-141.6 24.4 33.8 41.2 84.7 50 141.6h108zM177.2 18.4C105.8 39.6 47.8 92.1 19.3 160h108c8.7-56.9 25.5-107.8 49.9-141.6zM487.4 192H372.7c2.1 21 3.3 42.5 3.3 64s-1.2 43-3.3 64h114.6c5.5-20.5 8.6-41.8 8.6-64s-3.1-43.5-8.5-64zM120 256c0-21.5 1.2-43 3.3-64H8.6C3.2 212.5 0 233.8 0 256s3.2 43.5 8.6 64h114.6c-2-21-3.2-42.5-3.2-64zm39.5 96c14.5 89.3 48.7 152 88.5 152s74-62.7 88.5-152h-177zm159.3 141.6c71.4-21.2 129.4-73.7 158-141.6h-108c-8.8 56.9-25.6 107.8-50 141.6zM19.3 352c28.6 67.9 86.5 120.4 158 141.6-24.4-33.8-41.2-84.7-50-141.6h-108z"/></svg> <a href="http://ikosmidis.com">ikosmidis.com</a>  <svg style="height:0.8em;top:.04em;position:relative;fill:#969696;" viewBox="0 0 496 512"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg> <a href="https://github.com/ikosmidis">ikosmidis</a>  <svg style="height:0.8em;top:.04em;position:relative;fill:#969696;" viewBox="0 0 512 512"><path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"/></svg> <a href="https://twitter.com/ikosmidis_">ikosmidis_</a><br />
Reader in Data Science<br>University of Warwick & The Alan Turing Institute</div>
### <div class="line-block">e-Rum 2020<br>18 June 2020<br>Milan, Italy (in spirit)</div>

---

## Slide-deck version history

| Version | Last modified | Description                                            |
|---------|---------------|--------------------------------------------------------|
| v1.0    | 2020/06/19    | As presented at [eRum2020](https://2020.erum.io) on 2020/06/19 |
| v1.1    | 2020/06/20    | Fixed minor typos I noticed while presenting 😳        |

---

# Quasi likelihoods

---

## Incidence of leaf blotch on barley

.pull-left[
Percentage leaf area affected by *Rhynchosporium secalis* on 10 varieties of barley tested at 9 sites

```r
data("barley", package = "gnm")
str(barley)
## 'data.frame':	90 obs. of  3 variables:
##  $ y      : num  0.0005 0 0 0.001 0.0025 0.0005 0.005 0.013 0.015 0.015 ...
##  $ site   : Factor w/ 9 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ variety: Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
```
]

.of-note[Wedderburn, R W M (1974). Quasilikelihood functions, generalized linear 
models and the Gauss-Newton method. *Biometrika*, **61**<br><svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 496 512"><path d="M336.5 160C322 70.7 287.8 8 248 8s-74 62.7-88.5 152h177zM152 256c0 22.2 1.2 43.5 3.3 64h185.3c2.1-20.5 3.3-41.8 3.3-64s-1.2-43.5-3.3-64H155.3c-2.1 20.5-3.3 41.8-3.3 64zm324.7-96c-28.6-67.9-86.5-120.4-158-141.6 24.4 33.8 41.2 84.7 50 141.6h108zM177.2 18.4C105.8 39.6 47.8 92.1 19.3 160h108c8.7-56.9 25.5-107.8 49.9-141.6zM487.4 192H372.7c2.1 21 3.3 42.5 3.3 64s-1.2 43-3.3 64h114.6c5.5-20.5 8.6-41.8 8.6-64s-3.1-43.5-8.5-64zM120 256c0-21.5 1.2-43 3.3-64H8.6C3.2 212.5 0 233.8 0 256s3.2 43.5 8.6 64h114.6c-2-21-3.2-42.5-3.2-64zm39.5 96c14.5 89.3 48.7 152 88.5 152s74-62.7 88.5-152h-177zm159.3 141.6c71.4-21.2 129.4-73.7 158-141.6h-108c-8.8 56.9-25.6 107.8-50 141.6zM19.3 352c28.6 67.9 86.5 120.4 158 141.6-24.4-33.8-41.2-84.7-50-141.6h-108z"/></svg> [DOI: 10.2307/2334725](https://www.jstor.org/stable/2334725)]
]

---

---

## Beta regression: modelling percentages

.pull-left[
#### Beta density 
$$
f(y\,|\, \mu, \phi) = \frac{y^{\mu\phi - 1}(1-y)^{(1 - \mu)\phi - 1}}{B(\mu\phi, (1 - \mu)\phi)}
$$

| | |
| --- | --- |
| Mean | `$\mu$` |
| Variance |  `$\mu(1-\mu)/(1 + \phi)$` |

so `$\phi$` is a *precision* parameter

#### Mean and precision regression structures
$$
  g_1(\mu) = x^T\beta \quad \text{and} \quad g_2(\phi) = z^T\gamma
$$

]

.pull-right[
<img src="brquasi_eRum2020_kosmidis_files/figure-html/beta_density-1.png" width="432" />
]

---

## Estimation using <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 640 512"><path d="M425.7 256c-16.9 0-32.8-9-41.4-23.4L320 126l-64.2 106.6c-8.7 14.5-24.6 23.5-41.5 23.5-4.5 0-9-.6-13.3-1.9L64 215v178c0 14.7 10 27.5 24.2 31l216.2 54.1c10.2 2.5 20.9 2.5 31 0L551.8 424c14.2-3.6 24.2-16.4 24.2-31V215l-137 39.1c-4.3 1.3-8.8 1.9-13.3 1.9zm212.6-112.2L586.8 41c-3.1-6.2-9.8-9.8-16.7-8.9L320 64l91.7 152.1c3.8 6.3 11.4 9.3 18.5 7.3l197.9-56.5c9.9-2.9 14.7-13.9 10.2-23.1zM53.2 41L1.7 143.8c-4.6 9.2.3 20.2 10.1 23l197.9 56.5c7.1 2 14.7-1 18.5-7.3L320 64 69.8 32.1c-6.9-.8-13.5 2.7-16.6 8.9z"/></svg> `betareg`

Additive effects of site and variety for both mean and .yellow-background[for precision], 
.small[

```r
leaf_blotch_formula <- y ~ site + variety | `site + variety`
```
]

Maximum likelihood using default logit and log links for mean and precision
.small[

```r
library("betareg")
betareg(y ~ site + variety | site + variety, data = barley)
## Error in betareg(y ~ site + variety | site + variety, data = barley): invalid dependent variable, all observations must be in (0, 1)
```
]

---

## Quasi likelihoods

#### Mean regression and variance specification

`$\displaystyle g_1(\mu) = \eta = x^T\beta$` &nbsp; `%>%` &nbsp; `$\displaystyle {\rm var}(Y) = \phi v(\mu)$` &nbsp; (`%>%` &nbsp; Independence)

#### Quasi-likelihood equations
Then, estimates of `$\beta$` with optimal frequentist performance result by solving 
`$$\displaystyle \sum\frac{d\mu / d\eta}{\phi v(\mu)}(Y - \mu)x = 0$$`

---

## Quasi likelihoods in R

```r
m1 <- glm(y ~ site + variety,
   data = barley,
   family = quasibinomial())
```
]
]

```r
m2 <- glm(y ~ site + variety,
   data = barley,
   family = gnm::wedderburn())
```
]
]

---

# Finite sample bias and improvements

---

## Quality of the estimates: bias

Bootstrap with case resampling of size 10K
.small[

```r
set.seed(20200619)
boot_time <- system.time({
    boot_m1 <- car::Boot(m1, f = coef, R = 10000, method = "case")
})
```
]

---

## Quality of the estimates: bias

```r
boot_m1_summary <- summary(boot_m1)
print(boot_m1_summary, digits = 3)
## 
## Number of bootstrap replications R = 9993 
##             original bootBias bootSE bootMed
## (Intercept)   -8.055 -0.00515  0.533  -8.067
## siteB          1.639 -0.09291  0.510   1.616
## siteC          3.327  0.00782  0.509   3.329
## siteD          3.582 -0.07075  0.651   3.572
## siteE          3.583  0.03211  0.363   3.607
## siteF          3.893 -0.02130  0.561   3.909
## siteG          4.730  0.03165  0.517   4.752
## siteH          5.523  0.05129  0.432   5.556
## siteI          6.795  0.13748  0.482   6.907
## variety2       0.150 -0.14739  0.575   0.105
## variety3       0.689 -0.13242  0.570   0.638
## variety4       1.048 -0.09889  0.642   1.000
## variety5       1.615 -0.12684  0.715   1.490
## variety6       2.371 -0.17564  0.782   2.323
## variety7       2.571 -0.07867  0.618   2.503
## variety8       3.342 -0.02382  0.622   3.330
## variety9       3.500 -0.02987  0.544   3.497
## varietyX       4.253 -0.00472  0.530   4.248
```
]

---

## Reducing bias

Reduction of estimation bias has received great attention
since the early statistical literature

see, e.g.,

.of-note[Kosmidis I (2014). Bias in parametric estimation: reduction and
useful side-effects. *WIRE Computational Statistics*, **6**,
185-196<br><svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 496 512"><path d="M336.5 160C322 70.7 287.8 8 248 8s-74 62.7-88.5 152h177zM152 256c0 22.2 1.2 43.5 3.3 64h185.3c2.1-20.5 3.3-41.8 3.3-64s-1.2-43.5-3.3-64H155.3c-2.1 20.5-3.3 41.8-3.3 64zm324.7-96c-28.6-67.9-86.5-120.4-158-141.6 24.4 33.8 41.2 84.7 50 141.6h108zM177.2 18.4C105.8 39.6 47.8 92.1 19.3 160h108c8.7-56.9 25.5-107.8 49.9-141.6zM487.4 192H372.7c2.1 21 3.3 42.5 3.3 64s-1.2 43-3.3 64h114.6c5.5-20.5 8.6-41.8 8.6-64s-3.1-43.5-8.5-64zM120 256c0-21.5 1.2-43 3.3-64H8.6C3.2 212.5 0 233.8 0 256s3.2 43.5 8.6 64h114.6c-2-21-3.2-42.5-3.2-64zm39.5 96c14.5 89.3 48.7 152 88.5 152s74-62.7 88.5-152h-177zm159.3 141.6c71.4-21.2 129.4-73.7 158-141.6h-108c-8.8 56.9-25.6 107.8-50 141.6zM19.3 352c28.6 67.9 86.5 120.4 158 141.6-24.4-33.8-41.2-84.7-50-141.6h-108z"/></svg> [DOI: 10.1002/wics.1296](https://doi.org/10.1002/wics.1296)]

---

## R packages for bias reduction

Ample availability of R packages for specific model classes:

...

---

## No software solutions for `quasi` families
.small[

```r
library("brglm2")
update(m1, method = "brglm_fit")
## Error in brglm_fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, : `brglmFit` does not currently support the `quasi`, `quasipoisson` and `quasibinomial` families.
```
]

Standard bias reduction methods depend on

+ derivatives of the log-likelihood

+ expectation of products of those or simulation under the model

---

## What about the bootstrap?

Non-parametric bootstrap can correct finite-sample bias in quasi likelihood estimation. For reduced-bias estimates, simply do
.small[

```r
br_estimates <- with(boot_m1_summary, original - bootBias)
setNames(br_estimates, rownames(boot_m1_summary))
## (Intercept)       siteB       siteC       siteD       siteE       siteF 
##  -8.0494950   1.7319736   3.3186987   3.6529847   3.5509496   3.9145550 
##       siteG       siteH       siteI    variety2    variety3    variety4 
##   4.6983687   5.4714003   6.6571040   0.2974726   0.8218866   1.1470483 
##    variety5    variety6    variety7    variety8    variety9    varietyX 
##   1.7415489   2.5467951   2.6491935   3.3658578   3.5299168   4.2577307
```
]

#### But

Bootstrapping takes a while
.small[

```r
boot_time
##    user  system elapsed 
##  21.765   0.536  22.644
```
]

---

requires some under-the-hood conventions if the model involves factors
.small[

```r
boot_m1_summary$R
##  [1] 9993 9993 9993 9993 9993 9993 9993 9993 9993 9993 9993 9993 9993 9993 9993
## [16] 9993 9993 9993
# Not 10000???
```
]

or the choice of a residual for "more semi-parametric" bootstrap
.small[

```r
car::Boot(m1, f = coef, R = 10000, method = "residual")
## Error in Boot.glm(m1, f = coef, R = 10000, method = "residual"): Residual bootstrap is not implemented in the 'car' function 'Boot'.
##   Use the 'boot' function in the 'boot' package to write
##   your own version of residual bootstrap for a glm.
```
]

requires a nested bootstrap for computing standard errors for the improved estimates

---

Is there any way to reduce bias that

1. Depends only on the quasi likelihood equations?

2. Does not require the full specification of a data generating process?

3. Does not require resampling and/or ad-hoc conventions?

4. Has the same computational complexity as quasi likelihood estimation?

---

# Yes x 4

---

## Reduced-bias M-estimation: A general method

Under conditions, a reduced-bias version of the `$M$`-estimator `$\hat\beta \leftarrow \sum \psi(\beta) = 0$` is

`$$\tilde\beta \leftarrow \sum \psi(\beta) + D(\beta) = 0$$`

where

`$$D_t = - {\rm trace}\left\{j(\beta)^{-1} d_t(\beta) \right\} - \frac{1}{2} {\rm trace}\left[ j(\beta)^{-1} e(\beta) \left\{j(\beta)^{-1}\right\}^\top u_t(\beta) \right]$$`

and `$j(\beta)$`, `$u_t(\beta)$`, `$e_t(\beta)$`, `$d_t(\beta)$` are sums of products of the first two derivatives of `$\psi(\beta)$`.

.of-note[
Kosmidis I, Lunardon N (2020). Empirical bias-reducing adjustments to estimating functions.<br> <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 496 512"><path d="M336.5 160C322 70.7 287.8 8 248 8s-74 62.7-88.5 152h177zM152 256c0 22.2 1.2 43.5 3.3 64h185.3c2.1-20.5 3.3-41.8 3.3-64s-1.2-43.5-3.3-64H155.3c-2.1 20.5-3.3 41.8-3.3 64zm324.7-96c-28.6-67.9-86.5-120.4-158-141.6 24.4 33.8 41.2 84.7 50 141.6h108zM177.2 18.4C105.8 39.6 47.8 92.1 19.3 160h108c8.7-56.9 25.5-107.8 49.9-141.6zM487.4 192H372.7c2.1 21 3.3 42.5 3.3 64s-1.2 43-3.3 64h114.6c5.5-20.5 8.6-41.8 8.6-64s-3.1-43.5-8.5-64zM120 256c0-21.5 1.2-43 3.3-64H8.6C3.2 212.5 0 233.8 0 256s3.2 43.5 8.6 64h114.6c-2-21-3.2-42.5-3.2-64zm39.5 96c14.5 89.3 48.7 152 88.5 152s74-62.7 88.5-152h-177zm159.3 141.6c71.4-21.2 129.4-73.7 158-141.6h-108c-8.8 56.9-25.6 107.8-50 141.6zM19.3 352c28.6 67.9 86.5 120.4 158 141.6-24.4-33.8-41.2-84.7-50-141.6h-108z"/></svg> [arXiv:2001.03786](https://arxiv.org/abs/2001.03786)]

---

## <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 640 512"><path d="M425.7 256c-16.9 0-32.8-9-41.4-23.4L320 126l-64.2 106.6c-8.7 14.5-24.6 23.5-41.5 23.5-4.5 0-9-.6-13.3-1.9L64 215v178c0 14.7 10 27.5 24.2 31l216.2 54.1c10.2 2.5 20.9 2.5 31 0L551.8 424c14.2-3.6 24.2-16.4 24.2-31V215l-137 39.1c-4.3 1.3-8.8 1.9-13.3 1.9zm212.6-112.2L586.8 41c-3.1-6.2-9.8-9.8-16.7-8.9L320 64l91.7 152.1c3.8 6.3 11.4 9.3 18.5 7.3l197.9-56.5c9.9-2.9 14.7-13.9 10.2-23.1zM53.2 41L1.7 143.8c-4.6 9.2.3 20.2 10.1 23l197.9 56.5c7.1 2 14.7-1 18.5-7.3L320 64 69.8 32.1c-6.9-.8-13.5 2.7-16.6 8.9z"/></svg> `brquasi`

Implements a quasi Newton iteration to solve the adjusted quasi-likelihood equations

`glm` as an interface with `method = "brquasi_fit"` (`glm`'s default is `method = "glm.fit"`)

```r
library("brquasi")
system.time({
    m1_br <- update(m1, method = "brquasi_fit")
})
##    user  system elapsed 
##   0.447   0.009   0.471

## or even faster substract an estimate of the bias
system.time({
    m1_br <- update(m1, method = "brquasi_fit", type = "eRBM")
})
##    user  system elapsed 
##   0.011   0.000   0.011
```
]

see `?brquasi_control` for finer control

---

```r
class(m1_br)
## [1] "brquasiFit" "glm"        "lm"
```

Base methods for `glm` objects work out-of-the-box (at least they should!) and most are supported by theory

Direct integration into packages with `glm` methods

---

```r
broom::glance(m1)
## # A tibble: 1 x 7
##   null.deviance df.null logLik   AIC   BIC deviance df.residual
##           <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
## 1          40.8      89     NA    NA    NA     6.13          72
broom::glance(m1_br)
## # A tibble: 1 x 7
##   null.deviance df.null logLik   AIC   BIC deviance df.residual
##           <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
## 1          40.8      89     NA    NA    NA     6.15          72
```
]

---

```r
broom::tidy(m1_br)
## # A tibble: 18 x 5
##    term        estimate std.error statistic  p.value
##    <chr>          <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)   -8.04      1.41     -5.70  1.18e- 8
##  2 siteB          1.70      1.42      1.20  2.32e- 1
##  3 siteC          3.33      1.33      2.49  1.27e- 2
##  4 siteD          3.63      1.33      2.73  6.34e- 3
##  5 siteE          3.57      1.33      2.68  7.28e- 3
##  6 siteF          3.92      1.32      2.96  3.09e- 3
##  7 siteG          4.71      1.32      3.57  3.60e- 4
##  8 siteH          5.49      1.32      4.16  3.19e- 5
##  9 siteI          6.70      1.33      5.05  4.37e- 7
## 10 variety2       0.229     0.723     0.316 7.52e- 1
## 11 variety3       0.782     0.672     1.16  2.44e- 1
## 12 variety4       1.11      0.652     1.70  8.91e- 2
## 13 variety5       1.68      0.628     2.68  7.43e- 3
## 14 variety6       2.46      0.611     4.03  5.66e- 5
## 15 variety7       2.61      0.610     4.28  1.88e- 5
## 16 variety8       3.34      0.605     5.52  3.36e- 8
## 17 variety9       3.50      0.605     5.79  7.13e- 9
## 18 varietyX       4.24      0.608     6.98  2.99e-12
```
]

---

```r
broom::augment(m1_br)
## # A tibble: 90 x 10
##         y site  variety .fitted .se.fit   .resid    .hat .sigma .cooksd
##     <dbl> <fct> <fct>     <dbl>   <dbl>    <dbl>   <dbl>  <dbl>   <dbl>
##  1 0.0005 A     1         -8.04    1.41  9.18e-3 0.00721  0.294 4.51e-7
##  2 0      A     2         -7.81    1.39 -2.84e-2 0.00885  0.294 2.28e-6
##  3 0      A     3         -7.26    1.37 -3.75e-2 0.0148   0.294 6.71e-6
##  4 0.001  A     4         -6.93    1.36  8.29e-4 0.0201   0.294 9.09e-9
##  5 0.0025 A     5         -6.36    1.34  1.74e-2 0.0349   0.294 8.14e-6
##  6 0.0005 A     6         -5.58    1.33 -6.72e-2 0.0744   0.294 1.54e-4
##  7 0.005  A     7         -5.43    1.33  9.56e-3 0.0858   0.294 6.16e-6
##  8 0.013  A     8         -4.70    1.32  3.96e-2 0.174    0.294 2.85e-4
##  9 0.015  A     9         -4.54    1.32  4.09e-2 0.203    0.294 3.79e-4
## 10 0.015  A     X         -3.80    1.31 -4.98e-2 0.410    0.294 1.63e-3
## # … with 80 more rows, and 1 more variable: .std.resid <dbl>
```
]

---

```r
lmtest::waldtest(update(m1_br, . ~ . -site), m1_br)
## Wald test
## 
## Model 1: y ~ variety
## Model 2: y ~ site + variety
##   Res.Df Df      F    Pr(>F)    
## 1     80                        
## 2     72  8 18.529 9.609e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
]

---

## <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 512 512"><path d="M512 144v288c0 26.5-21.5 48-48 48H48c-26.5 0-48-21.5-48-48V144c0-26.5 21.5-48 48-48h88l12.3-32.9c7-18.7 24.9-31.1 44.9-31.1h125.5c20 0 37.9 12.4 44.9 31.1L376 96h88c26.5 0 48 21.5 48 48zM376 288c0-66.2-53.8-120-120-120s-120 53.8-120 120 53.8 120 120 120 120-53.8 120-120zm-32 0c0 48.5-39.5 88-88 88s-88-39.5-88-88 39.5-88 88-88 88 39.5 88 88z"/></svg>

[![Travis-CI Build
Status](https://travis-ci.org/ikosmidis/brquasi.svg?branch=master)](https://travis-ci.org/ikosmidis/brquasi)
[![CRAN\_Status\_Badge](https://www.r-pkg.org/badges/version/brquasi)](https://cran.r-project.org/package=brquasi)
[![Licence](https://img.shields.io/badge/licence-GPL--3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.en.html)

<br/>

```r
remotes::install_github("ikosmidis/brquasi")
```
]
]

.small[
`$\mathop{\rm solve} \left\{ \sum \psi + D = 0\right\}$`, `$D_t = -{\rm trace}\left\{j^{-1} d_t \right\} - \frac{1}{2} {\rm trace}\left[ j^{-1} e \left\{j^{-1}\right\}^\top u_t \right]$`
]

```r
library("brquasi")
glm(..., method = "brquasi_fit")
```
]
]

<br/>