#trying to solve question 'can you log-transform and specify Gamma in the same model' question
ToadsBd<-read.table(file.choose(),header=T)
list(ToadsBd)
#first see how well treatment group predicts Bd score with non-log transformed data
mod1<-glm(Bd~factor(group))
summary(mod1)
#massively overdispersed. Are the data non-normal?
shapiro.test(Bd)
W = 0.3652, p-value = 5.666e-13
#yes, definitely non-normal
#try log-transforming data and see if that helps 
plot(qqnorm(Bd),log="y")
#log plot straightens it out, almost, so yes log-transform helps
#try model again with log transformed Bd score
mod2<-glm(logBd~factor(group))
summary(mod2)
#a big improvement but still overdispersed
#other options - specify an error family? Looks like original data are Gamma distributed
#should test if variance increases or remains constant with mean on scale of the original, non-logged data
par(mfrow=c(2,2))
plot(mod1)
#can you tell this from a diagnostic plot? Not sure how. If not, how do you assess this?
#in the meantime, assume it does and try Gamma (using default link = reciprocal) with non-logged data
mod3<-glm(Bd~factor(group),family=Gamma)
summary(mod3)
#mod3 is a major improvement on mod1 and less dispersed than mod2 but has a much larger AIC than mod2
#is it valid to specify Gamma in a model where the data have been log-transformed?
#or does it have to be a choice between transformation or Gamma?
#if specify both, model is quite good, but it may not be valid. Please help!
mod4<-glm(logBd~factor(group),family=Gamma)
summary(mod4)
#residual deviance now well below df, not overdispersed and the effect of group on Bd is significant
#I would also like to include assessment of the effect of site, but this is a random effect requiring a mixed model
#I cannot find a mixed model that works with Gamma errors. What can I do?
