packages
##
library(ggbreak)
library(scales)
library(tidyverse)
library(googlesheets4)
library(googledrive)
library(plotly)
library(lme4)
library(fitdistrplus)
library(goft)
library(data.table)
Data
ss= "https://docs.google.com/spreadsheets/d/1dr6LCQJevHdeS08zrdkbeONfhxU7etsojqVvBU7ryfs/edit?usp=sharing"
hoja = 1
rango = "A1:N193"
Descriptive
Proportion of infested pigs
## infested n
## 1 Si 47
## 2 No 70
Number of “quistes”
Generalized Logistic Model
mod <- glm(data=df,
fertile.bin ~ size_cm * Sexo , family = "binomial")
summary(mod) # No interaction
##
## Call:
## glm(formula = fertile.bin ~ size_cm * Sexo, family = "binomial",
## data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.941051 0.789402 -3.726 0.000195 ***
## size_cm 0.431272 0.390205 1.105 0.269054
## SexoMacho 1.275366 0.892428 1.429 0.152976
## size_cm:SexoMacho -0.001491 0.429234 -0.003 0.997228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 213.71 on 191 degrees of freedom
## Residual deviance: 194.76 on 188 degrees of freedom
## AIC: 202.76
##
## Number of Fisher Scoring iterations: 5
#
mod <- glm(data=df,
fertile.bin ~ Sexo , family = "binomial")
summary(mod) # Sex association
##
## Call:
## glm(formula = fertile.bin ~ Sexo, family = "binomial", data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3224 0.4685 -4.957 7.16e-07 ***
## SexoMacho 1.5168 0.5039 3.010 0.00261 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 213.71 on 191 degrees of freedom
## Residual deviance: 201.84 on 190 degrees of freedom
## AIC: 205.84
##
## Number of Fisher Scoring iterations: 4
##
mod <- glm(data=df,
fertile.bin ~ size_cm + Sexo , family = "binomial")
summary(mod)
##
## Call:
## glm(formula = fertile.bin ~ size_cm + Sexo, family = "binomial",
## data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9391 0.5424 -5.418 6.02e-08 ***
## size_cm 0.4300 0.1626 2.645 0.00817 **
## SexoMacho 1.2728 0.5165 2.464 0.01373 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 213.71 on 191 degrees of freedom
## Residual deviance: 194.76 on 189 degrees of freedom
## AIC: 200.76
##
## Number of Fisher Scoring iterations: 5
Generalized Logistic Mixed Model
## Justifying Random effects inclution (Baseline glm Vs. baseine mixed-model )
m0.glm <- glm(data= df,
fertile.bin ~ 1, family = "binomial")
m0.glmer <- glmer(data=df,
fertile.bin ~ (1|code), family = "binomial")
##
aic.glmer <- AIC(logLik(m0.glmer))
aic.glm <- AIC(logLik(m0.glm))
aic.glmer ; aic.glm
## [1] 160.6177
## [1] 215.7116
#
null.id = -2 * logLik(m0.glm) + 2 * logLik(m0.glmer)
pchisq(as.numeric(null.id), df=1, lower.tail = FALSE)
## [1] 4.155105e-14
- Se justifica el uso del efecto aleatorio (cerdo), ya que el
AIC del modelo basal (nulo) con el efecto elatorio es menor (P<
0.001) que sin dicho efecto aleatorio. *
### SEX
mod <- glmer(data=df,
fertile.bin ~ (1|code) + Sexo, family = "binomial") # "Pig" as `random effect`
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fertile.bin ~ (1 | code) + Sexo
## Data: df
##
## AIC BIC logLik deviance df.resid
## 162.1 171.9 -78.0 156.1 189
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1186 -0.2690 -0.1283 -0.1030 3.7174
##
## Random effects:
## Groups Name Variance Std.Dev.
## code (Intercept) 5.794 2.407
## Number of obs: 192, groups: code, 39
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.1193 1.3264 -3.106 0.0019 **
## SexoMacho 0.8578 1.1719 0.732 0.4642
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SexoMacho -0.561
## Size-CM
mod <- glmer(data=df,
fertile.bin ~ (1|code) + size_cm , family = "binomial") # "Pig" as `random effect`
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fertile.bin ~ (1 | code) + size_cm
## Data: df
##
## AIC BIC logLik deviance df.resid
## 158.3 168.1 -76.2 152.3 189
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3304 -0.2820 -0.1102 -0.0700 3.1833
##
## Random effects:
## Groups Name Variance Std.Dev.
## code (Intercept) 6.386 2.527
## Number of obs: 192, groups: code, 39
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.8529 1.3814 -3.513 0.000443 ***
## size_cm 0.6644 0.3358 1.979 0.047857 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## size_cm -0.592
## `geom_smooth()` using formula = 'y ~ x'
Generalized Linear Model: GROUP comparison
mod <- glmer(data=df,
size_cm ~ as.factor(fertile.bin) + (1|code) , family = "Gamma")
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( inverse )
## Formula: size_cm ~ as.factor(fertile.bin) + (1 | code)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 385.6 398.6 -188.8 377.6 188
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25221 -0.45171 -0.01108 0.42757 2.48162
##
## Random effects:
## Groups Name Variance Std.Dev.
## code (Intercept) 0.1289 0.3590
## Residual 0.1656 0.4069
## Number of obs: 192, groups: code, 39
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 1.14513 0.13502 8.481 <2e-16 ***
## as.factor(fertile.bin)1 -0.04955 0.03756 -1.319 0.187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## as.fctr(.)1 -0.010
#
mod <- glmer(data=df,
size_cm ~ Sexo + (1|code), family = "Gamma")
summary(mod) # Sex association
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( inverse )
## Formula: size_cm ~ Sexo + (1 | code)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 387.3 400.3 -189.6 379.3 188
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.24258 -0.47387 -0.03492 0.39087 2.75095
##
## Random effects:
## Groups Name Variance Std.Dev.
## code (Intercept) 0.1303 0.3609
## Residual 0.1671 0.4088
## Number of obs: 192, groups: code, 39
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 1.16202 0.18130 6.409 1.46e-10 ***
## SexoMacho -0.04097 0.22134 -0.185 0.853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SexoMacho -0.666
