---
title: "Usage"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Usage}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(sclr)
```

## Model

See `vignette("sclr-math")` for model specification, log-likelihood, scores
and second derivatives.

## Model fitting

Model fitting is done with a function called `sclr`. It is used in 
the same way as other fitting functions like `lm`.

```{r}
# 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)
```

## Expected protection

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.

```{r}
# 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")])
```

## Protective titres

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`.

```{r}
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)
```