melt
Overview
melt provides a unified framework for data analysis with empirical likelihood methods. A collection of functions is available to perform multiple empirical likelihood tests and construct confidence intervals for various models in ‘R’. melt offers an easy-to-use interface and flexibility in specifying hypotheses and calibration methods, extending the framework to simultaneous inferences. The core computational routines are implemented with the ‘Eigen’ ‘C++’ library and ‘RcppEigen’ interface, with ‘OpenMP’ for parallel computation. Details of the testing procedures are provided in Kim, MacEachern, and Peruggia (2023). The package has a companion paper by Kim, MacEachern, and Peruggia (2024). This work was supported by the U.S. National Science Foundation under Grants No. SES-1921523 and DMS-2015552.
Installation
You can install the latest stable release of melt from CRAN.
install.packages("melt")You can install the development version of melt from GitHub or R-universe.
# install.packages("pak")
pak::pak("ropensci/melt")install.packages("melt", repos = "https://ropensci.r-universe.dev")Main functions
melt provides an intuitive API for performing the most common data analysis tasks:
el_mean()computes empirical likelihood for the mean.el_lm()fits a linear model with empirical likelihood.el_glm()fits a generalized linear model with empirical likelihood.confint()computes confidence intervals for model parameters.confreg()computes confidence region for model parameters.elt()tests a hypothesis with various calibration options.elmt()performs multiple testing simultaneously.
Usage
library(melt)
set.seed(971112)
## Test for the mean
data("precip")
(fit <- el_mean(precip, par = 30))
#>
#> Empirical Likelihood
#>
#> Model: mean
#>
#> Maximum EL estimates:
#> [1] 34.89
#>
#> Chisq: 8.285, df: 1, Pr(>Chisq): 0.003998
#> EL evaluation: converged
## Adjusted empirical likelihood calibration
elt(fit, rhs = 30, calibrate = "ael")
#>
#> Empirical Likelihood Test
#>
#> Hypothesis:
#> par = 30
#>
#> Significance level: 0.05, Calibration: Adjusted EL
#>
#> Statistic: 7.744, Critical value: 3.841
#> p-value: 0.005389
#> EL evaluation: converged
## Bootstrap calibration
elt(fit, rhs = 30, calibrate = "boot")
#>
#> Empirical Likelihood Test
#>
#> Hypothesis:
#> par = 30
#>
#> Significance level: 0.05, Calibration: Bootstrap
#>
#> Statistic: 8.285, Critical value: 3.84
#> p-value: 0.0041
#> EL evaluation: converged
## F calibration
elt(fit, rhs = 30, calibrate = "f")
#>
#> Empirical Likelihood Test
#>
#> Hypothesis:
#> par = 30
#>
#> Significance level: 0.05, Calibration: F
#>
#> Statistic: 8.285, Critical value: 3.98
#> p-value: 0.005318
#> EL evaluation: converged
## Linear model
data("mtcars")
fit_lm <- el_lm(mpg ~ disp + hp + wt + qsec, data = mtcars)
summary(fit_lm)
#>
#> Empirical Likelihood
#>
#> Model: lm
#>
#> Call:
#> el_lm(formula = mpg ~ disp + hp + wt + qsec, data = mtcars)
#>
#> Number of observations: 32
#> Number of parameters: 5
#>
#> Parameter values under the null hypothesis:
#> (Intercept) disp hp wt qsec
#> 29.04 0.00 0.00 0.00 0.00
#>
#> Lagrange multipliers:
#> [1] -260.167 -2.365 1.324 -59.781 25.175
#>
#> Maximum EL estimates:
#> (Intercept) disp hp wt qsec
#> 27.329638 0.002666 -0.018666 -4.609123 0.544160
#>
#> logL: -327.6 , logLR: -216.7
#> Chisq: 433.4, df: 4, Pr(>Chisq): < 2.2e-16
#> Constrained EL: converged
#>
#> Coefficients:
#> Estimate Chisq Pr(>Chisq)
#> (Intercept) 27.329638 443.208 < 2e-16 ***
#> disp 0.002666 0.365 0.54575
#> hp -0.018666 10.730 0.00105 **
#> wt -4.609123 439.232 < 2e-16 ***
#> qsec 0.544160 440.583 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cr <- confreg(fit_lm, parm = c("disp", "hp"), npoints = 200)
plot(cr)data("clothianidin")
fit2_lm <- el_lm(clo ~ -1 + trt, data = clothianidin)
summary(fit2_lm)
#>
#> Empirical Likelihood
#>
#> Model: lm
#>
#> Call:
#> el_lm(formula = clo ~ -1 + trt, data = clothianidin)
#>
#> Number of observations: 102
#> Number of parameters: 4
#>
#> Parameter values under the null hypothesis:
#> trtNaked trtFungicide trtLow trtHigh
#> 0 0 0 0
#>
#> Lagrange multipliers:
#> [1] -4.116e+06 -7.329e-01 -1.751e+00 -1.418e-01
#>
#> Maximum EL estimates:
#> trtNaked trtFungicide trtLow trtHigh
#> -4.479 -3.427 -2.800 -1.307
#>
#> logL: -918.9 , logLR: -447.2
#> Chisq: 894.4, df: 4, Pr(>Chisq): < 2.2e-16
#> EL evaluation: maximum iterations reached
#>
#> Coefficients:
#> Estimate Chisq Pr(>Chisq)
#> trtNaked -4.479 411.072 < 2e-16 ***
#> trtFungicide -3.427 59.486 1.23e-14 ***
#> trtLow -2.800 62.955 2.11e-15 ***
#> trtHigh -1.307 4.653 0.031 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(fit2_lm)
#> lower upper
#> trtNaked -5.002118 -3.9198229
#> trtFungicide -4.109816 -2.6069870
#> trtLow -3.681837 -1.9031795
#> trtHigh -2.499165 -0.1157222
## Generalized linear model
data("thiamethoxam")
fit_glm <- el_glm(visit ~ log(mass) + fruit + foliage + var + trt,
family = quasipoisson(link = "log"), data = thiamethoxam,
control = el_control(maxit = 100, tol = 1e-08, nthreads = 4)
)
summary(fit_glm)
#>
#> Empirical Likelihood
#>
#> Model: glm (quasipoisson family with log link)
#>
#> Call:
#> el_glm(formula = visit ~ log(mass) + fruit + foliage + var +
#> trt, family = quasipoisson(link = "log"), data = thiamethoxam,
#> control = el_control(maxit = 100, tol = 1e-08, nthreads = 4))
#>
#> Number of observations: 165
#> Number of parameters: 8
#>
#> Parameter values under the null hypothesis:
#> (Intercept) log(mass) fruit foliage varGZ trtSpray
#> -0.1098 0.0000 0.0000 0.0000 0.0000 0.0000
#> trtFurrow trtSeed phi
#> 0.0000 0.0000 1.4623
#>
#> Lagrange multipliers:
#> [1] 1319.19 210.54 -12.99 -24069.07 -318.90 -189.14 -53.35
#> [8] 262.32 -170.21
#>
#> Maximum EL estimates:
#> (Intercept) log(mass) fruit foliage varGZ trtSpray
#> -0.10977 0.24750 0.04654 -19.40632 -0.25760 0.06724
#> trtFurrow trtSeed
#> -0.03634 0.34790
#>
#> logL: -2272 , logLR: -1429
#> Chisq: 2859, df: 7, Pr(>Chisq): < 2.2e-16
#> Constrained EL: initialization failed
#>
#> Coefficients:
#> Estimate Chisq Pr(>Chisq)
#> (Intercept) -0.10977 0.090 0.764
#> log(mass) 0.24750 425.859 < 2e-16 ***
#> fruit 0.04654 29.024 7.15e-08 ***
#> foliage -19.40632 65.181 6.83e-16 ***
#> varGZ -0.25760 17.308 3.18e-05 ***
#> trtSpray 0.06724 0.860 0.354
#> trtFurrow -0.03634 0.217 0.641
#> trtSeed 0.34790 19.271 1.13e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Dispersion for quasipoisson family: 1.462288
## Test of no treatment effect
contrast <- c(
"trtNaked - trtFungicide", "trtFungicide - trtLow", "trtLow - trtHigh"
)
elt(fit2_lm, lhs = contrast)
#>
#> Empirical Likelihood Test
#>
#> Hypothesis:
#> trtNaked - trtFungicide = 0
#> trtFungicide - trtLow = 0
#> trtLow - trtHigh = 0
#>
#> Significance level: 0.05, Calibration: Chi-square
#>
#> Statistic: 26.6, Critical value: 7.815
#> p-value: 7.148e-06
#> Constrained EL: converged
## Multiple testing
contrast2 <- rbind(
c(0, 0, 0, 0, 0, 1, 0, 0),
c(0, 0, 0, 0, 0, 0, 1, 0),
c(0, 0, 0, 0, 0, 0, 0, 1)
)
elmt(fit_glm, lhs = contrast2)
#>
#> Empirical Likelihood Multiple Tests
#>
#> Overall significance level: 0.05
#>
#> Calibration: Multivariate chi-square
#>
#> Hypotheses:
#> Estimate Chisq Df
#> trtSpray = 0 0.06724 0.860 1
#> trtFurrow = 0 -0.03634 0.217 1
#> trtSeed = 0 0.34790 19.271 1Please note that this package is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.