R help - Apr 2012 - Specifying special poisson maximum likelihood

Hi everyone

I am stuck on specifying my own maximum likelihood function for a
special poisson model.

My poisson model is as follow: O ~ Pois(b*N + b*RR*E)

With
O = observed cases
b = constant (known)
N = number of unexposed persons (known)
E = number exposed persons (known)
RR = relative risk (value is assumed under a scenario, e.g. RR=2.0)

I used rpois to simulate the values of O for several years with known
values of b, N,  E and RR.

Hence, an example simulated data set would look like:

O <- c(173, 184, 154, 162, 158, 196, 203, 230, 240, 252, 252, 268, 255,
297, 265, 270, 265, 283)

N <-
c(1345246.66,1325945.58,1311490.16,1301155.99,1212731.00,1052598.42,889413.52,592314.47,31694.02,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00)

E <-
c(1387.337,7294.417,13634.838,24372.007,120172.455,290369.579,471271.484,796911.525,1388343.975,1445355.000,1464523.000,1482467.000,1492197.000,1494862.000,1489437.000,1478092.000,1468076.000,1458794.000)


In this dataset, I now want to estimate RR under the above model with
known values for E and N (the ones above) and the constant b should be
treated as nuisance parameter. In addition, no intercept should be
estimated. The model formula would thus be: O = b1*N + b2*E
Then, I would use a nonlinear combination of the estimates to calculate the RR:
RR = (b2/b1)

Can anybody help me specifying my own likelihood function to be used
with optim?

That would help me alot. Thank you.

Denis
_____________________________________________


Denis Aydin, MSc

Swiss Tropical and Public Health Institute
Socinstrasse 57, P.O. Box
4002 Basel, Switzerland

T +41 (0)61 284 86 09
F +41 (0)61 284 81 01
denis.aydin@unibas.ch

_____________________________________________

Postal address:

Denis Aydin, MSc
Swiss Tropical and Public Health Institute
Socinstrasse 57, P.O. Box
4002 Basel, Switzerland

Visitors address:

Eulerstrasse 68
4051 Basel, Switzerland

_____________________________________________

www.swisstph.ch

	[[alternative HTML version deleted]]

<Denis.Aydin <at> unibas.ch> writes:
> 
> Hi everyone
> 
> I am stuck on specifying my own maximum likelihood function for a
> special poisson model.
> 
> My poisson model is as follow: O ~ Pois(b*N + b*RR*E)
> 
> With
> O = observed cases
> b = constant (known)
> N = number of unexposed persons (known)
> E = number exposed persons (known)
> RR = relative risk (value is assumed under a scenario, e.g. RR=2.0)
> 
> I used rpois to simulate the values of O for several years with known
> values of b, N,  E and RR.
> 
> In this dataset, I now want to estimate RR under the above model with
> known values for E and N (the ones above) and the constant b should be
> treated as nuisance parameter. In addition, no intercept should be
> estimated. The model formula would thus be: O = b1*N + b2*E
> Then, I would use a nonlinear combination of the estimates to calculate the
RR: RR = (b2/b1)

  How about:

dat <- data.frame(O=c(173, 184, 154, 162,
                  158, 196, 203, 230, 240, 252, 252, 268, 255,
                  297, 265, 270, 265, 283),
                  N=c(1345246.66,1325945.58,1311490.16,
                  1301155.99,1212731.00,1052598.42,889413.52,
                  592314.47,31694.02,0.00,
                  0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00),
                  E=c(1387.337,7294.417,13634.838,24372.007,
                  120172.455,290369.579,471271.484,796911.525,
                  1388343.975,1445355.000,1464523.000,
                  1482467.000,1492197.000,
                  1494862.000,1489437.000,1478092.000,1468076.000,1458794.000))

## suppose  b=1

library(bbmle)
b <- 1
mle2(O~dpois(lambda=b*(N+exp(logRR)*E)),
     start=list(logRR=0),data=dat)


?