Getting Started with likelihood.model
This tutorial walks you through the fundamentals of building likelihood models for statistical inference.
Installation
Install the development version from GitHub:
if (!require(devtools)) {
install.packages("devtools")
}
devtools::install_github("queelius/likelihood.model")
You’ll also want the companion package for MLE:
devtools::install_github("queelius/algebraic.mle")
Tutorial 1: Your First Likelihood Model
Let’s fit a Weibull distribution to some simulated survival data.
Step 1: Generate Sample Data
library(likelihood.model)
set.seed(42)
# True parameters
true_shape <- 2
true_scale <- 10
# Generate survival times
n <- 100
times <- rweibull(n, shape = true_shape, scale = true_scale)
# Simulate right-censoring at time = 15
censor_time <- 15
observed_times <- pmin(times, censor_time)
event_indicator <- as.numeric(times <= censor_time)
# Create data frame
data <- data.frame(
time = observed_times,
status = event_indicator
)
Step 2: Specify the Likelihood Model
# Create a right-censored Weibull model
model <- likelihood_contr_model(
obs_type = "right_censored",
distribution = weibull_distribution(),
time_col = "time",
status_col = "status"
)
Step 3: Compute the Log-Likelihood
# Evaluate at a specific parameter vector
theta <- c(shape = 1.5, scale = 8)
ll <- loglik(model, theta, data)
print(ll)
Step 4: Maximum Likelihood Estimation
library(algebraic.mle)
# Fit the model
fit <- mle(model, data = data, start = c(shape = 1, scale = 5))
# View results
summary(fit)
Tutorial 2: Mixed Observation Types
Real datasets often contain multiple observation types. Here’s how to handle them.
Combining Exact and Censored Observations
# Some observations are exact, others are right-censored
data_mixed <- data.frame(
time = c(3.2, 5.1, 7.8, 10.0, 12.5, 15.0, 15.0),
status = c(1, 1, 1, 1, 1, 0, 0), # 0 = censored
obs_type = c("exact", "exact", "exact", "exact", "exact",
"right_censored", "right_censored")
)
# The model automatically handles both types based on status
model_mixed <- likelihood_contr_model(
obs_type = "right_censored", # Framework handles both
distribution = weibull_distribution()
)
Tutorial 3: Custom Distributions
You can define custom distributions for specialized applications.
# Define a custom distribution
my_distribution <- function() {
structure(
list(
name = "my_dist",
n_params = 2
),
class = "distribution"
)
}
# Implement required methods
pdf.my_dist <- function(dist, x, theta) {
# Your density function
}
cdf.my_dist <- function(dist, x, theta) {
# Your CDF
}
sf.my_dist <- function(dist, x, theta) {
# Survival function: 1 - CDF
}
Next Steps
- Explore the Examples for real-world applications
- Read the API Documentation for complete function reference
- Check out related blog posts for in-depth discussions