Following Sroka 2018, likelihood calculation for a nb glm model with a log odds link This provides odds ratios for count data!

nb_logodds.r

#Y is an Nx1 vector of the phenotype
#x is a NxN matrix of covariates
#b is beta (supplied by optim)
#d is a dispersion parameter (supplied by optim)
likfun1 <- function(y,x,theta)
{ 
 #the beta parameters will be the first N-1 inputs supplied by optim, dispersion the last
 b=theta[-length(theta)]
 d=theta[length(theta)]
 loglik1 <- sum ( log(gamma(y+d)) -log(gamma(y + 1)) - log(gamma(d)) +
 y*log((1 + exp(x%*%b))^(1/d) - 1) - 
 (1 + y/d)*log(1 + exp(x%*%b)) )
 return(-loglik1)
}

derivbeta <- function (y,x,theta)
{ 
 
 ymat <- cnv_fam_subset2$pheno - min(ymat) #have to set it to zero!
 xmat <- model.matrix (~ C1+C2+C3+C4+C5+LRR_SD+cnv,data=cnv_fam_subset2)
 
res1 <- optim(c(rep(0,8),1), likfun1, y=ymat,x=xmat ,hessian=TRUE,method="BFGS")
OI<-solve(res1$hessian)
se <- sqrt(diag(OI))
tv<-res1$par/se
pval<-2*pnorm(tv,lower.tail=F)
results<-cbind(res1$par,se,tv,pval)
results(colnames)<-c("b","se","t","p")
results(rownames)<-c("Const","C1")
print(results,digits=3)