Generic likelihood methods in R Peter Dalgaard Department of - - PowerPoint PPT Presentation

Generic likelihood methods in R

Peter Dalgaard

Department of Biostatistics University of Copenhagen

Wednesday coffee, February 2008

1 / 15

Overview

◮ Pre-history/motivation ◮ Current capabilities in R ◮ Ideas for extensions

2 / 15

Overview

◮ Pre-history/motivation ◮ Current capabilities in R ◮ Ideas for extensions

2 / 15

Overview

◮ Pre-history/motivation ◮ Current capabilities in R ◮ Ideas for extensions

2 / 15

Origin

◮ Ideas about generic likelihood software have been around

for a long time, e.g. in the “Lexical scoping” paper (Gentleman and Ihaka, JCGS 2000)

◮ The concrete occasion was a question on R-help by

Vincent Philion in July 2003

◮ Implementation of mle() in c.2003-2004 + fixups

3 / 15

Origin

◮ Ideas about generic likelihood software have been around

for a long time, e.g. in the “Lexical scoping” paper (Gentleman and Ihaka, JCGS 2000)

◮ The concrete occasion was a question on R-help by

Vincent Philion in July 2003

◮ Implementation of mle() in c.2003-2004 + fixups

3 / 15

Origin

◮ Ideas about generic likelihood software have been around

for a long time, e.g. in the “Lexical scoping” paper (Gentleman and Ihaka, JCGS 2000)

◮ The concrete occasion was a question on R-help by

Vincent Philion in July 2003

◮ Implementation of mle() in c.2003-2004 + fixups

3 / 15

Philion’s problem (lightly edited)

Hello and thank you for your interest in this problem. "real life data" would look like this: x y 28 0.03 21 0.1 11 0.3 15 ... 100 Where X is dose and Y is response. the relation is linear for log(resp) = b log(dose) + intcpt Response for dose 0 is a "control" = Ymax. So, What I want is the dose for 50 percent response. For instance, in example 1: Ymax = 28 (this is also an observation with Poisson error) So I want dose for response = 14 = approx. 0.3

4 / 15

Model

What Philion effectively suggested was

◮ Yi ∼ Pois(λ(xi)) ◮ λ(0) = ymax ◮ log λ(x) = α + β log x for x > 0 ◮ This is standard log-linear and can be fitted with glm() ◮ . . . but it is a strange model!

5 / 15

Model

What Philion effectively suggested was

◮ Yi ∼ Pois(λ(xi)) ◮ λ(0) = ymax ◮ log λ(x) = α + β log x for x > 0 ◮ This is standard log-linear and can be fitted with glm() ◮ . . . but it is a strange model!

5 / 15

Model

What Philion effectively suggested was

◮ Yi ∼ Pois(λ(xi)) ◮ λ(0) = ymax ◮ log λ(x) = α + β log x for x > 0 ◮ This is standard log-linear and can be fitted with glm() ◮ . . . but it is a strange model!

5 / 15

Model

What Philion effectively suggested was

◮ Yi ∼ Pois(λ(xi)) ◮ λ(0) = ymax ◮ log λ(x) = α + β log x for x > 0 ◮ This is standard log-linear and can be fitted with glm() ◮ . . . but it is a strange model!

5 / 15

Model

What Philion effectively suggested was

◮ Yi ∼ Pois(λ(xi)) ◮ λ(0) = ymax ◮ log λ(x) = α + β log x for x > 0 ◮ This is standard log-linear and can be fitted with glm() ◮ . . . but it is a strange model!

5 / 15

Alternative model

◮ Avoid strange behaviour as x → 0 ◮ E.g.: λ(x) = ymax/(1 + x/k) ◮ alias 1/λ(x) = 1/ymax + 1/kymax × x ◮ I.e., this is an inverse-linear Poisson model ◮ glm(y x, poisson("inverse")) +

reparameterization

6 / 15

Alternative model

◮ Avoid strange behaviour as x → 0 ◮ E.g.: λ(x) = ymax/(1 + x/k) ◮ alias 1/λ(x) = 1/ymax + 1/kymax × x ◮ I.e., this is an inverse-linear Poisson model ◮ glm(y x, poisson("inverse")) +

reparameterization

6 / 15

Alternative model

◮ Avoid strange behaviour as x → 0 ◮ E.g.: λ(x) = ymax/(1 + x/k) ◮ alias 1/λ(x) = 1/ymax + 1/kymax × x ◮ I.e., this is an inverse-linear Poisson model ◮ glm(y x, poisson("inverse")) +

reparameterization

6 / 15

Alternative model

◮ Avoid strange behaviour as x → 0 ◮ E.g.: λ(x) = ymax/(1 + x/k) ◮ alias 1/λ(x) = 1/ymax + 1/kymax × x ◮ I.e., this is an inverse-linear Poisson model ◮ glm(y x, poisson("inverse")) +

reparameterization

6 / 15

Alternative model

◮ Avoid strange behaviour as x → 0 ◮ E.g.: λ(x) = ymax/(1 + x/k) ◮ alias 1/λ(x) = 1/ymax + 1/kymax × x ◮ I.e., this is an inverse-linear Poisson model ◮ glm(y x, poisson("inverse")) +

reparameterization

6 / 15

Problem

◮ λ(x) = ymax/(1 + x/k) ◮ For large x this is proportional to x−1 ◮ . . . but not an arbitrary inverse power law (x−β)as in the

• riginal question

◮ Possible fix: λ(x) = ymax/(1 + x/k)β ◮ . . . but this is no longer a generalized linear model

7 / 15

Problem

◮ λ(x) = ymax/(1 + x/k) ◮ For large x this is proportional to x−1 ◮ . . . but not an arbitrary inverse power law (x−β)as in the

• riginal question

◮ Possible fix: λ(x) = ymax/(1 + x/k)β ◮ . . . but this is no longer a generalized linear model

7 / 15

Problem

◮ λ(x) = ymax/(1 + x/k) ◮ For large x this is proportional to x−1 ◮ . . . but not an arbitrary inverse power law (x−β)as in the

• riginal question

◮ Possible fix: λ(x) = ymax/(1 + x/k)β ◮ . . . but this is no longer a generalized linear model

7 / 15

Problem

◮ λ(x) = ymax/(1 + x/k) ◮ For large x this is proportional to x−1 ◮ . . . but not an arbitrary inverse power law (x−β)as in the

• riginal question

◮ Possible fix: λ(x) = ymax/(1 + x/k)β ◮ . . . but this is no longer a generalized linear model

7 / 15

Problem

◮ λ(x) = ymax/(1 + x/k) ◮ For large x this is proportional to x−1 ◮ . . . but not an arbitrary inverse power law (x−β)as in the

• riginal question

◮ Possible fix: λ(x) = ymax/(1 + x/k)β ◮ . . . but this is no longer a generalized linear model

7 / 15

This is silly!

◮ Why do we let the existence of standard model families

interfere with a sensible choice of model?

◮ When all you have is a hammer. . . ◮ But we have many other tools! ◮ E.g. quite powerful optimizers like optim ◮ Why not just write out the likelihood and maximize it?

8 / 15

This is silly!

◮ Why do we let the existence of standard model families

interfere with a sensible choice of model?

◮ When all you have is a hammer. . . ◮ But we have many other tools! ◮ E.g. quite powerful optimizers like optim ◮ Why not just write out the likelihood and maximize it?

8 / 15

This is silly!

◮ Why do we let the existence of standard model families

interfere with a sensible choice of model?

◮ When all you have is a hammer. . . ◮ But we have many other tools! ◮ E.g. quite powerful optimizers like optim ◮ Why not just write out the likelihood and maximize it?

8 / 15

This is silly!

◮ Why do we let the existence of standard model families

interfere with a sensible choice of model?

◮ When all you have is a hammer. . . ◮ But we have many other tools! ◮ E.g. quite powerful optimizers like optim ◮ Why not just write out the likelihood and maximize it?

8 / 15

This is silly!

◮ Why do we let the existence of standard model families

interfere with a sensible choice of model?

◮ When all you have is a hammer. . . ◮ But we have many other tools! ◮ E.g. quite powerful optimizers like optim ◮ Why not just write out the likelihood and maximize it?

8 / 15

The mle function

◮ Basically just a wrapper for optim("BFGS" method by

default.

◮ Supply − log L and the routine does the rest (hopefully). ◮ Nice methods for summary display, likelihood profiling

( √ LRT, possibly signed), and approximate confidence intervals

9 / 15

The mle function

◮ Basically just a wrapper for optim("BFGS" method by

default.

◮ Supply − log L and the routine does the rest (hopefully). ◮ Nice methods for summary display, likelihood profiling

( √ LRT, possibly signed), and approximate confidence intervals

9 / 15

The mle function

◮ Basically just a wrapper for optim("BFGS" method by

default.

◮ Supply − log L and the routine does the rest (hopefully). ◮ Nice methods for summary display, likelihood profiling

( √ LRT, possibly signed), and approximate confidence intervals

9 / 15

How to use it

> library(stats4) > mll <- function(ymax, k, la) + with(e1, +

• sum(dpois(response,

+ ymax / (1 + dose/k)^exp(la), + log=TRUE))) > fit <- mle(mll, start=list(ymax=28, k=.3, la=0))

10 / 15

Output

> summary(fit) Maximum likelihood estimation Call: mle(minuslogl = mll, start = list(ymax = 28, k = 0.3, la = 0)) Coefficients: Estimate Std. Error ymax 24.4735208 4.7498215 k 0.1884905 0.2735927 la

• 0.2285072

0.4696129

• 2 log L: 33.80755

11 / 15

Profiling

> par(mfrow=c(3,1)) > plot(profile(fit, del=.05))

15 20 25 30 35 40 0.0 1.0 2.0 ymax z 0.0 0.5 1.0 1.5 0.0 1.0 2.0 k z −1.0 −0.5 0.0 0.5 1.0 0.0 1.0 2.0 la z

12 / 15

Confidence intervals

> confint(profile(fit, del=.05, maxsteps=500)) 2.5 % 97.5 % ymax 17.35158789 36.428935 k 0.01226197 3.562976 la

• 0.94484692

1.131715

13 / 15

Design issues

◮ Implemented using the S4 class system. Mainly to get my

feet wet.

◮ The argument to mle is a function of the arguments only ◮ This is the cleanest design (and, sort of, in accordance

with the strict likelihood principle)

◮ Some people have been requesting a more traditional

interface with a data= argument, etc. but that is less general (e.g. likelihoods can depend on multiple data sets), so I have been resisting that, but that didn’t stop Ben Bolker (bbmle package). . .

◮ Ben has also added several apparently useful features,

though

14 / 15

Design issues

◮ Implemented using the S4 class system. Mainly to get my

feet wet.

◮ The argument to mle is a function of the arguments only ◮ This is the cleanest design (and, sort of, in accordance

with the strict likelihood principle)

◮ Some people have been requesting a more traditional

interface with a data= argument, etc. but that is less general (e.g. likelihoods can depend on multiple data sets), so I have been resisting that, but that didn’t stop Ben Bolker (bbmle package). . .

◮ Ben has also added several apparently useful features,

though

14 / 15

Design issues

◮ Implemented using the S4 class system. Mainly to get my

feet wet.

◮ The argument to mle is a function of the arguments only ◮ This is the cleanest design (and, sort of, in accordance

with the strict likelihood principle)

◮ Some people have been requesting a more traditional

interface with a data= argument, etc. but that is less general (e.g. likelihoods can depend on multiple data sets), so I have been resisting that, but that didn’t stop Ben Bolker (bbmle package). . .

◮ Ben has also added several apparently useful features,

though

14 / 15

Design issues

◮ Implemented using the S4 class system. Mainly to get my

feet wet.

◮ The argument to mle is a function of the arguments only ◮ This is the cleanest design (and, sort of, in accordance

with the strict likelihood principle)

◮ Some people have been requesting a more traditional

interface with a data= argument, etc. but that is less general (e.g. likelihoods can depend on multiple data sets), so I have been resisting that, but that didn’t stop Ben Bolker (bbmle package). . .

◮ Ben has also added several apparently useful features,

though

14 / 15

Design issues

◮ Implemented using the S4 class system. Mainly to get my

feet wet.

◮ The argument to mle is a function of the arguments only ◮ This is the cleanest design (and, sort of, in accordance

with the strict likelihood principle)

◮ Some people have been requesting a more traditional

interface with a data= argument, etc. but that is less general (e.g. likelihoods can depend on multiple data sets), so I have been resisting that, but that didn’t stop Ben Bolker (bbmle package). . .

◮ Ben has also added several apparently useful features,

though

14 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15

Ideas/TODO list

◮ Likelihood “Meccano” (or -algebra)

◮ Simplified specification in special cases, e.g.

mll <- PoissonLike(response ~ ymax / (1 + dose/k) ^ exp(la), ...)

◮ Combination of multiple independent likelihoods (with

common parameters)

◮ Self-starting models as in nls ◮ Multilevel models (AGQ, etc.) ◮ Improved approximations (cf. hoa package) ◮ Fixing the documentation (!!!) ◮ Multivariate profiling

15 / 15