[R] Specifying special poisson maximum likelihood

[Ben Bolker]

 <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)


?