I have a dataset with 5 replicates (batch) of a 6 treatments (trt), each treatment has 5 associated treatment times (time) with the treatment time for controls being 0. Several treatments in batch 1 failed and were dropped from the dataset. From each treatment x time combination there is an associated number of "bad" results (bad) out of a total sample (total). The total sample is 10 for almost all of the treatments, two treatments had smaller samples due to sampling errors.
I am trying to fit a negative binomial model with bad as the response variable and trt and time as predictors. The data look like this:
library(dplyr)
# create dataset
data.frame(
batch=c(rep("A",11),rep("B",26),rep("C",26),rep("D",26),rep("E",26)),
trt=c("ctrl",rep("trt1",5),rep("trt5",5),rep(c("ctrl",rep("trt1",5),rep("trt2",5),rep("trt3",5),rep("trt4",5),rep("trt5",5)),4)),
time=c(0,c(2:6)*5,c(1:5)*5,rep(c(0,rep(c(2:6)*5,3),rep(c(1:5)*5,2)),4)),
bad=c(1,0,1,0,0,0,0,2,0,0,0,3,3,0,2,0,1,2,3,2,1,1,1,1,3,3,0,0,1,2,1,5,1,2,5,2,1,0,0,2,1,3,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,2,0,1,0,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,1,3,1,3,0,2,1,1,3,1,1,1,4,1,3,3,0,0,0,3,1,1,1,3),
total=c(rep(10,21),9,rep(10,11),8,rep(10,81))
)%>%
# factors for model
mutate(trt=factor(trt,
levels=c("ctrl","trt1","trt2","trt3","trt4","trt5")),
time=factor(time,
levels=c(0:6)*5))->df
I used the MASS package to fit a model as shown here:
# fit model
MASS::glm.nb(bad~time*trt+offset(log(total)),
data=df)->mod
summary(mod)
To check the fit of the model I used DHARMa:
# check model fit
DHARMa::simulateResiduals(mod)%>%plot()
and got the resulting plot, which seems to imply a good fit:
I then used emmeans to compare the EMMs to the means from the raw dataset, but found that the standard errors from the model were very high. There are results that seem significant from the graph, but the p-values are relatively high (I'm guessing from the high standard error?)
se<-function(x,na=T)# x is a numeric list, na functions the same as na.rm (defaults to T)
{
if(na==T)
{x<-na.omit(x)}
return(sd(x)/sqrt(length(x)))
}
# check model compared to data
# get means and standard errors for dataset
df%>%
# standardize by total sample since some were less than 10
mutate(bad.std=bad/total*10)%>%
# means and SEs
group_by(trt,time)%>%
summarise(means=mean(bad.std),
ses=se(bad.std),
.groups="keep")->y
# get emmeans and standard errors from model
emmeans::emmeans(mod,
trt.vs.ctrl~trt*time)$emmeans%>%
as.data.frame()%>%
# drop NA cols
subset(!is.na(emmean))%>%
# back transform from log()
mutate(mean.trans=exp(emmean),
se.trans=exp(SE))%>%
# combine with raw dataset
inner_join(y,by=c("trt","time"))%>%
# reorganize for plotting
mutate(trans=paste(mean.trans,se.trans),
og=paste(means,ses))%>%
# declutter columns
select(trt,time,trans,og)%>%
pivot_longer(cols=c(trans,og),
names_to="dataset")%>%
rowwise()%>%
mutate(means=str_split(value," ")[[1]][1]%>%as.numeric(),
ses=str_split(value," ")[[1]][2]%>%as.numeric())%>%
# plot to compare
ggplot(aes(x=time,
fill=dataset))+
geom_col(aes(y=means),
position=position_dodge())+
geom_errorbar(aes(ymin=means-ses,
ymax=means+ses),
position=position_dodge())+
facet_grid(.~trt,
scales="free")
I would like to compare the treatments with the controls using Dunnett's test:
emmeans::emmeans(mod,
trt.vs.ctrl~trt*time)$contrasts
and compare treatments at different treatment times:
emmeans::emmeans(mod,
pairwise~trt|time)$contrasts
Is there something wrong with the model or am I using the emmeans package improperly? Looking at the bar plot of the data it would seem that some treatment times in trt3 and trt4 are significant compared to the control, but the standard error from the model is too high to show any significance. What am I missing?