knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(sclr)
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) # Two-titre fit to included simulated data fit2 <- sclr(status ~ logHI + logNI, two_titre_data) summary(fit2)
The predict
method will return the point estimate and a confidence
interval of
$\beta_0 + \beta_1X_1 + ... + \beta_kX_k$
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")]) # 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")])
To 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) prot_lvls2 <- data.frame(logNI = log(c(0.1, 10, 40))) protHI2 <- get_protection_level(fit2, "logHI", prot_lvls2) print(protHI2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.