# [R] Fitting a Bivariate Poisson log-normal distribution

From: Famoye, Felix <famoy1kf_at_cmich.edu>
Date: Sun, 27 Jul 2008 15:36:41 -0400

I wrote an R program to fit a bivariate Poisson log-normal model. The model is first proposed by Aitchison and Ho (1989) in Biometrika, Volume 76 pages 643-653. The model involves using a double integration. The following is the program and the data that I used. My major problem was that R was running for a long time (more than 3 hours) with no results. I would like to place some printing at some places to see if R was doing something. Unfortunately, there was nothing printed when I included printing in the function "f". Note that I have used the package "adapt" for double integration. Can someone provide some help on this problem? I would like to know if R is doing anything at all in fitting the model. Thank you for your time.

#This program will be used to estimate the parameters of
Poisson-lognormal model.

y1 = yd[[1]]
y2 = yd[[2]]
n = length(y1)
x0= rep(1, n)

```y1  = cbind(y1)
y2  = cbind(y2)
yy=cbind(y1,y2)
xx = cbind(x0,yd[[3]])
av1=mean(y1)
```

av2=mean(y2)
cov=var(yy)
mu=cbind(av1, av2)
```rho=cov[1,2]/sqrt(cov[1,1]*cov[2,2])
sd1=sqrt(cov[1,1])
sd2=sqrt(cov[2,2])
oth=rbind(sd1,sd2,rho)
```

z = vector(length=2) #This is a column vector with two rows

#To find the inital estimates for the regression coefficients

```para1 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(y1+0.5))
para2 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(y2+0.5))
parab = rbind(para1, para2, oth)
```

parab

#---------------------------------------------------------------------
# Fitting Poisson lognorma model
#f is the function that defines the (negative) log likelihood
f <- function(parab,y1,y2,xx)
{
vv1=(parab[5])^2
vv2=(parab[6])^2
para.1= cbind(parab[1:2])
para.2= cbind(parab[3:4])
mu.1 = xx%*%para.1 - vv1/2
#mu.1 = exp(mu.1)
mu.2 = xx%*%para.2 - vv2/2
#mu.2 = exp(mu.2)

mu.v=cbind(mu.1,mu.2)

```cov[1,1]=vv1
cov[2,2]=vv2
cov[1,2]=parab[7]*parab[5]*parab[6]
cov[2,1]=cov[1,2]
```

n=length(y1)
va=rep(0.0, n)
for (i in 1:n)
{
fred=function(z)
{
```tz=log(z)-mu.v[i,]
qq=(t(tz)%*%ginv(cov))%*%tz
fq=(exp(-qq/2))/(2*pi*z[1]*z[2]*sqrt(det(cov)))
p1=((z[1]^y1[i])*exp(-z[1]))/gamma(y1[i]+1)
p2=((z[2]^y2[i])*exp(-z[2]))/gamma(y2[i]+1)
```
int=fq*p1*p2
}
}
-sum(va)
}

I.p = nlminb(start=parab,objective=f,y1=y1,y2=y2,xx=xx) print("Poisson Log-normal Regression Model") I.p

0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
1 0 0
1 0 0
1 0 0
1 0 0
1 0 0
1 0 0
1 0 0
2 0 0
2 0 0
2 0 0
2 0 0
2 0 0
2 0 0
3 0 0
3 0 0
3 0 0
3 0 0
3 0 0
3 0 0
4 0 0
4 0 0
4 0 0
5 0 0
5 0 0
7 0 0
8 0 0
9 0 0
9 0 0
12 0 0
13 0 0
14 0 0
16 0 0
20 0 0
0 1 0
0 1 0
0 1 0
0 1 0
0 1 0
0 1 0
0 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
1 1 0
2 1 0
2 1 0
2 1 0
2 1 0
2 1 0
2 1 0
2 1 0
3 1 0
3 1 0
3 1 0
4 1 0
4 1 0
5 1 0
5 1 0
5 1 0
9 1 0
10 1 0
10 1 0
15 1 0
0 2 0
0 2 0
0 2 0
0 2 0
0 2 0
0 2 0
0 2 0
0 2 0
0 2 0
1 2 0
1 2 0
1 2 0
1 2 0
1 2 0
1 2 0
1 2 0
1 2 0
1 2 0
2 2 0
2 2 0
2 2 0
3 2 0
4 2 0
4 2 0
5 2 0
6 2 0
6 2 0
6 2 0
7 2 0
7 2 0
12 2 0
15 2 0
26 2 0
0 3 0
0 3 0
0 3 0
0 3 0
0 3 0
0 3 0
0 3 0
0 3 0
0 3 0
1 3 0
1 3 0
1 3 0
1 3 0
1 3 0
2 3 0
2 3 0
4 3 0
5 3 0
7 3 0
10 3 0
0 4 0
0 4 0
0 4 0
0 4 0
1 4 0
1 4 0
1 4 0
1 4 0
3 4 0
3 4 0
0 5 0
0 6 0
0 6 0
1 6 0
3 6 0
0 7 0
0 7 0
1 7 0
0 8 0
1 8 0
0 9 0
0 9 0
0 9 0
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
1 0 1
2 0 1
2 0 1
2 0 1
2 0 1
2 0 1
2 0 1
3 0 1
3 0 1
3 0 1
3 0 1
3 0 1
3 0 1
3 0 1
3 0 1
4 0 1
4 0 1
4 0 1
4 0 1
4 0 1
4 0 1
4 0 1
4 0 1
5 0 1
5 0 1
7 0 1
7 0 1
8 0 1
9 0 1
9 0 1
11 0 1
15 0 1
15 0 1
19 0 1
20 0 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
0 1 1
1 1 1
1 1 1
1 1 1
2 1 1
2 1 1
2 1 1
2 1 1
2 1 1
4 1 1
4 1 1
5 1 1
5 1 1
5 1 1
6 1 1
8 1 1
8 1 1
8 1 1
8 1 1
9 1 1
23 1 1
23 1 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
0 2 1
1 2 1
1 2 1
1 2 1
1 2 1
1 2 1
1 2 1
1 2 1
1 2 1
1 2 1
1 2 1
2 2 1
2 2 1
2 2 1
3 2 1
3 2 1
3 2 1
4 2 1
4 2 1
6 2 1
0 3 1
0 3 1
0 3 1
0 3 1
0 3 1
0 3 1
0 3 1
1 3 1
1 3 1
1 3 1
1 3 1
1 3 1
2 3 1
2 3 1
2 3 1
3 3 1
6 3 1
0 4 1
0 4 1
0 4 1
0 4 1
0 4 1
1 4 1
1 4 1
1 4 1
1 4 1
2 4 1
5 4 1
0 5 1
0 5 1
0 5 1
1 5 1
1 5 1
0 6 1
0 6 1
1 7 1
1 7 1
1 8 1
0 10 1
0 10 1

Felix Famoye
Department of Mathematics
Central Michigan University
Mt. Pleasant, Michigan 48859-0001

e-mail: felix.famoye_at_cmich.edu
web site: http://www.cst.cmich.edu/users/famoy1kf/ voice: (989)774-5497, fax: (989)774-2414

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 27 Jul 2008 - 20:06:49 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 27 Jul 2008 - 21:33:08 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.