[R] betabinomial model

[Zaihra T]

   Hi,

   can anyone help me fit betabinomial model to the following dataset where
   each iD is a cluster in itself , if i use package aod 's betabinom model it
   gives an estimate of zero to phi(the correlation coeficient ) and if i fix
   it to the anova type estimate obtained from icc( in package aod) then it
   says system is exactly singular. And when i try to fit my loglikelihood by
   writing a function directly for loglikelihood and try to optimise using
   optim  it  gives  warnings  as well as error  and if i use NLMINB then
   convergence is false.

   Can  someone  save my soul..... pl!!!!!!!!z below is my program script
   ..............

   #ovarian cancer data and parity among 769 probands and their mothers

   #cancer in proband

   y1<-c(1,1,1,1,1,1,1,rep(1,29),rep(1,36),rep(1,41),r!
   ep(1,85),rep(1,105),rep(1,84),

   1,rep(0,22),rep(0,33),rep(0,38),rep(0,50),rep(0,103),rep(0,135))

   #cancer in mother

   y2<-c(1,1,1,1,1,1,1,rep(0,29),rep(0,36),rep(0,41),rep(0,85),rep(0,105),rep(0
   ,84),

   1,rep(0,22),rep(0,33),rep(0,38),rep(0,50),rep(0,103),rep(0,135))

   # total events in each cluster

   y<-y1+y2

   n<-rep(2,769)

   z<-rep(1,length(y))

   n1<-length(y)

   # number of childbirths in mother (0 for 1-2,1 for 3+)

   z1<-c(0,0,1,1,1,1,1,rep(0,29),rep(0,36),rep(0,41),rep(1,85),rep(1,105),rep(1
   ,84),

   1,rep(0,22),rep(0,33),rep(0,38),rep(1,50),rep(1,103),rep(1,135))

   # of child births in mother three catag 2 dummy variables z2=1(1-2) else
   0,z3=1(o birth)else0so if z2

   # z3 both r zero its 3 births catageroy

   z2<-c(0,1,0,0,0,1,0,rep(0,29),rep(1,36),rep(0,41),rep(0,85),rep(1,105),rep(0
   ,84),

   0,rep(0,22),rep(1,33),rep(0,38),rep(0,50),rep(1,103),rep(0,135))

   z3<-c(1,0,1,1,1,0,0,r!
   ep(1,29),rep(0,36),rep(0,41),rep(1,85),rep(0,105),rep(0,84),

   1,r ep(1,22),rep(0,33),rep(0,38),rep(1,50),rep(0,103),rep(0,135))


   cancer<-as.data.frame(cbind(y1,y2,y,z,z1,z2,n))

   cancer

   # fitting betabinomial model accounting for icc-----------------

   beta.binomial<- betabin(cbind(y, n - y) ~ z1 +z2+z3, ~ 1, data =cancer)

   beta.binomial

   # calculating ICC using REML/ML method-------------

   intra.clus<-icc(n,y,data=cancer)

   intra.clus

   a1<-0

   a2<-0

   #loglikelihood function for ordinary betabinomial model------------

   b.bin<-function(th){

   beta0<-th[1]

   beta1<-th[2]

   beta2<-th[3]

   beta3<-th[4]

   rho<-th[5]

   x<-cbind(z,z1,z2,z3)

   beta<-c(beta0,beta1,beta2,beta3)

   p<-exp(x%*%beta)

   q<-(1+x%*%beta)

   for(i in 1:n1)

   {

   for (j in 0:y[i]-1)

   {

   a1<-a1+log(ifelse(y[i]-1<0,1,p[i]+rho*q*j))

   #a1<-a1+log(ifelse(y[i]-1<0 ||p[i]+rho*q*j<=0,1,p[i]+rho*q*j! ))

   }

   }

   for(i in 1:n1)

   {

   for (j in 0:1-y[i])

   {

   a2<-a2+log(ifelse(1-y[i]<0,1,1+rho*q*j))

   #a2<-a2+log(ifelse(1-y[i]<0 ||1+rho*q*j<=0 ,1,1+rho*q*j))

   }

   }

   out<--(a1+a2-n1*log(ifelse(1+rho<=0,1,1+rho)))

   }

   d<--b.bin(c(-1.2628,-0.0936,0.3485,0.6155,-.2918))

   +sum(log(factorial(2)/factorial(y)*factorial(2-y)))

   d

   betabinom<- optim(c(-1.2628,-0.0936,0.3485,0.6155,-.2918),b.bin, method =
   "BFGS")

   betabinom<-nlminb(start=c(-1.2628,-0.0936,0.3485,0.6155,-.2918),b.bin)

   betabinom


   beta.binomial<-  betabin(cbind(y,  n  -  y)  ~  z1  +z2+z3,  ~ 1, data
   =cancer,fixpar=list(5,-.2918))

   beta.binomial