The package 'VGAM' (all caps) has a betabinomial model fitting facility in it.
(I doubt this will save your soul, by the way, but it may save your butt. I think that is probably your most immediate concern, though.) -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Zaihra T Sent: Thursday, 20 March 2008 2:52 AM To: [EMAIL PROTECTED] Subject: [R] betabinomial model 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),r ep(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),r ep(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),r ep(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 ______________________________________________ R-help@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. ______________________________________________ R-help@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.
