R Under development (unstable) (2020-03-24 r78039) -- "Unsuffered Consequences" Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > pkgname <- "stats4" > source(file.path(R.home("share"), "R", "examples-header.R")) > options(warn = 1) > library('stats4') > > base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') > base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv') > cleanEx() > nameEx("mle") > ### * mle > > flush(stderr()); flush(stdout()) > > ### Name: mle > ### Title: Maximum Likelihood Estimation > ### Aliases: mle > ### Keywords: models > > ### ** Examples > > ## Avoid printing to unwarranted accuracy > od <- options(digits = 5) > > ## Simulated EC50 experiment with count data > x <- 0:10 > y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) > > ## Easy one-dimensional MLE: > nLL <- function(lambda) -sum(stats::dpois(y, lambda, log = TRUE)) > fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y)) > > ## sanity check --- notice that "nobs" must be input > ## (not guaranteed to be meaningful for any likelihood) > stopifnot(nobs(fit0) == length(y)) > > > # For 1D, this is preferable: > fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y), + method = "Brent", lower = 1, upper = 20) > > ## This needs a constrained parameter space: most methods will accept NA > ll <- function(ymax = 15, xhalf = 6) { + if(ymax > 0 && xhalf > 0) + -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)) + else NA + } > (fit <- mle(ll, nobs = length(y))) Call: mle(minuslogl = ll, nobs = length(y)) Coefficients: ymax xhalf 24.9931 3.0571 > mle(ll, fixed = list(xhalf = 6)) Call: mle(minuslogl = ll, fixed = list(xhalf = 6)) Coefficients: ymax xhalf 19.288 6.000 > > ## Alternative using bounds on optimization > ll2 <- function(ymax = 15, xhalf = 6) + -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)) > mle(ll2, lower = rep(0, 2)) Call: mle(minuslogl = ll2, lower = rep(0, 2)) Coefficients: ymax xhalf 24.9994 3.0558 > > AIC(fit) [1] 61.208 > BIC(fit) [1] 62.004 > > summary(fit) Maximum likelihood estimation Call: mle(minuslogl = ll, nobs = length(y)) Coefficients: Estimate Std. Error ymax 24.9931 4.2244 xhalf 3.0571 1.0348 -2 log L: 57.208 > logLik(fit) 'log Lik.' -28.604 (df=2) > vcov(fit) ymax xhalf ymax 17.8459 -3.7206 xhalf -3.7206 1.0708 > plot(profile(fit), absVal = FALSE) > confint(fit) Profiling... 2.5 % 97.5 % ymax 17.8845 34.6194 xhalf 1.6616 6.4792 > > ## Use bounded optimization > ## The lower bounds are really > 0, > ## but we use >=0 to stress-test profiling > (fit2 <- mle(ll2, lower = c(0, 0))) Call: mle(minuslogl = ll2, lower = c(0, 0)) Coefficients: ymax xhalf 24.9994 3.0558 > plot(profile(fit2), absVal = FALSE) > > ## A better parametrization: > ll3 <- function(lymax = log(15), lxhalf = log(6)) + -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)) > (fit3 <- mle(ll3)) Call: mle(minuslogl = ll3) Coefficients: lymax lxhalf 3.2189 1.1170 > plot(profile(fit3), absVal = FALSE) > exp(confint(fit3)) Profiling... 2.5 % 97.5 % lymax 17.8815 34.6186 lxhalf 1.6615 6.4794 > > # Regression tests for bounded cases (this was broken in R 3.x) > fit4 <- mle(ll, lower = c(0, 4)) # has max on boundary > confint(fit4) Profiling... 2.5 % 97.5 % ymax 17.446 26.5081 xhalf NA 6.9109 > > ## direct check that fixed= and constraints work together > mle(ll, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary Call: mle(minuslogl = ll, fixed = list(ymax = 23), lower = c(0, 4)) Coefficients: ymax xhalf 23 4 > > ## Linear regression using MLE > x <- 1:10 > y <- c(0.48, 2.24, 2.22, 5.15, 4.64, 5.53, 7, 8.8, 7.67, 9.23) > > LM_mll <- function(formula, data = environment(formula)) + { + y <- model.response(model.frame(formula, data)) + X <- model.matrix(formula, data) + b0 <- numeric(NCOL(X)) + names(b0) <- colnames(X) + function(b=b0, sigma=1) + -sum(dnorm(y, X %*% b, sigma, log=TRUE)) + } > > mll <- LM_mll(y ~ x) > > summary(lm(y~x)) # for comparison -- notice variance bias in MLE Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -0.937 -0.500 -0.211 0.278 1.273 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0927 0.5376 0.17 0.87 x 0.9461 0.0866 10.92 4.4e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.787 on 8 degrees of freedom Multiple R-squared: 0.937, Adjusted R-squared: 0.929 F-statistic: 119 on 1 and 8 DF, p-value: 4.39e-06 > summary(mle(mll, lower=c(-Inf,-Inf, 0.01))) Maximum likelihood estimation Call: mle(minuslogl = mll, lower = c(-Inf, -Inf, 0.01)) Coefficients: Estimate Std. Error b.(Intercept) 0.092667 0.480869 b.x 0.946061 0.077499 sigma 0.703919 0.157400 -2 log L: 21.357 > summary(mle(mll, lower=list(sigma = 0.01))) # alternative specification Maximum likelihood estimation Call: mle(minuslogl = mll, lower = list(sigma = 0.01)) Coefficients: Estimate Std. Error b.(Intercept) 0.092667 0.480869 b.x 0.946061 0.077499 sigma 0.703919 0.157400 -2 log L: 21.357 > > confint(mle(mll, lower=list(sigma = 0.01))) Profiling... 2.5 % 97.5 % b.(Intercept) -0.94831 1.1336 b.x 0.77829 1.1138 sigma 0.48017 1.1755 > plot(profile(mle(mll, lower=list(sigma = 0.01)))) > > Binom_mll <- function(x, n) + { + force(x); force(n) ## beware lazy evaluation + function(p=.5) -dbinom(x, n, p, log=TRUE) + } > > ## Likelihood functions for different x. > ## This code goes wrong, if force(x) is not used in Binom_mll: > > curve(Binom_mll(0, 10)(p), xname="p", ylim=c(0, 10)) > mll_list <- list(10) > for (x in 1:10) + mll_list[[x]] <- Binom_mll(x, 10) > for (mll in mll_list) + curve(mll(p), xname="p", add=TRUE) > > mll <- Binom_mll(4,10) > mle(mll, lower = 1e-16, upper = 1-1e-16) # limits must be inside (0,1) Call: mle(minuslogl = mll, lower = 1e-16, upper = 1 - 1e-16) Coefficients: p 0.4 > > ## Boundary case: This works, but fails if limits are set closer to 0 and 1 > mll <- Binom_mll(0, 10) > mle(mll, lower=.005, upper=.995) Call: mle(minuslogl = mll, lower = 0.005, upper = 0.995) Coefficients: p 0.005 > > ## Not run: > ##D ## We can use limits closer to the boundaries if we use the > ##D ## drop-in replacement optimr() from the optimx package. > ##D > ##D mle(mll, lower = 1e-16, upper = 1-1e-16, optim=optimx::optimr) > ## End(Not run) > > > options(od) > > > > cleanEx() > nameEx("update-methods") > ### * update-methods > > flush(stderr()); flush(stdout()) > > ### Name: update-methods > ### Title: Methods for Function 'update' in Package 'stats4' > ### Aliases: update-methods update,ANY-method update,mle-method > ### Keywords: methods > > ### ** Examples > > x <- 0:10 > y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) > ll <- function(ymax = 15, xhalf = 6) + -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)) > fit <- mle(ll) Warning in stats::dpois(y, lambda = ymax/(1 + x/xhalf), log = TRUE) : NaNs produced > ## note the recorded call contains ..1, a problem with S4 dispatch > update(fit, fixed = list(xhalf = 3)) Call: mle(minuslogl = ll, fixed = ..1) Coefficients: ymax xhalf 25.19609 3.00000 > > > > ### *