Beyond GLM: The potential for a generic likelihood toolbox Peter - - PowerPoint PPT Presentation

Beyond GLM: The potential for a generic likelihood toolbox

Peter Dalgaard

Department of Biostatistics University of Copenhagen

Royal Statistical Society, Nottingham, September 2008

1 / 14

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Introduction

◮ Developments in statistics are driven by its tools ◮ Tools that we have available are much stronger than what

we currently use them for

◮ It is feasible to use features of the R language to work with

models at a level that transcends the current practice

◮ (However, this is all still at a very preliminary stage) ◮ Overview

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

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Generalized Linear Models

◮ 3 dozen years ago generalized linear models were born ◮ This extended the capabilities of regression analysis and

brought major modes of analysis into the same framework

◮ Multiple regression ◮ Logistic regression ◮ Poisson (rate) regression ◮ (Cox regression)

◮ Encapsulated in GLIM and Genstat software (c.1974)

3 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Some observations about GLMs

◮ Reuse of prior technology/ideas

◮ Weighted least squares ◮ Analysis of variance/deviance ◮ Wilkinson-Rogers formulas

◮ Innovation

◮ Use of W-R formulas by software (Genstats "TREAT"

directive)

◮ Focusing

◮ Esp. on “What can we do with linear predictors”

◮ Delimiting

◮ Tendency to “forget” other options ◮ Partially false distinctions between models that “can be

fitted with standard software” and “difficult models”

4 / 14

Postulate

◮ Useful as they are, there are cases where GLM modelling

actually inhibits proper analysis

◮ We should become better at working with non-standard

models

◮ As a starting point, we could consider

5 / 14

Postulate

◮ Useful as they are, there are cases where GLM modelling

actually inhibits proper analysis

◮ We should become better at working with non-standard

models

◮ As a starting point, we could consider

5 / 14

Postulate

◮ Useful as they are, there are cases where GLM modelling

actually inhibits proper analysis

◮ We should become better at working with non-standard

models

◮ As a starting point, we could consider

5 / 14

Illustrative example: Philion (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

5 / 14

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!

6 / 14

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!

6 / 14

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!

6 / 14

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!

6 / 14

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!

6 / 14

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")) + reparametrization

7 / 14

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")) + reparametrization

7 / 14

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")) + reparametrization

7 / 14

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")) + reparametrization

7 / 14

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")) + reparametrization

7 / 14

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

8 / 14

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

8 / 14

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

8 / 14

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

8 / 14

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

8 / 14

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?

9 / 14

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?

9 / 14

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?

9 / 14

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?

9 / 14

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?

9 / 14

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

10 / 14

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

10 / 14

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

10 / 14

How to use it

> library(stats4) > library(ISwR); e1 <- subset(philion,experiment==1) > 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 / 14

Output

> summary(fit) Maximum likelihood estimation Call: mle(minuslogl = mll, start = list(ymax = 28, k = 0.3, la = 0)) Coefficients: Estimate Std. Error ymax 26.1709859 5.1073174 k 0.1935226 0.1551403 la

• 0.2548393

0.2514938

• 2 log L: 33.14818

10 / 14

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.2 0.4 0.6 0.8 1.0 0.0 1.0 2.0 k z −0.5 0.0 0.5 0.0 1.0 2.0 la z

10 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Further advantages to working with likelihoods

◮ Penalties / Prior information ◮ Truncation ◮ Censoring ◮ Mixtures ◮ Combined experiments ◮ Misclassification ◮ Random parameters ◮ Modified likelihoods

(of course, not all likelihoods are equally tractable. . . )

11 / 14

Model specification

◮ Direct specification of minus log likelihood is somewhat

error-prone

◮ Automatic generation of likelihoods would be an attractive

• ption

◮ R allows computation on the language ◮ A preliminary sketch: How to generate the likelihood from a

generalized linear model

12 / 14

Model specification

◮ Direct specification of minus log likelihood is somewhat

error-prone

◮ Automatic generation of likelihoods would be an attractive

• ption

◮ R allows computation on the language ◮ A preliminary sketch: How to generate the likelihood from a

generalized linear model

12 / 14

Model specification

◮ Direct specification of minus log likelihood is somewhat

error-prone

◮ Automatic generation of likelihoods would be an attractive

• ption

◮ R allows computation on the language ◮ A preliminary sketch: How to generate the likelihood from a

generalized linear model

12 / 14

Model specification

◮ Direct specification of minus log likelihood is somewhat

error-prone

◮ Automatic generation of likelihoods would be an attractive

• ption

◮ R allows computation on the language ◮ A preliminary sketch: How to generate the likelihood from a

generalized linear model

12 / 14

A direct approach

## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18,17,15,20,10,20,25,13,12)

• utcome <- gl(3,1,9)

treatment <- gl(3,3) glm.D93 <- glm(counts ~ outcome + treatment, family=poisson()) summary(glm.D93) ## Same thing, using mle() X <- model.matrix(counts ~ outcome + treatment) mll <- function(a=0, o2=0, o3=0, t2=0, t3=0){ beta <- c(a,o2,o3,t2,t3) lp <- X %*% beta

• sum(dpois(counts, exp(lp), log=TRUE))

} mle.D93 <- mle(mll) summary(mle.D93)

12 / 14

More elaborate version

formula <- counts ~ outcome + treatment family.mll <- function(obs, pred)

• sum(dpois(obs, exp(pred), log=TRUE))

# generic code # ------------ X <- model.matrix(formula) fr <- model.frame(formula) Y <- model.response(fr) L <- rep(list(0),ncol(X)); names(L) <- colnames(X) mll <- function(){ cl <- match.call() cl[[1]] <- ‘c‘ beta <- eval.parent(cl) lp <- X %*% beta family.mll(Y, lp) } formals(mll) <- L #------------ summary(mle(mll))

12 / 14

Likelihood algebra

◮ Given that we have ways of constructing likelihoods, we

might also want to have ways of modifiying them:

◮ reparametrization ◮ joint likelihoods ◮ parameter constraints ◮ Longer term: generic handling of multilevel models (AGQ,

etc.)

◮ Even the simpler items above lead to some interesting

issues when you try to implement them generically

13 / 14

Likelihood algebra

◮ Given that we have ways of constructing likelihoods, we

might also want to have ways of modifiying them:

◮ reparametrization ◮ joint likelihoods ◮ parameter constraints ◮ Longer term: generic handling of multilevel models (AGQ,

etc.)

◮ Even the simpler items above lead to some interesting

issues when you try to implement them generically

13 / 14

Likelihood algebra

◮ Given that we have ways of constructing likelihoods, we

might also want to have ways of modifiying them:

◮ reparametrization ◮ joint likelihoods ◮ parameter constraints ◮ Longer term: generic handling of multilevel models (AGQ,

etc.)

◮ Even the simpler items above lead to some interesting

issues when you try to implement them generically

13 / 14

Likelihood algebra

◮ Given that we have ways of constructing likelihoods, we

might also want to have ways of modifiying them:

◮ reparametrization ◮ joint likelihoods ◮ parameter constraints ◮ Longer term: generic handling of multilevel models (AGQ,

etc.)

◮ Even the simpler items above lead to some interesting

issues when you try to implement them generically

13 / 14

Likelihood algebra

◮ Given that we have ways of constructing likelihoods, we

might also want to have ways of modifiying them:

◮ reparametrization ◮ joint likelihoods ◮ parameter constraints ◮ Longer term: generic handling of multilevel models (AGQ,

etc.)

◮ Even the simpler items above lead to some interesting

issues when you try to implement them generically

13 / 14

Likelihood algebra

◮ Given that we have ways of constructing likelihoods, we

might also want to have ways of modifiying them:

◮ reparametrization ◮ joint likelihoods ◮ parameter constraints ◮ Longer term: generic handling of multilevel models (AGQ,

etc.)

◮ Even the simpler items above lead to some interesting

issues when you try to implement them generically

13 / 14

Example: Reparametrization

Suppose you want to define

Reparam(f, list(b=quote(exp(d))), list(d=0))

so that given

f <- function(a=0,b=1,c=0){....}

the return value is

function(a=0,c=0,d=0) f(a, exp(d), c)

How do you actually do that?

13 / 14

Example: Reparametrization

Suppose you want to define

Reparam(f, list(b=quote(exp(d))), list(d=0))

so that given

f <- function(a=0,b=1,c=0){....}

the return value is

function(a=0,c=0,d=0) f(a, exp(d), c)

How do you actually do that?

13 / 14

Example: Reparametrization

Suppose you want to define

Reparam(f, list(b=quote(exp(d))), list(d=0))

so that given

f <- function(a=0,b=1,c=0){....}

the return value is

function(a=0,c=0,d=0) f(a, exp(d), c)

How do you actually do that?

13 / 14

Example: Reparametrization

Suppose you want to define

Reparam(f, list(b=quote(exp(d))), list(d=0))

so that given

f <- function(a=0,b=1,c=0){....}

the return value is

function(a=0,c=0,d=0) f(a, exp(d), c)

How do you actually do that?

13 / 14

Example: Reparametrization

Here’s one way

Reparam <- function(f, redef, xpars) { FUN <- substitute(f) l <- formals(f) redefined <- match(names(redef), names(l)) newformals <- c(l[-redefined], xpars) args <- lapply(names(l), as.name) names(args) <- names(l) args[redefined] <- redef cl <- as.call(c(FUN,args)) ff <- f formals(ff) <- newformals body(ff) <- cl ff }

13 / 14

Summary

◮ Standard modelling tools are great, but can be limiting ◮ It is for multiple reasons desirable to work with models

through their likelihood functions and develop generic methods for analyzing and manipulating them

◮ The structure of R as a functional programming language

makes this feasible in a fairly straightforward way

14 / 14

Summary

◮ Standard modelling tools are great, but can be limiting ◮ It is for multiple reasons desirable to work with models

through their likelihood functions and develop generic methods for analyzing and manipulating them

◮ The structure of R as a functional programming language

makes this feasible in a fairly straightforward way

14 / 14

Summary

◮ Standard modelling tools are great, but can be limiting ◮ It is for multiple reasons desirable to work with models

through their likelihood functions and develop generic methods for analyzing and manipulating them

◮ The structure of R as a functional programming language

makes this feasible in a fairly straightforward way

14 / 14