I have activity data for 6 individuals (ID
) obtained using two different formulas (RMS.X16
and VeDBA.X16
), and I want to test which of those formulas fit better to my response variable RMS.V13AP
.
Here you can see the relationship of RMS.X16
and VeDBA.X16
with RMS.V13AP
.
Below you can see also a histogram of my response variable RMS.V13AP
(upper left). In the graph appears also the histogram of my response variable with different transformations.
As you can see, RMS.V13AP
is positive and continuous, so I decided to make a GLM using a Gamma distribution with a log link. I included ID
as a fixed effect and not as a random effect because I want to test if the relationship changes among individuals. To answer the main question (differences between formulas for predicting RMS.V13AP
), I did this:
mod1 <- glm(RMS.V13AP ~ VeDBA.X16+ID,data =
FormulaValidation.57s, family=Gamma(link=log))
mod2 <- glm(RMS.V13AP ~ RMS.X16+ID, data = FormulaValidation.57s,
family=Gamma(link=log))
AIC(mod1,mod2)
df AIC
mod1 8 3831.017
mod2 8 3812.568
ED.mod1 <- 100* (1-(mod1$deviance/mod1$null.deviance)) # Explained deviance
ED.mod1
[1] 62.3574
ED.mod2 <- 100* (1-(mod2$deviance/mod2$null.deviance))
ED.mod2
[1] 62.54075
It seems that the formula RMS
fits better my data, however differences between RMS
and VeDBA
are very low according to the explained deviance.
My doubts arise when I make some diagnostic plots. First, I show the plots I usually find in literature:
mod1.diag <- glm.diag(mod1)
mod2.diag <- glm.diag(mod2)
glm.diag.plots(mod1,mod1.diag)
glm.diag.plots(mod2,mod2.diag)
After some searching, I found also this approach to explore residuals (simulations):
library(DHARMa)
sim_nbz <- simulate(mod2, nsim =1000)
str(sim_nbz)
sim_nbz = do.call(cbind, sim_nbz)
head(sim_nbz)
sim_res_nbz = createDHARMa(simulatedResponse = sim_nbz,
observedResponse =
FormulaValidation.57s$RMS.V13AP,
fittedPredictedResponse =
predict(mod2),
integerResponse = FALSE)
plotSimulatedResiduals(sim_res_nbz)
The way I see it, there is clear evidence of residual patterns but I don't know what should be next step.
Should I care about those residual patterns? Should I transform my response variable and try to use another type of regression?