# Fit a M_tbh model to the Perom data:
(pdata <- aux.posbernoulli.t(with(Perom, cbind(y1, y2, y3, y4, y5, y6))))
Perom <- data.frame(Perom,
bei = 0, # Add this
pdata$cap.hist1) # Incorporate these
head(Perom) # Augmented with behavioural effect indicator variables
tail(Perom)
Run the code above in your browser using DataLab