Tutorial 10.5a - Logistic regression and proportional and percentage data

We perform the logistic regression using the glm() function. The most important (=commonly used) parameters/arguments for logistic regression are:

• formula: the linear model relating the response variable to the linear predictor
• family: specification of the error distribution (and link function). Can be specified as either a string or as one of the family functions (which allows for the explicit declaration of the link function). For examples:
  • family="binomial" or equivalently family=binomial(link="logit")
  • family=binomial(link="probit")
• data: the data frame containing the data
• weights: vector of weights to be used in the fitting process. These are particularly useful when dealing with grouped binary data from unequal sample sizes (see below).

I will demonstrate logistic regression with a range of possible link functions (each of which yield different parameter interpretations):

• logit $$log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmL <- glm(y ~ x, data = dat, family = "binomial")
  

• probit $$\phi^{-1}\left(\pi\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmP <- glm(y ~ x, data = dat, family = binomial(link = "probit"))
  

• complimentary log-log $$log(-log(1-\pi))=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmC <- glm(y ~ x, data = dat, family = binomial(link = "cloglog"))
  

We perform the logistic regression using the glmmTMB() function. The most important (=commonly used) parameters/arguments for logistic regression are:

• formula: the linear model relating the response variable to the linear predictor
• family: specification of the error distribution (and link function). Can be specified as either a string or as one of the family functions (which allows for the explicit declaration of the link function). For examples:
  • family="binomial" or equivalently family=binomial(link="logit")
  • family=binomial(link="probit")
• data: the data frame containing the data
• weights: vector of weights to be used in the fitting process. These are particularly useful when dealing with grouped binary data from unequal sample sizes (see below).

I will demonstrate logistic regression with a range of possible link functions (each of which yield different parameter interpretations):

• logit $$log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  library(glmmTMB)
  dat.glmmTMBL <- glmmTMB(y ~ x, data = dat, family = "binomial")
  

• probit $$\phi^{-1}\left(\pi\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmmTMBP <- glmmTMB(y ~ x, data = dat, family = binomial(link = "probit"))
  

• complimentary log-log $$log(-log(1-\pi))=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmmTMBC <- glmmTMB(y ~ x, data = dat, family = binomial(link = "cloglog"))
  

We perform the logistic regression using the gamlss() function. The most important (=commonly used) parameters/arguments for logistic regression are:

• formula: the linear model relating the response variable to the linear predictor
• family: specification of the error distribution (and link function). This is a different specification from glm() and supports a much wider variety of families, including mixtures. For examples:
  • family=BI() or equivalently family=BI(mu.link="logit")
  • family=BI(mu.link="probit")
• data: the data frame containing the data
• weights: vector of weights to be used in the fitting process. These are particularly useful when dealing with grouped binary data from unequal sample sizes (see below).

I will demonstrate logistic regression with a range of possible link functions (each of which yield different parameter interpretations):

• logit $$log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  library(gamlss)
  dat.gamlssL <- gamlss(y ~ x, data = dat, family = BI())
  

  GAMLSS-RS iteration 1: Global Deviance = 11.6514 
  GAMLSS-RS iteration 2: Global Deviance = 11.6514 
  

• probit $$\phi^{-1}\left(\pi\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.gamlssP <- gamlss(y ~ x, data = dat, family = BI(mu.link = "probit"))
  

  GAMLSS-RS iteration 1: Global Deviance = 11.4522 
  GAMLSS-RS iteration 2: Global Deviance = 11.4522 
  

• complimentary log-log $$log(-log(1-\pi))=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.gamlssC <- gamlss(y ~ x, data = dat, family = BI(mu.link = "cloglog"))
  

  GAMLSS-RS iteration 1: Global Deviance = 12.1082 
  GAMLSS-RS iteration 2: Global Deviance = 12.1082 
  

Lets explore the diagnostics - particularly the residuals

plot(dat.glmL)

Unfortunately, unlike with linear models (Gaussian family), the expected distribution of data (residuals) varies over the range of fitted values for numerous (often competing) ways that make diagnosing (and attributing causes thereof) miss-specified generalized linear models from standard residual plots very difficult. The use of standardized (Pearson) residuals or deviance residuals can partly address this issue, yet they still do not offer completely consistent diagnoses across all issues (miss-specified model, over-dispersion, zero-inflation).

An alternative approach is to use simulated data from the fitted model to calculate an empirical cumulative density function from which residuals are are generated as values corresponding to the observed data along the density function. The (rather Bayesian) rationale is that if the model is correctly specified, then the observed data can be considered as a random draw from the fitted model. If this is the case, then residuals calculated from empirical cumulative density functions based on data simulated from the fitted model, should be totally uniform (flat) across the range of the linear predictor (fitted values) regardless of the model (Binomial, Poisson, linear, quadratic, hierarchical or otherwise). This uniformity can be explored by examining qq-plots (in which the trend should match a straight line) and plots of residual against the fitted values and each individual predictor (in which the noise should be uniform around zero across the range of x-values).

To illustrate this, lets generate 10 simulated data sets from our fitted model. This will generate a matrix with 10 columns and as many rows as there were in the original data. Think of it as 10 attempts to simulate the original data from the model.

dat.sim <- simulate(dat.glmL, n = 10)
dat.sim

   sim_1 sim_2 sim_3 sim_4 sim_5 sim_6 sim_7 sim_8 sim_9 sim_10
1      0     0     0     0     0     0     0     0     0      0
2      0     0     0     0     0     0     0     0     0      0
3      0     0     0     0     0     0     0     0     0      0
4      0     0     0     0     0     0     0     0     0      0
5      0     0     0     0     0     0     0     0     0      0
6      0     0     1     0     0     0     1     0     0      0
7      0     0     0     0     0     0     0     0     0      0
8      0     0     0     0     0     1     0     0     0      0
9      1     1     0     0     0     1     0     1     0      0
10     0     1     1     1     0     1     0     1     0      0
11     1     1     0     1     0     0     1     0     0      1
12     1     1     1     0     0     0     1     0     0      0
13     0     1     0     0     1     1     0     0     1      1
14     0     1     1     1     1     0     1     1     1      1
15     1     1     0     1     1     1     0     1     1      1
16     1     0     1     1     1     1     1     1     1      1
17     1     1     1     1     1     1     1     1     1      0
18     1     1     1     1     1     1     1     1     1      1
19     1     1     1     1     1     1     1     1     1      1
20     1     1     1     1     1     1     1     1     1      1

dat.sim <- simulate(dat.glmL, n = 250)
par(mfrow = c(5, 4), mar = c(3, 3, 1, 1))
resid <- NULL
for (i in 1:nrow(dat.sim)) {
    e = ecdf(data.matrix(dat.sim[i, ] + runif(250, -0.5, 0.5)))
    plot(e, main = i, las = 1)
    resid[i] <- e(dat$y[i] + runif(250, -0.5, 0.5))
}

resid

 [1] 0.836 0.296 0.760 0.236 0.872 0.232 0.040 0.240 0.756 0.356 0.132 0.864 0.972 0.528 0.056 0.248
[17] 0.264 0.076 0.800 0.640

plot(resid ~ fitted(dat.glmL))

The DHARMa package provides a number of convenient routines to explore standardized residuals simulated from fitted models based on the concepts outlined above. Along with generating simulated residuals, simple qq plots and residual plots are available. By default, the residual plots include quantile regression lines (0.25, 0.5 and 0.75), each of which should be straight and flat.

• in overdispersed models, the qq trend will deviate substantially from a straight line
• non-linear models will display trends in the residuals

library(DHARMa)

dat.sim <- simulateResiduals(dat.glmL)
dat.sim

Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
 
Scaled residual values: 0.492 0.664 0.216 0.952 0 0.84 0.592 0.66 0.988 0.48 0.456 0.9 0.596 0.692 0.028 0.988 0.944 0.456 0.436 0.176 ...

plotSimulatedResiduals(dat.sim)

Conclusions: there is no obvious patterns in the qq-plot or residuals, or at least there are no obvious trends remaining that would be indicative of overdispersion or non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmL, type = "pearson")^2)
  1 - pchisq(dat.resid, dat.glmL$df.resid)
  

  [1] 0.8571451
  

• Deviance ($G^2$) - similar to the $\chi^2$ test above, yet uses deviance

  1 - pchisq(dat.glmL$deviance, dat.glmL$df.resid)
  

  [1] 0.8647024
  

• The DHARMa package also has a routine for running a Kologorov-Smirnov test test to explore overall uniformity of the residuals as a goodness-of-fit test on the scaled residuals.

  testUniformity(dat.sim)
  

  	One-sample Kolmogorov-Smirnov test
  
  data:  simulationOutput$scaledResiduals
  D = 0.236, p-value = 0.2153
  alternative hypothesis: two-sided
  

The DHARMa package also provides elegant ways to explore overdispersion, zero-inflation (and spatial and temporal autocorrelation). For these methods, the model is repeatedly refit with the simulated data to yield bootstrapped refitted residuals. The test for overdispersion compares the approximate deviances of the observed model with the those of the simulated models.

testOverdispersion(simulateResiduals(dat.glmL, refit = T))

	DHARMa nonparametric overdispersion test via comparison to simulation under H0 = fitted
	model

data:  simulateResiduals(dat.glmL, refit = T)
dispersion = 1.1716, p-value = 0.388
alternative hypothesis: overdispersion

testZeroInflation(simulateResiduals(dat.glmL, refit = T))

	DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 =
	fitted model

data:  simulateResiduals(dat.glmL, refit = T)
ratioObsExp = 0.99171, p-value = 0.404
alternative hypothesis: more

In this demonstration, we fitted three logistic regressions (one for each link function). We could explore the residual plots of each of these for the purpose of comparing fit. We can also compare the fit of each of these three models via AIC (or deviance).

AIC(dat.glmL, dat.glmP, dat.glmC)

         df      AIC
dat.glmL  2 15.65141
dat.glmP  2 15.45223
dat.glmC  2 16.10823

Additional methods available can be explored via the following (substituting mod for the name of the fitted model).
methods(class=class(mod))

Lets explore the diagnostics - particularly the residuals

library(ggplot2)
ggplot(data = NULL) + geom_point(aes(y = residuals(dat.glmmTMBL, type = "pearson"), x = fitted(dat.glmmTMBL)))

qq.line = function(x) {
    # following four lines from base R's qqline()
    y <- quantile(x[!is.na(x)], c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1L] - slope * x[1L]
    return(c(int = int, slope = slope))
}
QQline = qq.line(resid(dat.glmmTMBL, type = "pearson"))
ggplot(data = NULL, aes(sample = resid(dat.glmmTMBL, type = "pearson"))) + stat_qq() + geom_abline(intercept = QQline[1],
    slope = QQline[2])

Unfortunately, unlike with linear models (Gaussian family), the expected distribution of data (residuals) varies over the range of fitted values for numerous (often competing) ways that make diagnosing (and attributing causes thereof) miss-specified generalized linear models from standard residual plots very difficult. The use of standardized (Pearson) residuals or deviance residuals can partly address this issue, yet they still do not offer completely consistent diagnoses across all issues (miss-specified model, over-dispersion, zero-inflation).

An alternative approach is to use simulated data from the fitted model to calculate an empirical cumulative density function from which residuals are are generated as values corresponding to the observed data along the density function. The (rather Bayesian) rationale is that if the model is correctly specified, then the observed data can be considered as a random draw from the fitted model. If this is the case, then residuals calculated from empirical cumulative density functions based on data simulated from the fitted model, should be totally uniform (flat) across the range of the linear predictor (fitted values) regardless of the model (Binomial, Poisson, linear, quadratic, hierarchical or otherwise). This uniformity can be explored by examining qq-plots (in which the trend should match a straight line) and plots of residual against the fitted values and each individual predictor (in which the noise should be uniform around zero across the range of x-values).

To illustrate this, lets generate 10 simulated data sets from our fitted model. This will generate a matrix with 10 columns and as many rows as there were in the original data. Think of it as 10 attempts to simulate the original data from the model.

dat.sim <- simulate(dat.glmmTMBL, n = 10)
dat.sim

   sim_1.1 sim_1.2 sim_2.1 sim_2.2 sim_3.1 sim_3.2 sim_4.1 sim_4.2 sim_5.1 sim_5.2 sim_6.1 sim_6.2
1        0       1       0       1       0       1       0       1       0       1       0       1
2        0       1       0       1       0       1       0       1       0       1       0       1
3        0       1       0       1       0       1       0       1       0       1       0       1
4        0       1       0       1       0       1       0       1       0       1       0       1
5        0       1       0       1       0       1       0       1       0       1       0       1
6        0       1       0       1       0       1       0       1       0       1       0       1
7        0       1       0       1       0       1       0       1       0       1       0       1
8        0       1       0       1       0       1       0       1       0       1       0       1
9        0       1       0       1       0       1       1       0       1       0       1       0
10       0       1       1       0       0       1       0       1       0       1       0       1
11       0       1       0       1       0       1       0       1       1       0       0       1
12       0       1       0       1       1       0       0       1       1       0       0       1
13       0       1       0       1       1       0       0       1       1       0       1       0
14       1       0       1       0       0       1       0       1       1       0       1       0
15       0       1       1       0       1       0       1       0       1       0       1       0
16       1       0       1       0       1       0       1       0       1       0       1       0
17       1       0       1       0       1       0       1       0       1       0       1       0
18       1       0       1       0       1       0       1       0       1       0       1       0
19       1       0       1       0       1       0       1       0       1       0       1       0
20       1       0       1       0       1       0       1       0       1       0       1       0
   sim_7.1 sim_7.2 sim_8.1 sim_8.2 sim_9.1 sim_9.2 sim_10.1 sim_10.2
1        0       1       0       1       0       1        0        1
2        0       1       0       1       0       1        0        1
3        0       1       0       1       0       1        0        1
4        0       1       0       1       0       1        0        1
5        0       1       0       1       0       1        0        1
6        0       1       0       1       0       1        0        1
7        0       1       0       1       0       1        0        1
8        0       1       0       1       0       1        0        1
9        0       1       1       0       0       1        0        1
10       0       1       1       0       0       1        0        1
11       0       1       1       0       1       0        0        1
12       0       1       1       0       1       0        0        1
13       0       1       1       0       0       1        0        1
14       1       0       0       1       1       0        1        0
15       1       0       1       0       1       0        0        1
16       1       0       1       0       1       0        1        0
17       1       0       1       0       1       0        1        0
18       1       0       1       0       1       0        1        0
19       1       0       1       0       1       0        1        0
20       1       0       1       0       1       0        1        0

For binomial glmmTMB models, there are two columns (one for success and one for failure) for each simulation. We only need one of these (the success field), so we will strip out every second column.

Now for each row of these simulated data, we calculate the empirical cumulative density function and use this function to predict new y-values (=residuals) corresponding to each observed y-value. Actually, 10 simulated samples is totally inadequate, we should use at least 250. I initially used 10, just so we could explore the output. We will re-simulate 250 times. Note also, for integer responses (including binary), uniform random noise is added to both the simulated and observed data so that we can sensibly explore zero inflation. The resulting residuals will be on a scale from 0 to 1 and therefore the residual plot should be centered around a y-value of 0.5.

dat.sim <- simulate(dat.glmmTMBL, n = 250)[, seq(1, 250 * 2, by = 2)]
par(mfrow = c(5, 4), mar = c(3, 3, 1, 1))
resid <- NULL
for (i in 1:nrow(dat.sim)) {
    e = ecdf(data.matrix(dat.sim[i, ] + runif(250, -0.5, 0.5)))
    plot(e, main = i, las = 1)
    resid[i] <- e(dat$y[i] + runif(250, -0.5, 0.5))
}

resid

 [1] 0.440 0.460 0.232 0.340 0.776 0.704 0.728 0.244 0.904 0.600 0.304 0.472 0.596 0.272 0.136 0.216
[17] 0.832 0.828 0.692 0.992

plot(resid ~ fitted(dat.glmmTMBL))

Conclusions: there is no obvious patterns in the qq-plot or residuals, or at least there are no obvious trends remaining that would be indicative of overdispersion or non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmmTMBL, type = "pearson")^2)
  1 - pchisq(dat.resid, df.residual(dat.glmmTMBL))
  

  [1] 0.8571451
  

• Deviance ($G^2$) - similar to the $\chi^2$ test above, yet uses deviance

  1 - pchisq(as.numeric(-2 * logLik(dat.glmmTMBL)), df.residual(dat.glmmTMBL))
  

  [1] 0.8647024
  

In this demonstration, we fitted three logistic regressions (one for each link function). We could explore the residual plots of each of these for the purpose of comparing fit. We can also compare the fit of each of these three models via AIC (or deviance).

AIC(dat.glmmTMBL, dat.glmmTMBP, dat.glmmTMBC)

             df      AIC
dat.glmmTMBL  2 15.65141
dat.glmmTMBP  2 15.45223
dat.glmmTMBC  2 16.10823

Lets explore the diagnostics - particularly the residuals

plot(dat.gamlssL)

******************************************************************
	 Summary of the Randomised Quantile Residuals
                           mean   =  0.1849072 
                       variance   =  1.272725 
               coef. of skewness  =  0.0372293 
               coef. of kurtosis  =  2.229781 
Filliben correlation coefficient  =  0.9830794 
******************************************************************

Unfortunately, unlike with linear models (Gaussian family), the expected distribution of data (residuals) varies over the range of fitted values for numerous (often competing) ways that make diagnosing (and attributing causes thereof) miss-specified generalized linear models from standard residual plots very difficult. The use of standardized (Pearson) residuals or deviance residuals can partly address this issue, yet they still do not offer completely consistent diagnoses across all issues (miss-specified model, over-dispersion, zero-inflation).

Conclusions: there is no obvious patterns in the qq-plot or residuals, or at least there are no obvious trends remaining that would be indicative of overdispersion or non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.gamlssL, type = "weighted", what = "mu")^2)
  1 - pchisq(dat.resid, dat.gamlssL$df.resid)
  

  [1] 0.8571451
  

• Deviance ($G^2$) - similar to the $\chi^2$ test above, yet uses deviance

  1 - pchisq(deviance(dat.gamlssL), dat.gamlssL$df.resid)
  

  [1] 0.8647024
  

In this demonstration, we fitted three logistic regressions (one for each link function). We could explore the residual plots of each of these for the purpose of comparing fit. We can also compare the fit of each of these three models via AIC (or deviance).

AIC(dat.gamlssL, dat.gamlssP, dat.gamlssC)

            df      AIC
dat.gamlssP  2 15.45223
dat.gamlssL  2 15.65141
dat.gamlssC  2 16.10823

summary(dat.glmL)

Call:
glm(formula = y ~ x, family = "binomial", data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.97157  -0.33665  -0.08191   0.30035   1.59628  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -6.9899     3.1599  -2.212   0.0270 *
x             0.6559     0.2936   2.234   0.0255 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 27.526  on 19  degrees of freedom
Residual deviance: 11.651  on 18  degrees of freedom
AIC: 15.651

Number of Fisher Scoring iterations: 6

exp(coef(dat.glmL))

(Intercept)           x 
0.000921113 1.926780616 

Conclusions:

• Intercept: when x is equal to zero, the odds of a success is 9 × 10-4 times greater than the odds of failure. Or the corollary, the odds of a failure is 1085.6432 times greater than the odds of success.
• Slope: for every 1 unit increase in x, the ratio odds of success to odds of failure changes by a factor of 1.9268. That is the odds ratio of success to failure nearly doubles with every 1 unit increase in x. From a inference testing perspective, we would reject the null hypothesis of no relationship.

summary(dat.glmmTMBL)

 Family: binomial  ( logit )
Formula:          y ~ x
Data: dat

     AIC      BIC   logLik deviance df.resid 
    15.7     17.6     -5.8     11.7       18 


Conditional model:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -6.9899     3.1599  -2.212   0.0270 *
x             0.6559     0.2936   2.234   0.0255 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

exp(fixef(dat.glmmTMBL)$cond)

 (Intercept)            x 
0.0009211129 1.9267806990 

Conclusions:

• Intercept: when x is equal to zero, the odds of a success is 9 × 10-4 times greater than the odds of failure. Or the corollary, the odds of a failure is 1085.6432 times greater than the odds of success.
• Slope: for every 1 unit increase in x, the ratio odds of success to odds of failure changes by a factor of 1.9268. That is the odds ratio of success to failure nearly doubles with every 1 unit increase in x. From a inference testing perspective, we would reject the null hypothesis of no relationship.

summary(dat.gamlssL)

******************************************************************
Family:  c("BI", "Binomial") 

Call:  gamlss(formula = y ~ x, family = BI(), data = dat) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  logit
Mu Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -6.9899     3.1598  -2.212   0.0401 *
x             0.6559     0.2936   2.234   0.0384 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  20 
Degrees of Freedom for the fit:  2
      Residual Deg. of Freedom:  18 
                      at cycle:  2 
 
Global Deviance:     11.65141 
            AIC:     15.65141 
            SBC:     17.64287 
******************************************************************

exp(coef(dat.gamlssL))

(Intercept)           x 
0.000921113 1.926780616 

Conclusions:

• Intercept: when x is equal to zero, the odds of a success is 9 × 10-4 times greater than the odds of failure. Or the corollary, the odds of a failure is 1085.6432 times greater than the odds of success.
• Slope: for every 1 unit increase in x, the ratio odds of success to odds of failure changes by a factor of 1.9268. That is the odds ratio of success to failure nearly doubles with every 1 unit increase in x. From a inference testing perspective, we would reject the null hypothesis of no relationship.

exp(coef(dat.glmL)[2])

       x 
1.926781 

# Or for arbitrary unit changes
library(oddsratio)
## 1 unit change
or_glm(dat, dat.glmL, inc = list(x = 1))

  predictor oddsratio CI_low (2.5 %) CI_high (97.5 %) increment
1         x     1.927          1.278            4.551         1

## 5 units change
or_glm(dat, dat.glmL, inc = list(x = 5))

  predictor oddsratio CI_low (2.5 %) CI_high (97.5 %) increment
1         x    26.556          3.407         1951.831         5

## Alternatively, we could express this as relative risk RR <- OR / (1 - P0 + (P0 * OR)) where P0 is
## the proportion of incidence in response
sjstats::or_to_rr(exp(coef(dat.glmL)[2]), mean(model.frame(dat.glmL)[[1]]))

       x 
1.359711 

sjstats::odds_to_rr(dat.glmL)

                    RR     lower.ci  upper.ci
(Intercept) 0.00167349 1.814187e-07 0.1349919
x           1.35971129 1.135840e+00 1.7517501

• the odds ratio is 1.9267806. For a one unit increase in x, the odds of success nearly double.
• the associated relative risk is 1.3597113. For a one unit increase in x, the probability of success increases by 36%. Every one unit change in x results in a 36% increase in risk of y being 1.

exp(fixef(dat.glmmTMBL)$cond[2])

       x 
1.926781 

## Alternatively, we could express this as relative risk RR <- OR / (1 - P0 + (P0 * OR)) where P0 is
## the proportion of incidence in response
sjstats::or_to_rr(exp(fixef(dat.glmmTMBL)$cond[2]), mean(model.frame(dat.glmmTMBL)[[1]]))

       x 
1.359711 

• the odds ratio is 1.9267806. For a one unit increase in x, the odds of success nearly double.
• the associated relative risk is 1.3597113. For a one unit increase in x, the probability of success increases by 36%. Every one unit change in x results in a 36% increase in risk of y being 1.

exp(coef(dat.gamlssL)[2])

       x 
1.926781 

## Alternatively, we could express this as relative risk RR <- OR / (1 - P0 + (P0 * OR)) where P0 is
## the proportion of incidence in response
sjstats::or_to_rr(exp(coef(dat.gamlssL)[2]), mean(model.frame(dat.gamlssL)[[1]]))

       x 
1.359711 

• the odds ratio is 1.9267806. For a one unit increase in x, the odds of success nearly double.
• the associated relative risk is 1.3597113. For a one unit increase in x, the probability of success increases by 36%. Every one unit change in x results in a 36% increase in risk of y being 1.

We perform the logistic regression using the glm() function. Clearly the number of successes is also dependent on the number of trials. Larger numbers of trials might be expected to yeild higher numbers of successes. Model weights can be used to account for differences in the total numbers of trials within each group. The most important (=commonly used) parameters/arguments for logistic regression are:

• formula: the linear model relating the response variable to the linear predictor
• family: specification of the error distribution (and link function). Can be specified as either a string or as one of the family functions (which allows for the explicit declaration of the link function). For examples:
  • family="binomial" or equivalently family=binomial(link="logit")
  • family=binomial(link="probit")
• data: the data frame containing the data
• weights: vector of weights to be used in the fitting process. These are particularly useful when dealing with grouped binary data from unequal sample sizes.

I will demonstrate logistic regression with a range of possible link functions (each of which yield different parameter interpretations):

• logit $$log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmL <- glm(cbind(success, failure) ~ x, data = dat, family = "binomial")
  # OR
  dat <- within(dat, total <- success + failure)
  dat.glmL <- glm(I(success/total) ~ x, data = dat, family = "binomial", weight = total)
  

• probit $$\phi^{-1}\left(\pi\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmP <- glm(cbind(success, failure) ~ x, data = dat,
      family = binomial(link = "probit"))
  # OR
  dat <- within(dat, total <- success + failure)
  dat.glmP <- glm(I(success/total) ~ x, data = dat, family = binomial(link = "probit"),
      weight = total)
  

• complimentary log-log $$log(-log(1-\pi))=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmC <- glm(cbind(success, failure) ~ x, data = dat,
      family = binomial(link = "cloglog"))
  # OR
  dat <- within(dat, total <- success + failure)
  dat.glmC <- glm(I(success/total) ~ x, data = dat, family = binomial(link = "cloglog"),
      weight = total)
  

We perform the logistic regression using the glmmTMB() function. The most important (=commonly used) parameters/arguments for logistic regression are:

• formula: the linear model relating the response variable to the linear predictor
• family: specification of the error distribution (and link function). Can be specified as either a string or as one of the family functions (which allows for the explicit declaration of the link function). For examples:
  • family="binomial" or equivalently family=binomial(link="logit")
  • family=binomial(link="probit")
• data: the data frame containing the data
• weights: vector of weights to be used in the fitting process. These are particularly useful when dealing with grouped binary data from unequal sample sizes (see below).

I will demonstrate logistic regression with a range of possible link functions (each of which yield different parameter interpretations):

• logit $$log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  library(glmmTMB)
  dat.glmmTMBL <- glmmTMB(cbind(success, failure) ~ x,
      data = dat, family = "binomial")
  

• probit $$\phi^{-1}\left(\pi\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmmTMBP <- glmmTMB(cbind(success, failure) ~ x,
      data = dat, family = binomial(link = "probit"))
  

• complimentary log-log $$log(-log(1-\pi))=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.glmmTMBC <- glmmTMB(cbind(success, failure) ~ x,
      data = dat, family = binomial(link = "cloglog"))
  

We perform the logistic regression using the gamlss() function. The most important (=commonly used) parameters/arguments for logistic regression are:

• formula: the linear model relating the response variable to the linear predictor
• family: specification of the error distribution (and link function). This is a different specification from glm() and supports a much wider variety of families, including mixtures. For examples:
  • family=BI() or equivalently family=BI(mu.link="logit")
  • family=BI(mu.link="probit")
• data: the data frame containing the data
• weights: vector of weights to be used in the fitting process. These are particularly useful when dealing with grouped binary data from unequal sample sizes (see below).

I will demonstrate logistic regression with a range of possible link functions (each of which yield different parameter interpretations):

• logit $$log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  library(gamlss)
  dat.gamlssL <- gamlss(cbind(success, failure) ~ x, data = dat,
      family = BI())
  

  GAMLSS-RS iteration 1: Global Deviance = 38.6222 
  GAMLSS-RS iteration 2: Global Deviance = 38.6222 
  

• probit $$\phi^{-1}\left(\pi\right)=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.gamlssP <- gamlss(cbind(success, failure) ~ x, data = dat,
      family = BI(mu.link = "probit"))
  

  GAMLSS-RS iteration 1: Global Deviance = 38.6267 
  GAMLSS-RS iteration 2: Global Deviance = 38.6267 
  

• complimentary log-log $$log(-log(1-\pi))=\beta_0+\beta_1x_1+\epsilon_i, \hspace{1cm}\epsilon\sim Binom(\lambda)$$

  dat.gamlssC <- gamlss(cbind(success, failure) ~ x, data = dat,
      family = BI(mu.link = "cloglog"))
  

  GAMLSS-RS iteration 1: Global Deviance = 38.5955 
  GAMLSS-RS iteration 2: Global Deviance = 38.5955 
  

plot(resid(dat.glmL) ~ fitted(dat.glmL))
lines(lowess(resid(dat.glmL) ~ fitted(dat.glmL)))

dat.sim <- simulateResiduals(dat.glmL)
plotSimulatedResiduals(dat.sim)

Conclusions: there is no obvious patterns in the residuals, or at least there are no obvious trends remaining that would be indicative of non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmL, type = "pearson")^2)
  1 - pchisq(dat.resid, dat.glmL$df.resid)
  

  [1] 0.6763891
  

• Deviance ($G^2$) - similar to the $chi^2$ test above, yet uses deviance

  1 - pchisq(dat.glmL$deviance, dat.glmL$df.resid)
  

  [1] 0.6646471
  

• using the DHARMa package

  testUniformity(dat.sim)
  

  	One-sample Kolmogorov-Smirnov test
  
  data:  simulationOutput$scaledResiduals
  D = 0.18, p-value = 0.8473
  alternative hypothesis: two-sided
  

In the absence of any evidence of a goodness of the fit of the model, there is unlikely to be any issues with overdispersion. However, it is worth estimating the dispersion parameter nonetheless.

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmL, type = "pearson")^2)
  dat.resid/dat.glmL$df.resid
  

  [1] 0.7174332
  

• Deviance ($G^2$) - similar to the $chi^2$ test above, yet uses deviance

  dat.glmL$deviance/dat.glmL$df.resid
  

  [1] 0.7305604
  

• using the DHARMa package

  testOverdispersion(simulateResiduals(dat.glmL, refit = T))
  

  Error in ecdf(out$refittedResiduals[i, ] + runif(out$nSim, -0.5, 0.5)): 'x' must have 1 or more non-missing values
  

  testZeroInflation(simulateResiduals(dat.glmL, refit = T))
  

  Error in ecdf(out$refittedResiduals[i, ] + runif(out$nSim, -0.5, 0.5)): 'x' must have 1 or more non-missing values
  

In this demonstration, we fitted three logistic regressions (one for each link function). We could explore the residual plots of each of these for the purpose of comparing fit. We can also compare the fit of each of these three models via AIC (or deviance).

AIC(dat.glm, dat.glmP, dat.glmC)

Error in AIC(dat.glm, dat.glmP, dat.glmC): object 'dat.glm' not found

library(ggplot2)
ggplot(data = NULL) + geom_point(aes(y = residuals(dat.glmmTMBL,
    type = "pearson"), x = fitted(dat.glmmTMBL)))

qq.line = function(x) {
    # following four lines from base R's qqline()
    y <- quantile(x[!is.na(x)], c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1L] - slope * x[1L]
    return(c(int = int, slope = slope))
}
QQline = qq.line(resid(dat.glmmTMBL, type = "pearson"))
ggplot(data = NULL, aes(sample = resid(dat.glmmTMBL, type = "pearson"))) +
    stat_qq() + geom_abline(intercept = QQline[1], slope = QQline[2])

dat.sim <- simulate(dat.glmmTMBL, n = 250)[, seq(1, 250 * 2,
    by = 2)]
par(mfrow = c(4, 3), mar = c(3, 3, 1, 1))
resid <- NULL
for (i in 1:nrow(dat.sim)) {
    e = ecdf(data.matrix(dat.sim[i, ] + runif(250, -0.5, 0.5)))
    plot(e, main = i, las = 1)
    resid[i] <- e(dat$success[i] + runif(250, -0.5, 0.5))
}
resid

 [1] 0.756 0.364 0.172 0.596 0.572 0.920 0.040 0.492 0.344 0.632

plot(resid ~ fitted(dat.glmmTMBL))

Conclusions: there is no obvious patterns in the residuals, or at least there are no obvious trends remaining that would be indicative of non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmmTMBL, type = "pearson")^2)
  1 - pchisq(dat.resid, df.residual(dat.glmmTMBL))
  

  [1] 0.6763891
  

In the absence of any evidence of a goodness of the fit of the model, there is unlikely to be any issues with overdispersion. However, it is worth estimating the dispersion parameter nonetheless.

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmmTMBL, type = "pearson")^2)
  dat.resid/df.residual(dat.glmmTMBL)
  

  [1] 0.7174332
  

In this demonstration, we fitted three logistic regressions (one for each link function). We could explore the residual plots of each of these for the purpose of comparing fit. We can also compare the fit of each of these three models via AIC (or deviance).

AIC(dat.glmmTMBL, dat.glmmTMBP, dat.glmmTMBC)

             df      AIC
dat.glmmTMBL  2 42.62224
dat.glmmTMBP  2 42.62673
dat.glmmTMBC  2 42.59548

plot(dat.gamlssL)

******************************************************************
	 Summary of the Randomised Quantile Residuals
                           mean   =  0.0461986 
                       variance   =  0.4918538 
               coef. of skewness  =  -0.04136291 
               coef. of kurtosis  =  2.284877 
Filliben correlation coefficient  =  0.9844053 
******************************************************************

Unfortunately, unlike with linear models (Gaussian family), the expected distribution of data (residuals) varies over the range of fitted values for numerous (often competing) ways that make diagnosing (and attributing causes thereof) miss-specified generalized linear models from standard residual plots very difficult. The use of standardized (Pearson) residuals or deviance residuals can partly address this issue, yet they still do not offer completely consistent diagnoses across all issues (miss-specified model, over-dispersion, zero-inflation).

Conclusions: there is no obvious patterns in the qq-plot or residuals, or at least there are no obvious trends remaining that would be indicative of overdispersion or non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.gamlssL, type = "weighted", what = "mu")^2)
  1 - pchisq(dat.resid, dat.gamlssL$df.resid)
  

  [1] 0.6763891
  

In this demonstration, we fitted three logistic regressions (one for each link function). We could explore the residual plots of each of these for the purpose of comparing fit. We can also compare the fit of each of these three models via AIC (or deviance).

AIC(dat.gamlssL, dat.gamlssP, dat.gamlssC)

            df      AIC
dat.gamlssC  2 42.59548
dat.gamlssL  2 42.62224
dat.gamlssP  2 42.62673

summary(dat.glmL)

Call:
glm(formula = I(success/total) ~ x, family = "binomial", data = dat, 
    weights = total)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.45415  -0.32619   0.07158   0.45284   1.43145  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.01854    1.08802   3.693 0.000221 ***
x           -0.25622    0.06803  -3.766 0.000166 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 21.0890  on 9  degrees of freedom
Residual deviance:  5.8445  on 8  degrees of freedom
AIC: 42.622

Number of Fisher Scoring iterations: 3

exp(coef(dat.glmL))

(Intercept)           x 
 55.6200133   0.7739699 

library(broom)
tidy(dat.glmL)

# A tibble: 2 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    4.02     1.09        3.69 0.000221
2 x             -0.256    0.0680     -3.77 0.000166

glance(dat.glmL)

# A tibble: 1 x 7
  null.deviance df.null logLik   AIC   BIC deviance df.residual
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
1          21.1       9  -19.3  42.6  43.2     5.84           8

Conclusions:

• Intercept: when x is equal to zero, the odds of a success is 55.62 times greater than the odds of failure.
• Slope: for every 1 unit increase in x, the ratio odds of success to odds of failure changes by a factor of 0.774. That is the odds ratio of success to failure declines at a rate of nearly 3/4 for every 1 unit increase in x. From a inference testing perspective, we would reject the null hypothesis of no relationship.

summary(dat.glmmTMBL)

 Family: binomial  ( logit )
Formula:          cbind(success, failure) ~ x
Data: dat

     AIC      BIC   logLik deviance df.resid 
    42.6     43.2    -19.3     38.6        8 


Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.01854    1.08804   3.693 0.000221 ***
x           -0.25622    0.06803  -3.766 0.000166 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

exp(fixef(dat.glmmTMBL)$cond)

(Intercept)           x 
 55.6200323   0.7739699 

library(broom)
library(broom.mixed)  # devtools::install_github('bbolker/broom.mixed')
tidy(dat.glmmTMBL)

  effect component group        term   estimate  std.error statistic      p.value
1  fixed      cond fixed (Intercept)  4.0185434 1.08803755  3.693387 0.0002212871
2  fixed      cond fixed           x -0.2562223 0.06803094 -3.766261 0.0001657106

glance(dat.glmmTMBL)

  sigma    logLik      AIC      BIC df.residual
1     1 -19.31112 42.62224 43.22741           8

Conclusions:

• Intercept: when x is equal to zero, the odds of a success is 55.62 times greater than the odds of failure.
• Slope: for every 1 unit increase in x, the ratio odds of success to odds of failure changes by a factor of 0.774. That is the odds ratio of success to failure declines at a rate of nearly 3/4 for every 1 unit increase in x. From a inference testing perspective, we would reject the null hypothesis of no relationship.

summary(dat.gamlssL)

******************************************************************
Family:  c("BI", "Binomial") 

Call:  gamlss(formula = cbind(success, failure) ~ x, family = BI(),      data = dat) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  logit
Mu Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  4.01854    1.08805   3.693   0.0061 **
x           -0.25622    0.06803  -3.766   0.0055 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  10 
Degrees of Freedom for the fit:  2
      Residual Deg. of Freedom:  8 
                      at cycle:  2 
 
Global Deviance:     38.62224 
            AIC:     42.62224 
            SBC:     43.22741 
******************************************************************

exp(coef(dat.gamlssL))

(Intercept)           x 
 55.6200134   0.7739699 

library(broom)
tidy(dat.gamlssL)

  parameter        term   estimate std.error statistic     p.value
1        mu (Intercept)  4.0185431 1.0880241  3.693432 0.006099988
2        mu           x -0.2562222 0.0680295 -3.766340 0.005494302

glance(dat.gamlssL)

# A tibble: 1 x 6
     df logLik   AIC   BIC deviance df.residual
  <int>  <dbl> <dbl> <dbl>    <dbl>       <dbl>
1     0  -19.3  42.6  43.2     38.6          8.

Conclusions:

• Intercept: when x is equal to zero, the odds of a success is 55.62 times greater than the odds of failure.
• Slope: for every 1 unit increase in x, the ratio odds of success to odds of failure changes by a factor of 0.774. That is the odds ratio of success to failure declines at a rate of nearly 3/4 for every 1 unit increase in x. From a inference testing perspective, we would reject the null hypothesis of no relationship.

A measure of the strength of the relationship can be obtained according to: $$quasi R^2 = 1-\left(\frac{deviance}{null~deviance}\right)$$

1 - (dat.glmL$deviance/dat.glmL$null)

[1] 0.722866

We might also be interested in the LD50 - the value of $x$ where the probability switches from favoring 1 to favoring 0. LD50 is calculated as: $$LD50 = - \frac{intercept}{slope}$$

-dat.glmL$coef[1]/dat.glmL$coef[2]

(Intercept) 
   15.68382 

A measure of the strength of the relationship can be obtained according to: $$quasi R^2 = 1-\left(\frac{deviance}{null~deviance}\right)$$

1 - exp((2/nrow(dat)) * (logLik(update(dat.glmmTMBL, ~1))[1] - logLik(dat.glmmTMBL)[1]))

[1] 0.7822599

# This does not seem to work for proportional data.  The formula above relies on deviance being
# residual deviance that is -2*log(L/L_0) rather than -2*log(L)
1 - (as.numeric(-2 * logLik(dat.glmmTMBL))/as.numeric(-2 * logLik(update(dat.glmmTMBL, ~1))))

[1] 0.2830044

MuMIn::r.squaredGLMM(dat.glmmTMBL)  # only seems to work when it is a mixed model

Error in reformulate(vapply(lapply(.findbars(form), "[[", 2L), deparse, : 'termlabels' must be a character vector of length at least one

We might also be interested in the LD50 - the value of $x$ where the probability switches from favoring 1 to favoring 0. LD50 is calculated as: $$LD50 = - \frac{intercept}{slope}$$

-fixef(dat.glmmTMBL)$cond[1]/fixef(dat.glmmTMBL)$cond[2]

(Intercept) 
   15.68382 

A measure of the strength of the relationship can be obtained according to: $$quasi R^2 = 1-\left(\frac{deviance}{null~deviance}\right)$$

Rsq(dat.gamlssL)

[1] 0.7822599

We might also be interested in the LD50 - the value of $x$ where the probability switches from favoring 1 to favoring 0. LD50 is calculated as: $$LD50 = - \frac{intercept}{slope}$$

-coef(dat.gamlssL)[1]/coef(dat.gamlssL)[2]

(Intercept) 
   15.68382 

plot(resid(dat.glm) ~ fitted(dat.glm))
lines(lowess(resid(dat.glm) ~ fitted(dat.glm)))

dat.sim <- simulateResiduals(dat.glm)
plotSimulatedResiduals(dat.sim)

Conclusions: there is no obvious patterns in the residuals, or at least there are no obvious trends remaining that would be indicative of non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glm, type = "pearson")^2)
  1 - pchisq(dat.resid, dat.glm$df.resid)
  

  [1] 0.9999987
  

• Deviance ($G^2$) - similar to the $chi^2$ test above, yet uses deviance

  1 - pchisq(dat.glm$deviance, dat.glm$df.resid)
  

  [1] 0.9999983
  

• using the DHARMa package

  testUniformity(dat.sim)
  

  	One-sample Kolmogorov-Smirnov test
  
  data:  simulationOutput$scaledResiduals
  D = 0.168, p-value = 0.5681
  alternative hypothesis: two-sided
  

In the absence of any evidence of a goodness of the fit of the model, there is unlikely to be any issues with overdispersion. However, it is worth estimating the dispersion parameter nonetheless.

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glm, type = "pearson")^2)
  dat.resid/dat.glm$df.resid
  

  [1] 0.1131458
  

• Deviance ($G^2$) - similar to the $chi^2$ test above, yet uses deviance

  dat.glm$deviance/dat.glm$df.resid
  

  [1] 0.1172321
  

• using the DHARMa package

  testOverdispersion(simulateResiduals(dat.glm, refit = T))
  

  	DHARMa nonparametric overdispersion test via comparison to simulation under H0 = fitted
  	model
  
  data:  simulateResiduals(dat.glm, refit = T)
  dispersion = 0.14353, p-value = 0.94
  alternative hypothesis: overdispersion
  

  testZeroInflation(simulateResiduals(dat.glm, refit = T))
  

  

  	DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 =
  	fitted model
  
  data:  simulateResiduals(dat.glm, refit = T)
  ratioObsExp = 0.18484, p-value = 1
  alternative hypothesis: more
  

plot(resid(dat.betareg) ~ fitted(dat.betareg))
lines(lowess(resid(dat.betareg) ~ fitted(dat.betareg)))

Conclusions: there is no obvious patterns in the residuals, or at least there are no obvious trends remaining that would be indicative of non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.betareg, type = "pearson")^2)
  1 - pchisq(dat.resid, dat.betareg$df.resid)
  

  [1] 0.2952528
  

In the absence of any evidence of a goodness of the fit of the model, there is unlikely to be any issues with overdispersion. However, it is worth estimating the dispersion parameter nonetheless.

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.betareg, type = "pearson")^2)
  dat.resid/dat.betareg$df.resid
  

  [1] 1.152881
  

library(ggplot2)
ggplot(data = NULL) + geom_point(aes(y = residuals(dat.glmmTMB,
    type = "pearson"), x = fitted(dat.glmmTMB)))

qq.line = function(x) {
    # following four lines from base R's qqline()
    y <- quantile(x[!is.na(x)], c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1L] - slope * x[1L]
    return(c(int = int, slope = slope))
}
QQline = qq.line(resid(dat.glmmTMB, type = "pearson"))
ggplot(data = NULL, aes(sample = resid(dat.glmmTMB, type = "pearson"))) +
    stat_qq() + geom_abline(intercept = QQline[1], slope = QQline[2])

dat.sim <- simulate(dat.glmmTMB, n = 250)
par(mfrow = c(5, 4), mar = c(3, 3, 1, 1))
resid <- NULL
for (i in 1:nrow(dat.sim)) {
    e = ecdf(data.matrix(dat.sim[i, ] + runif(250, -0.5, 0.5)))
    plot(e, main = i, las = 1)
    resid[i] <- e(dat$y[i] + runif(250, -0.5, 0.5))
}

resid

 [1] 0.572 0.664 0.472 0.204 0.532 0.748 0.152 0.256 0.368 0.736 0.840 0.524 0.492 0.700 0.808 0.336
[17] 1.000 0.616 0.596 0.328

par(mfrow = c(1, 1))
plot(resid ~ fitted(dat.glmmTMB))

Conclusions: there is no obvious patterns in the residuals, or at least there are no obvious trends remaining that would be indicative of non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmmTMB, type = "pearson")^2)
  1 - pchisq(dat.resid, df.residual(dat.glmmTMB))
  

  [1] 0.999998
  

In the absence of any evidence of a goodness of the fit of the model, there is unlikely to be any issues with overdispersion. However, it is worth estimating the dispersion parameter nonetheless.

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.glmmTMB, type = "pearson")^2)
  dat.resid/df.residual(dat.glmmTMB)
  

  [1] 0.1098094
  

plot(dat.gamlssL)

******************************************************************
	 Summary of the Randomised Quantile Residuals
                           mean   =  0.0461986 
                       variance   =  0.4918538 
               coef. of skewness  =  -0.04136291 
               coef. of kurtosis  =  2.284877 
Filliben correlation coefficient  =  0.9844053 
******************************************************************

Conclusions: there is no obvious patterns in the qq-plot or residuals, or at least there are no obvious trends remaining that would be indicative of overdispersion or non-linearity.

We can also explore the goodness of the fit of the model via:

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.gamlssL, type = "weighted", what = "mu")^2)
  1 - pchisq(dat.resid, dat.gamlssL$df.resid)
  

  [1] 0.6763891
  

In the absence of any evidence of a goodness of the fit of the model, there is unlikely to be any issues with overdispersion. However, it is worth estimating the dispersion parameter nonetheless.

• Pearson's $\chi^2$ residuals - explores whether there are any significant patterns remaining in the residuals

  dat.resid <- sum(resid(dat.gamlss, type = "weighted", what = "mu")^2)
  dat.resid/df.residual(dat.gamlss)
  

  [1] 0.9089873
  

summary(dat.glm)

Call:
glm(formula = y ~ x, family = "binomial", data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.70852  -0.17005  -0.03006   0.25817   0.58884  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -4.8626     2.0981  -2.318   0.0205 *
x             0.4588     0.1936   2.369   0.0178 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 14.1513  on 19  degrees of freedom
Residual deviance:  2.1102  on 18  degrees of freedom
AIC: 13.985

Number of Fisher Scoring iterations: 5

binomial()$linkinv(coef(dat.glm))

(Intercept)           x 
0.007671185 0.612720691 

library(broom)
tidy(dat.glm)

# A tibble: 2 x 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   -4.86      2.10      -2.32  0.0205
2 x              0.459     0.194      2.37  0.0178

glance(dat.glm)

# A tibble: 1 x 7
  null.deviance df.null logLik   AIC   BIC deviance df.residual
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
1          14.2      19  -4.99  14.0  16.0     2.11          18

Conclusions:

• Intercept: when x is equal to zero, the estimated percent cover is 0.0077.
• Slope: for every 1 unit increase in x, the percentage cover changes by 0.6127.

summary(dat.betareg)

Call:
betareg(formula = y2 ~ x, data = dat)

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-1.7681 -0.5719 -0.2479  0.7006  2.2933 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.30109    0.50582  -8.503   <2e-16 ***
x            0.41680    0.04713   8.843   <2e-16 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)   
(phi)    9.499      3.163   3.003  0.00268 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 27.05 on 3 Df
Pseudo R-squared: 0.8771
Number of iterations: 22 (BFGS) + 3 (Fisher scoring) 

dat.betareg$link$mean$linkinv(coef(dat.betareg))

(Intercept)           x       (phi) 
 0.01337254  0.60271817  0.99992507 

library(broom)
tidy(dat.betareg)

# A tibble: 3 x 6
  component term        estimate std.error statistic  p.value
  <chr>     <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 mean      (Intercept)   -4.30     0.506      -8.50 1.84e-17
2 mean      x              0.417    0.0471      8.84 9.31e-19
3 precision (phi)          9.50     3.16        3.00 2.68e- 3

glance(dat.betareg)

# A tibble: 1 x 6
  pseudo.r.squared df.null logLik   AIC   BIC df.residual
             <dbl>   <dbl>  <dbl> <dbl> <dbl>       <int>
1            0.877     18.   27.0 -48.1 -45.1          17

Conclusions:

• Intercept: when x is equal to zero, the estimated percent cover is 0.0134.
• Slope: for every 1 unit increase in x, the percentage cover changes by 0.6027.

summary(dat.glmmTMB)

 Family: beta  ( logit )
Formula:          y2 ~ x
Data: dat

     AIC      BIC   logLik deviance df.resid 
   -48.1    -45.1     27.0    -54.1       17 


Overdispersion parameter for beta family ():  9.5 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.30109    0.51895  -8.288   <2e-16 ***
x            0.41680    0.04792   8.699   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

family(dat.glmmTMB)$linkinv(fixef(dat.glmmTMB)$cond)

(Intercept)           x 
 0.01337254  0.60271817 

library(broom)
library(broom.mixed)  # devtools::install_github('bbolker/broom.mixed')
tidy(dat.glmmTMB)

  effect component group        term   estimate  std.error statistic      p.value
1  fixed      cond fixed (Intercept) -4.3010894 0.51894614 -8.288123 1.150497e-16
2  fixed      cond fixed           x  0.4168038 0.04791626  8.698587 3.360419e-18

glance(dat.glmmTMB)

    sigma   logLik       AIC       BIC df.residual
1 9.49893 27.04501 -48.09003 -45.10283          17

Conclusions:

• Intercept: when x is equal to zero, the estimated percent cover is 0.0077.
• Slope: for every 1 unit increase in x, the percentage cover changes by 0.6127.

summary(dat.gamlss)

******************************************************************
Family:  c("BEINF", "Beta Inflated") 

Call:  gamlss(formula = y ~ x, family = BEINF(), data = dat) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  logit
Mu Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.03415    0.65895  -6.122 1.95e-05 ***
x            0.38129    0.06634   5.748 3.85e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
Sigma link function:  logit
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -0.7222     0.2464  -2.931   0.0103 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
Nu link function:  log 
Nu Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -1.9459     0.7559  -2.574   0.0212 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
Tau link function:  log 
Tau Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -1.2528     0.5669   -2.21   0.0431 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  20 
Degrees of Freedom for the fit:  5
      Residual Deg. of Freedom:  15 
                      at cycle:  5 
 
Global Deviance:     7.92312 
            AIC:     17.92312 
            SBC:     22.90178 
******************************************************************

binomial()$linkinv(coef(dat.gamlss))  # logit inverse

(Intercept)           x 
 0.01739277  0.59418538 

library(broom)
tidy(dat.gamlss)

  parameter        term   estimate  std.error statistic      p.value
1        mu (Intercept) -4.0341548 0.62893393 -6.414274 4.879024e-06
2        mu           x  0.3812949 0.06322293  6.030959 1.057823e-05
3     sigma (Intercept) -0.7222413 0.22820075 -3.164938 5.099291e-03
4        nu (Intercept) -1.9459101 0.74535591 -2.610713 1.718756e-02
5       tau (Intercept) -1.2527630 0.55901699 -2.241011 3.716592e-02

glance(dat.gamlss)

# A tibble: 1 x 6
     df logLik   AIC   BIC deviance df.residual
  <int>  <dbl> <dbl> <dbl>    <dbl>       <dbl>
1     0  -3.96  17.9  22.9     7.92         15.

Conclusions:

• Intercept: when x is equal to zero, the estimated percent cover is 0.0174.
• Slope: for every 1 unit increase in x, the percentage cover changes by 0.5942.

A measure of the strength of the relationship can be obtained according to: $$quasi R^2 = 1-\left(\frac{deviance}{null~deviance}\right)$$

1 - exp((2/nrow(dat)) * (logLik(update(dat.glmmTMB, ~1))[1] - logLik(dat.glmmTMB)[1]))

[1] 0.8679458

MuMIn::r.squaredGLMM(dat.glmmTMB)  # only seems to work when it is a mixed model

Error in reformulate(vapply(lapply(.findbars(form), "[[", 2L), deparse, : 'termlabels' must be a character vector of length at least one