# NOT RUN {
library(lodr)
## Using example dataset provided in lodr package: lod_data_ex
## 3 covariates: x1, x2, x3 with x2 and x3 subject to a lower limit of
## detection of 0
# Replace values marked as under limit of detection using 0 with NA,
# add column of ones for intercept
lod_data_with_int <-
as.matrix(cbind("Intercept"=rep(1, dim(lod_data_ex)[1]), lod_data_ex))
lod_data_ex_edit <-
data.frame(apply(lod_data_with_int, MARGIN = 2, FUN=function(x){ifelse(x==0, NA, x)}))
# Fit linear model to dataset with only subjects without covariates under
# limit of detection to get initial estimate for the regression parameters.
beta_inital_est <- coef(lm(y~x1+x2+x3, data=lod_data_ex_edit))
# Get initial estimates of mean vector and covariance matrix for covariates and variance of outcome,
# again using data from subjects without covariates under limit of detection
mean_x_inital <- colMeans(lod_data_ex_edit[,c(-1,-2)], na.rm = TRUE)
sigma_x_inital <- cov(lod_data_ex_edit[,c(-1,-2)], use="pairwise.complete.obs")
sigma_2_y_inital <- sigma(lm(y~x1+x2+x3, data=lod_data_ex_edit))^2
# Fit model, report regression parameter estimates from last iteration
LOD_matrix <- cbind(c(NA, NA, -100, -100), c(NA, NA, 0, 0))
## no_of_samples set to 100 for computational speed/illustration purposes only.
## At least 250 is recommended.
fit_object <-
LOD_fit(y_data=lod_data_ex_edit[,2],
x_data=as.matrix(lod_data_ex_edit[,-2]),
mean_x_preds=mean_x_inital, beta=beta_inital_est, sigma_2_y=sigma_2_y_inital,
sigma_x_preds=sigma_x_inital, no_of_samples=100,
threshold=0.001, max_iterations=100, LOD_u_l=LOD_matrix)
fit_object$beta_estimate_last_iteration
# }
Run the code above in your browser using DataLab