See vignette("sclr-math") for model specification,
log-likelihood, scores and second derivatives.
Model fitting is done with a function called sclr. It is
used in the same way as other fitting functions like
lm.
# One-titre fit to included simulated data
fit1 <- sclr(status ~ logHI, one_titre_data)
summary(fit1)
#> Call: status ~ logHI
#>
#> Parameter estimates
#> theta beta_0 beta_logHI
#> -0.03497876 -5.42535734 2.14877741
#>
#> 95% confidence intervals
#> 2.5 % 97.5 %
#> theta -0.1350572 0.06509969
#> beta_0 -6.4417802 -4.40893449
#> beta_logHI 1.8146909 2.48286390
#>
#> Log likelihood: -2469.765
# Two-titre fit to included simulated data
fit2 <- sclr(status ~ logHI + logNI, two_titre_data)
summary(fit2)
#> Call: status ~ logHI + logNI
#>
#> Parameter estimates
#> theta beta_0 beta_logHI beta_logNI
#> -0.1231823 -5.6536038 2.2098460 1.0409995
#>
#> 95% confidence intervals
#> 2.5 % 97.5 %
#> theta -0.2395860 -0.006778714
#> beta_0 -6.6252070 -4.682000648
#> beta_logHI 1.8981507 2.521541191
#> beta_logNI 0.8738325 1.208166538
#>
#> Log likelihood: -1858.46The predict method will return the point estimate and a
confidence interval of β0 + β1X1 + ... + βkXk
where k is the number of
covariates. It will also apply the inverse logit transformation to these
estimates and interval bounds to get the point estimate and the interval
for the probability of protection on the original scale.
# One-titre fit
preddata1 <- data.frame(logHI = seq(0, 8, length.out = 101))
pred1 <- predict(fit1, preddata1)
head(pred1[, c("logHI", "prot_l", "prot_point", "prot_u")])
#> # A tibble: 6 × 4
#> logHI prot_l prot_point prot_u
#> <dbl> <dbl> <dbl> <dbl>
#> 1 0 0.00159 0.00438 0.0120
#> 2 0.08 0.00194 0.00520 0.0139
#> 3 0.16 0.00236 0.00617 0.0160
#> 4 0.24 0.00288 0.00732 0.0185
#> 5 0.32 0.00351 0.00868 0.0213
#> 6 0.4 0.00427 0.0103 0.0246
# Two-titre fit
preddata2 <- data.frame(logHI = seq(0, 8, length.out = 101), logNI = 1)
pred2 <- predict(fit2, preddata2)
head(pred2[, c("logHI", "logNI", "prot_l", "prot_point", "prot_u")])
#> # A tibble: 6 × 5
#> logHI logNI prot_l prot_point prot_u
#> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0 1 0.00432 0.00983 0.0222
#> 2 0.08 1 0.00527 0.0117 0.0258
#> 3 0.16 1 0.00643 0.0139 0.0299
#> 4 0.24 1 0.00785 0.0166 0.0347
#> 5 0.32 1 0.00957 0.0197 0.0402
#> 6 0.4 1 0.0117 0.0235 0.0466To get the estimated titre (and the confidence interval) that
corresponds to a particular protection level (eg. 50%), use the
get_protection_level function. Its interface is similar to
that of predict.
protHI1 <- get_protection_level(fit1, "logHI", lvl = 0.5)
print(protHI1)
#> # A tibble: 3 × 3
#> logHI prot_prob est
#> <dbl> <dbl> <chr>
#> 1 2.39 0.5 low bound
#> 2 2.52 0.5 point
#> 3 2.64 0.5 upper bound
prot_lvls2 <- data.frame(logNI = log(c(0.1, 10, 40)))
protHI2 <- get_protection_level(fit2, "logHI", prot_lvls2)
print(protHI2)
#> logNI logHI prot_prob est
#> 1 -2.302585 3.4024635 0.5 low bound
#> 2 2.302585 1.3380153 0.5 low bound
#> 3 3.688879 0.6568784 0.5 low bound
#> 4 -2.302585 3.6430566 0.5 point
#> 5 2.302585 1.4736837 0.5 point
#> 6 3.688879 0.8206373 0.5 point
#> 7 -2.302585 3.8723188 0.5 upper bound
#> 8 2.302585 1.5883861 0.5 upper bound
#> 9 3.688879 0.9605298 0.5 upper bound