Let’s consider the following data:

library(tidyverse)

n <- 100
x <- c(.2, .4, .6, .8, 1)
k <- c(10, 26, 73, 94, 97)
prop <- k / n

dat <- tibble(x, k, n, prop)

dat |> 
  ggplot(aes(x, prop)) +
  geom_point()

Let’s fit a psychometric function to this data maximizing the likelihood function.

create_log_lik_fun <- function(x, k, n) {
  function(p) {
    phi <- 1 / (1 + exp(-p[1] - p[2] * x))
    sum(k * log(phi) + (n - k) * log(1 - phi))
  }
}

log_lik_fun <- create_log_lik_fun(dat$x, dat$k, dat$n)

neg_log_lik_fun <- function(p) -log_lik_fun(p)

fit <- optim(c(0, 0), neg_log_lik_fun)

p_mle <- fit$par

dat |> 
  ggplot(aes(x, prop)) +
  geom_point() +
  geom_line(data = tibble(x = seq(0, 1, .01)) |> 
              mutate(prop = 1 / (1 + exp(-p_mle[1] - p_mle[2] * x))))

Goodness of fit

To assess the goodness of fit, we can use the ratio of the likelihood of the proposed model evaluated in the best MLE parameters relative to the likelihood of the saturated model evaluated in the best MLE parameters (see MLE for the binomial distribution and psychometric functions).

\[\Lambda = \frac{L_{\text{proposed}}}{L_{\text{saturated}}}\] The log likelihood of the proposed model is

log_lik_proposed <- log_lik_fun(p_mle)

log_lik_proposed
## [1] -186.1255

The log likelihood of the saturated model is

create_log_lik_saturated_fun <- function(x, k, n) {
  function(p) {
    # p <- pmax(pmin(p, 0.9999), 0.0001) # Very important if optimization is performed to avoid log(0) and log(1) !!
    sum(k * log(p) + (n - k) * log(1 - p))
  }
}

log_lik_saturated_fun <- create_log_lik_saturated_fun(dat$x, dat$k, dat$n)

log_lik_saturated <- log_lik_saturated_fun(dat$prop)

log_lik_saturated
## [1] -184.3108

We calculate de deviance and estimate the p-value using Wilks theorem

deviance <- -2 * (log_lik_proposed - log_lik_saturated)

deviance
## [1] 3.629222
pchisq(deviance, length(dat$prop) - 2, lower.tail = FALSE)
## [1] 0.3043851

We do not have evidence to reject the proposed model.

Calculating deviance using glm

The deviance is called Residual deviance

model <- glm(cbind(k, n - k) ~ x, data = dat, family = binomial(logit))

model
## 
## Call:  glm(formula = cbind(k, n - k) ~ x, family = binomial(logit), 
##     data = dat)
## 
## Coefficients:
## (Intercept)            x  
##      -4.080        8.216  
## 
## Degrees of Freedom: 4 Total (i.e. Null);  3 Residual
## Null Deviance:       304.4 
## Residual Deviance: 3.629     AIC: 27.87

Assessing whether the predictor is necessary

To assess whether the predictor is necessary, we can use the ratio of the likelihood of the null model (model with just an intercept) relative to the saturated model.

\[\Lambda = \frac{L_{\text{null}}}{L_{\text{saturated}}}\]

create_log_lik_null_fun <- function(k, n) {
  function(p) {
    sum(k * log(p) + (n - k) * log(1 - p))
  }
}

log_lik_null_fun <- create_log_lik_null_fun(dat$k, dat$n)

neg_log_lik_null_fun <- function(p) -log_lik_null_fun(p)

fit_null <- optim(.5, neg_log_lik_null_fun)
## Warning in optim(0.5, neg_log_lik_null_fun): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
p_null <- fit_null$par

log_lik_null <- log_lik_null_fun(p_null)

Then the null deviance is

deviance_null <- -2 * (log_lik_null - log_lik_saturated)

deviance_null
## [1] 304.39

which is the result that we obtain using glm

model
## 
## Call:  glm(formula = cbind(k, n - k) ~ x, family = binomial(logit), 
##     data = dat)
## 
## Coefficients:
## (Intercept)            x  
##      -4.080        8.216  
## 
## Degrees of Freedom: 4 Total (i.e. Null);  3 Residual
## Null Deviance:       304.4 
## Residual Deviance: 3.629     AIC: 27.87

To obtain the p value, we can use Wilks theorem

pchisq(deviance_null, length(dat$prop) - 1, lower.tail = FALSE)
## [1] 1.224051e-64