model_global<-lm(dt$liver_residuals ~ dt$treatments + dt$Tenebrio_sum_week4)
summary(model_global)
##
## Call:
## lm(formula = dt$liver_residuals ~ dt$treatments + dt$Tenebrio_sum_week4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.032280 -0.009541 -0.000973 0.007450 0.035571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.015383 0.014139 -1.088 0.2850
## dt$treatmentsDose 1 0.008926 0.008981 0.994 0.3280
## dt$treatmentsDose 2 0.015331 0.008957 1.712 0.0970 .
## dt$treatmentsDose 3 0.019092 0.008666 2.203 0.0351 *
## dt$treatmentsDose 4 0.010413 0.008632 1.206 0.2368
## dt$treatmentsDose 5 0.017585 0.008958 1.963 0.0587 .
## dt$Tenebrio_sum_week4 0.001426 0.005423 0.263 0.7943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01551 on 31 degrees of freedom
## Multiple R-squared: 0.1756, Adjusted R-squared: 0.01599
## F-statistic: 1.1 on 6 and 31 DF, p-value: 0.3845
# Obtain estimated marginal means
emmeans_object <- emmeans(model_global, ~ treatments)
# Perform pairwise post-hoc comparisons with Tukey adjustment
pairwise_comparisons <- contrast(emmeans_object, method = "pairwise", adjust = "tukey")
# Summarize and view the results
summary(pairwise_comparisons)
## contrast estimate SE df t.ratio p.value
## Control - Dose 1 -0.00893 0.00898 31 -0.994 0.9164
## Control - Dose 2 -0.01533 0.00896 31 -1.712 0.5348
## Control - Dose 3 -0.01909 0.00867 31 -2.203 0.2651
## Control - Dose 4 -0.01041 0.00863 31 -1.206 0.8304
## Control - Dose 5 -0.01758 0.00896 31 -1.963 0.3854
## Dose 1 - Dose 2 -0.00640 0.00898 31 -0.713 0.9789
## Dose 1 - Dose 3 -0.01017 0.00875 31 -1.162 0.8511
## Dose 1 - Dose 4 -0.00149 0.00866 31 -0.172 1.0000
## Dose 1 - Dose 5 -0.00866 0.00899 31 -0.963 0.9261
## Dose 2 - Dose 3 -0.00376 0.00867 31 -0.434 0.9979
## Dose 2 - Dose 4 0.00492 0.00863 31 0.570 0.9923
## Dose 2 - Dose 5 -0.00225 0.00896 31 -0.252 0.9998
## Dose 3 - Dose 4 0.00868 0.00832 31 1.043 0.8995
## Dose 3 - Dose 5 0.00151 0.00866 31 0.174 1.0000
## Dose 4 - Dose 5 -0.00717 0.00863 31 -0.831 0.9594
##
## P value adjustment: tukey method for comparing a family of 6 estimates
trt_vs_ctrl_comparisons <- contrast(emmeans_object, method = "trt.vs.ctrl")
summary(trt_vs_ctrl_comparisons)
## contrast estimate SE df t.ratio p.value
## Dose 1 - Control 0.00893 0.00898 31 0.994 0.7422
## Dose 2 - Control 0.01533 0.00896 31 1.712 0.3205
## Dose 3 - Control 0.01909 0.00867 31 2.203 0.1343
## Dose 4 - Control 0.01041 0.00863 31 1.206 0.6143
## Dose 5 - Control 0.01758 0.00896 31 1.963 0.2108
##
## P value adjustment: dunnettx method for 5 tests
summary(model_global)
##
## Call:
## lm(formula = dt$liver_residuals ~ dt$treatments + dt$Tenebrio_sum_week4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.032280 -0.009541 -0.000973 0.007450 0.035571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.015383 0.014139 -1.088 0.2850
## dt$treatmentsDose 1 0.008926 0.008981 0.994 0.3280
## dt$treatmentsDose 2 0.015331 0.008957 1.712 0.0970 .
## dt$treatmentsDose 3 0.019092 0.008666 2.203 0.0351 *
## dt$treatmentsDose 4 0.010413 0.008632 1.206 0.2368
## dt$treatmentsDose 5 0.017585 0.008958 1.963 0.0587 .
## dt$Tenebrio_sum_week4 0.001426 0.005423 0.263 0.7943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01551 on 31 degrees of freedom
## Multiple R-squared: 0.1756, Adjusted R-squared: 0.01599
## F-statistic: 1.1 on 6 and 31 DF, p-value: 0.3845
#1. Check the homogeneity of the variance
#Plot predicted values vs residual values
par(mar=c(4,4,.5,.5))
plot(resid(model_global) ~ fitted(model_global),
xlab = 'Predicted values',
ylab = 'Normalized residuals')
abline(h = 0, lty = 2)
# Homogeneous dispersion of the residuals means that the assumption is respected.
#2. Check the independence of the model residuals with each covariate
# In order to check the independence of the model residuals we need to plot residuals vs each covariate of the model
par(mfrow = c(1,3), mar=c(4,4,.5,.5))
plot(resid(model_global) ~ dt$Tenebrio_sum_week4,
xlab = "Length", ylab = "Normalized residuals")
abline(h = 0, lty = 2)
# Homogeneous dispersion of the residuals around 0 means no pattern of residuals depending on the variable, therefore
# the assumption is respected!
# Note: The clusters are due to the data structure, where fish of only 5 size classes (large, small, and three groups in
# between) were captured.
#3. Check the normality of the model residuals
# Check the normality of the model residuals as residuals following a normal distribution indicate that the model is
# not biased.
hist(resid(model_global))
# The residuals are normal! This means our model is not biased.
##### Model diagnosis #####
#Check the homogeneity of the variance
#In order to check the homogeneity of the variance, we can plot predicted values vs residual values.
simulateResiduals(model_global, plot = TRUE) # QQ and Quantile test plots look good
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.964 0.996 0.512 0.696 0.876 0.156 0.48 0.556 0.384 0.26 0.244 0.308 0.744 0.752 0.2 0.704 0.46 0.54 0.984 0.424 ...
#Note: Homogeneous dispersion of the residuals means that the assumption is respected.
#Check for dispersion of the data
testDispersion(model_global,plot=TRUE) #GOOD
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.85943, p-value = 0.584
## alternative hypothesis: two.sided
#Note Over/underdispersion means that the observed data is more / less dispersed than expected under the fitted model
#Checking zero inflation
testZeroInflation(model_global) # Because data was transformed this test is NA
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
