Hi all

Really appreciate a hand here. I am trying to model some data with a bionomial outcome. I believe that I need to use a negative bionomial distribution as there were a lot of samples where a large number of zeros were present i.e. none sick and when modelled using a bionomial distribution in the lme4 package the residuals were extremely high. Hence the attempted use of gamlss package.

I have had difficulty coding the model for the gamlss package, perhaps I have done something wrong. Also I would like to include the denominator in the outcome as the sample size varied per group i.e. in the lme4 package I coded it as: glmer(cbind(dat$NSick,dat$Ntest); but couldn't seem to do this in gamlss.

The data below is "made up data", but reflects the analysis I am doing i.e. data clustered within farm (hence the need for a random effect), Lots of zero outcome (and the need to use a negative bionomial equation).


Sorry in advance for this torturous code; although I use R I am not greatly proficient and am not a statistician.


Breed=as.factor(c(rep("Merino",10),rep("Romney",20),rep("Perindale",0),
rep("Merino",10),rep("Romney",10),rep("Perindale",10),rep("Merino",15),
rep("Romney",5),rep("Perindale",10),rep("Merino",10),
rep("Romney",10),rep("Perindale",10)))
Farm=c(rep("one",30),rep("two",30),rep("three",30),rep("four",30))
NSick=as.numeric(c(c(rep(1,5)),c(rep(0,25)),c(rep(1,3)),c(rep(0,16)),c(rep(1,3)),
c(rep(0,18)),c(rep(1,10)),c(rep(0,15)),c(rep(1,3)),c(rep(0,19)),c(rep(1,3))))
Ntest=as.numeric(c(c(rep(20,5)),c(rep(18,25)),c(rep(20,3)),c(rep(15,16)),c(rep(6,3)),
c(rep(19,18)),c(rep(20,10)),c(rep(20,15)),c(rep(18,3)),c(rep(20,19)),c(rep(17,3))))

dat=data.frame(Farm,Breed, NSick, Ntest)
dat

# Bionomial mixed model

library(lme4)

sick.glmm <- glmer(cbind(dat$NSick,dat$Ntest) ~ Breed + (1 | Farm), data = dat,
family = binomial)
summary(sick.glmm)


# Random-intercept negative binomial mixed model

library(gamlss)

rand.in.mx <- gamlss(NSick ~ Breed + (1 | Farm),
                                                            data = dat,
                                                            family = NBI
                                                            )
# Code seems to fall over and couldn't seem to add a denominator ie Ntest


Kind regards

Andrew McFadden | Ministry for Primary industries | New Zealand


This email message and any attachment(s) is intended solely for the addressee(s)
named above. The information it contains is confidential and may be legally
privileged.  Unauthorised use of the message, or the information it contains,
may be unlawful. If you have received this message by mistake please call the
sender immediately on 64 4 8940100 or notify us by return email and erase the
original message and attachments. Thank you.

The Ministry for Primary Industries accepts no responsibility for changes
made to this email or to any attachments after transmission from the office.

	[[alternative HTML version deleted]]

