MULTI-LEVEL ANALYSIS OF STRESS RESPONSES IN AMPHIBIAN OF EXTREME HABITATS

Resultados Salida de Campo 2019 - Afloramientos Rocosos Reserva Natural Bojonawi

Pregunta de investigación:

El nivel de desecación de las charcas afecta la morfología y el desempeño locomotor de renacuajos de Leptodactylus lithonaetes

Resultados:

Estos resultados están divididos en dos secciones:

  1. Los resultados morfológicos.

  2. Los resultados en el desempeño locomotor.

Resultados morfológicos

Se utilizaron 7 posturas diferentes. Cada postura se dividió en 3 tratamientos (Agua Baja AB,Agua Estable AE ,Agua Desecación AD), para un total de 21 cajas. En cada caja se introdujeron 10 renacuajos. Al final del experimento, de cada caja se tomaron entre 3 y 5 renacuajos y a cada renacuajo se le tomaron tres veces las siguientes medidas morfológicas:

  • Altura de la cabeza

  • Altura de la cola

  • Longitud de la cola

  • Perímetro del músculo

  • Longitud de la cabeza

*Eliminé la longitud total por estar relacionada directamente con la longitud de la cola y la longitud de la cabeza

knitr::include_graphics("dis_morf.png")

Analisis para evaluar la repetibilidad de las medidas tomadas en cada renacuajo:

La Base de datos

datos_morfo <- read.table("error_prueba_rep.txt",h=T)
head(datos_morfo)
##   Id_Random_Foto Tratamiento Replica Caja Individuo Medicion Altura_Cabeza
## 1             23          AB      R1    1         1        1         2.167
## 2             23          AB      R1    1         1        2         2.129
## 3             23          AB      R1    1         1        3         2.145
## 4            112          AB      R1    1         2        1         2.437
## 5            112          AB      R1    1         2        2         2.536
## 6            112          AB      R1    1         2        3         2.587
##   Longitud_Total Altura_Cola Longitud_Cola Perim_Musculo1 Longitud_Cabeza
## 1         24.217       1.653        17.420         36.499           6.577
## 2         24.611       1.741        17.326         37.880           6.661
## 3         23.892       1.691        17.248         38.119           6.600
## 4         20.958       1.268        14.001         30.291           6.998
## 5         21.001       1.336        14.207         30.549           7.177
## 6         21.011       1.341        14.029         31.187           7.144
  1. Para evaluar la repetibilidad entre las tres medidas tomadas para cada animal se realizaron modelos mixtos NULOS sin varible independiente. A partir de ese modelo se calculó el valor ICC (intraclass correlation coefficient) que permite determinar la "estabilidad en la medida" (Tomado de Finch et al., 2014: Multilevel modeling using R-p44).

ICC>0.90: Excelente

ICC entre 0,75 y 0,89: Bueno

ICC entre 0,50 y 0,74: Moderado

ICC<0.50: Malo

Se usó el paquete performance para ese cálculo de ICC.

Librerias

library(Matrix)
library(lme4)
library(lmerTest)
library(performance)

Long_cabeza:

mixtolong_ca<-lmer(Longitud_Cabeza~1+(1|Id_Random_Foto),data=datos_morfo)

summary(mixtolong_ca)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Longitud_Cabeza ~ 1 + (1 | Id_Random_Foto)
##    Data: datos_morfo
## 
## REML criterion at convergence: 396.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8609 -0.2246 -0.0074  0.1973  4.0627 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 0.9752   0.9875  
##  Residual                   0.0478   0.2186  
## Number of obs: 336, groups:  Id_Random_Foto, 112
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.24794    0.09407 111.00000   66.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtolong_ca, by_group = TRUE)
## # ICC by Group
## 
## Group          |   ICC
## ----------------------
## Id_Random_Foto | 0.953

Altura_Cabeza:

mixtoalt_cabeza<-lmer(Altura_Cabeza~1+(1|Id_Random_Foto),data=datos_morfo)
summary(mixtoalt_cabeza)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Altura_Cabeza ~ 1 + (1 | Id_Random_Foto)
##    Data: datos_morfo
## 
## REML criterion at convergence: -203.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9044 -0.3624 -0.0375  0.3089  3.5067 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 0.08257  0.2873  
##  Residual                   0.01102  0.1050  
## Number of obs: 336, groups:  Id_Random_Foto, 112
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.07796    0.02775 110.99999   74.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtoalt_cabeza, by_group = TRUE)
## # ICC by Group
## 
## Group          |   ICC
## ----------------------
## Id_Random_Foto | 0.882

Longitud_Total

mixtolong_total<-lmer(Longitud_Total~1+(1|Id_Random_Foto),data=datos_morfo)
summary(mixtolong_total)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Longitud_Total ~ 1 + (1 | Id_Random_Foto)
##    Data: datos_morfo
## 
## REML criterion at convergence: 1138.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8329 -0.1778 -0.0015  0.1830  4.2879 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 9.1205   3.0200  
##  Residual                   0.4337   0.6585  
## Number of obs: 336, groups:  Id_Random_Foto, 112
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  21.0012     0.2876 111.0000   73.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtolong_total, by_group = TRUE)
## # ICC by Group
## 
## Group          |   ICC
## ----------------------
## Id_Random_Foto | 0.955

Altura_Cola

mixtoalt_cola<-lmer(Altura_Cola~1+(1|Id_Random_Foto),data=datos_morfo)
summary(mixtoalt_cola)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Altura_Cola ~ 1 + (1 | Id_Random_Foto)
##    Data: datos_morfo
## 
## REML criterion at convergence: -483.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1861 -0.3434  0.0208  0.3360  3.3589 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 0.06755  0.25991 
##  Residual                   0.00353  0.05941 
## Number of obs: 336, groups:  Id_Random_Foto, 112
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.33883    0.02477 111.00000   54.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtoalt_cola, by_group = TRUE)
## # ICC by Group
## 
## Group          |   ICC
## ----------------------
## Id_Random_Foto | 0.950

Longitud_Cola

mixtolong_cola<-lmer(Longitud_Cola~1+(1|Id_Random_Foto),data=datos_morfo)
summary(mixtolong_cola)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Longitud_Cola ~ 1 + (1 | Id_Random_Foto)
##    Data: datos_morfo
## 
## REML criterion at convergence: 918.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5035 -0.2443 -0.0051  0.1887  4.3348 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 4.6640   2.1596  
##  Residual                   0.2264   0.4758  
## Number of obs: 336, groups:  Id_Random_Foto, 112
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  14.7535     0.2057 111.0000   71.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtolong_cola)
## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.954
##   Conditional ICC: 0.954

Perímetro_Musculo

mixtoperim_mus<-lmer(Perim_Musculo1~1+(1|Id_Random_Foto),data=datos_morfo)
summary(mixtoperim_mus)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Perim_Musculo1 ~ 1 + (1 | Id_Random_Foto)
##    Data: datos_morfo
## 
## REML criterion at convergence: 1476.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1177 -0.4424  0.0020  0.4175  3.4042 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 22.65    4.759   
##  Residual                    1.25    1.118   
## Number of obs: 336, groups:  Id_Random_Foto, 112
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  32.0149     0.4538 111.0000   70.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtoperim_mus, by_group = TRUE)
## # ICC by Group
## 
## Group          |   ICC
## ----------------------
## Id_Random_Foto | 0.948

Dada la repetibilidad encontrada en las medidas para todas las variables, usé la mediana de cada individuo y con ese valor realicé el resto de los análisis.

Use la mediana porque la mediana no es afectada por los valores extremos y se recomienda mas cuando los valores no tienen distribuución normal.

datos_morfo_med <- read.table("morfologia_mediana_expcasa.txt",h=T)

Figuras exploratorias

explicación pirate_plots

library(ggplot2)
library(viridis)
library(viridisLite)

Altura Cabeza

#para volver una variable categorica:
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)
#para cambiar el orden de los niveles tanto en las figuras como en los analisis y que aparezca primero el tratamiento Agua Estable (AE)
datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AE")
ggplot(datos_morfo_med, aes(x=Postura, y=Altura_Cabeza, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "Altura Cabeza (mm)")

Altura Cola

ggplot(datos_morfo_med, aes(x=Postura, y=Altura_Cola, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "Altura Cola (mm)")

Longitud total

datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AE")
ggplot(datos_morfo_med, aes(x=Postura, y=Longitud_Total, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "Longitud total (mm)")

Longitud Cola

ggplot(datos_morfo_med, aes(x=Postura, y=Longitud_Cola, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "Longitud Cola (mm)")

Perímetro Músculo

ggplot(datos_morfo_med, aes(x=Postura, y=Perimetro_Musculo, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "Perímetro músculo (mm)")

Longitud_Cabeza

ggplot(datos_morfo_med, aes(x=Postura, y=Longitud_Cabeza, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "postura",y = "Longitud Cabeza (mm)")

Test de normalidad y Homogenidad de varianzas

library(tidyverse)
library(rstatix)
library(ggpubr)

Normalidad:

Altura Cabeza
# Density plot
ggdensity(datos_morfo_med$Altura_Cabeza, fill = "lightgray")

# QQ plot
ggqqplot(datos_morfo_med$Altura_Cabeza)

#Shapiro_test

datos_morfo_med %>%group_by(Tratamiento) %>%
  shapiro_test(Altura_Cabeza)
## # A tibble: 3 x 4
##   Tratamiento variable      statistic       p
##   <fct>       <chr>             <dbl>   <dbl>
## 1 AE          Altura_Cabeza     0.963 0.229  
## 2 AB          Altura_Cabeza     0.988 0.958  
## 3 AD          Altura_Cabeza     0.911 0.00458
Longitud_Total
# Density plot
ggdensity(datos_morfo_med$Longitud_Total, fill = "lightgray")

# QQ plot
ggqqplot(datos_morfo_med$Longitud_Total)

#Shapiro_test
datos_morfo_med %>%group_by(Tratamiento) %>%
  shapiro_test(Longitud_Total)
## # A tibble: 3 x 4
##   Tratamiento variable       statistic        p
##   <fct>       <chr>              <dbl>    <dbl>
## 1 AE          Longitud_Total     0.969 0.340   
## 2 AB          Longitud_Total     0.964 0.318   
## 3 AD          Longitud_Total     0.870 0.000345
Altura_Cola
# Density plot
ggdensity(datos_morfo_med$Altura_Cola, fill = "lightgray")

# QQ plot
ggqqplot(datos_morfo_med$Altura_Cola)

#Shapiro_test
datos_morfo_med %>%group_by(Tratamiento) %>%
  shapiro_test(Altura_Cola)
## # A tibble: 3 x 4
##   Tratamiento variable    statistic      p
##   <fct>       <chr>           <dbl>  <dbl>
## 1 AE          Altura_Cola     0.975 0.518 
## 2 AB          Altura_Cola     0.978 0.701 
## 3 AD          Altura_Cola     0.939 0.0365
Longitud_Cola
# Density plot
ggdensity(datos_morfo_med$Longitud_Cola, fill = "lightgray")

# QQ plot
ggqqplot(datos_morfo_med$Longitud_Cola)

#Shapiro_test
datos_morfo_med %>%group_by(Tratamiento) %>%
  shapiro_test(Longitud_Cola)
## # A tibble: 3 x 4
##   Tratamiento variable      statistic        p
##   <fct>       <chr>             <dbl>    <dbl>
## 1 AE          Longitud_Cola     0.979 0.653   
## 2 AB          Longitud_Cola     0.946 0.0948  
## 3 AD          Longitud_Cola     0.870 0.000335
Perímetro_Musculo
# Density plot
ggdensity(datos_morfo_med$Perimetro_Musculo, fill = "lightgray")

# QQ plot
ggqqplot(datos_morfo_med$Perimetro_Musculo)

#Shapiro_test
datos_morfo_med %>%group_by(Tratamiento) %>%
  shapiro_test(Perimetro_Musculo)
## # A tibble: 3 x 4
##   Tratamiento variable          statistic        p
##   <fct>       <chr>                 <dbl>    <dbl>
## 1 AE          Perimetro_Musculo     0.976 0.550   
## 2 AB          Perimetro_Musculo     0.952 0.140   
## 3 AD          Perimetro_Musculo     0.879 0.000564
Longitud_Cabeza
# Density plot
ggdensity(datos_morfo_med$Longitud_Cabeza, fill = "lightgray")

# QQ plot
ggqqplot(datos_morfo_med$Longitud_Cabeza)

#Shapiro_test
datos_morfo_med %>%group_by(Tratamiento) %>%
  shapiro_test(Longitud_Cabeza)
## # A tibble: 3 x 4
##   Tratamiento variable        statistic       p
##   <fct>       <chr>               <dbl>   <dbl>
## 1 AE          Longitud_Cabeza     0.914 0.00585
## 2 AB          Longitud_Cabeza     0.946 0.0905 
## 3 AD          Longitud_Cabeza     0.921 0.00933

Homogenidad de varianza (de cada variable por separado):

library(car)

Levene’s test: A robust alternative to the Bartlett’s test that is less sensitive to departures from normality.

Fligner-Killeen’s test: a non-parametric test which is very robust against departures from normality.

#1. Las variables independientes deben ser categoricas
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)
datos_morfo_med$Postura<- as.factor(datos_morfo_med$Postura)
datos_morfo_med$Caja<- as.factor(datos_morfo_med$Caja)
Longitud_Total
# Levene's test with multiple independent variables
leveneTest( Longitud_Total~Tratamiento*Postura*Caja, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group 20   0.731 0.7845
##       91
Longitud_Cola
# Levene's test with multiple independent variables
leveneTest( Longitud_Cola~Tratamiento, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  2.3386 0.1013
##       109
Longitud_Cabeza
# Levene's test with multiple independent variables
leveneTest( Longitud_Cabeza~Tratamiento, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.8604 0.4259
##       109
Altura_Cabeza
# Levene's test with multiple independent variables
leveneTest( Altura_Cabeza~Tratamiento, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2   1.276 0.2833
##       109
Altura_Cola
# Levene's test with multiple independent variables
leveneTest( Altura_Cola~Tratamiento, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   2  4.0223 0.02064 *
##       109                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Perímetro_Músculo
# Levene's test with multiple independent variables
leveneTest(Perimetro_Musculo~Tratamiento, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  2.2158  0.114
##       109

Análisis de componentes pricipales

Correlación entre variables

morfo<-datos_morfo_med[,8:12] #primero los individuos activos y luego las variables
head(morfo)
##   Altura_Cabeza Altura_Cola Longitud_Cola Perimetro_Musculo Longitud_Cabeza
## 1         2.145       1.691        17.326            37.880           6.600
## 2         2.536       1.336        14.029            30.549           7.144
## 3         2.328       1.495        17.072            37.172           8.039
## 4         2.172       1.760        17.075            37.207           6.721
## 5         2.062       1.496        15.598            34.448           6.555
## 6         2.283       1.455        15.270            33.554           7.736
La Correlación
corrl<-round(cor(morfo),3)
Las figuras de las correlaciones con coeficientes
library(corrplot)
library(ggcorrplot)

ggcorrplot(corrl, hc.order = TRUE, type = "lower",
           lab = TRUE)

El PCA

Los paquetes:

para hacer el test de barlett y el analisis de componentes: library(psych)

para graficar que variables pesan mas en que PC library(corrplot)

for ggplot2-based visualization: library(factoextra)

library(factoextra)
library(FactoMineR)
library(psych)
El test de bartlett:

Ho:the sample correlation matrix came from a multivariate normal population in which the variables of interest are independent.

Rejection of the hypothesis is taken as an indication that the data are appropriate for analysis.

(CHARLES D. DZIUBAN, 1974) Bartlett's test of sphericity tests the hypothesis that your correlation matrix is an identity matrix,which would indicate that your variables are unrelated and therefore unsuitable for structure detection. Small values (less than 0.05) of the significance level indicate that a factor analysis may be useful with your data.

cortest.bartlett(corrl,n=112)
## $chisq
## [1] 878.6884
## 
## $p.value
## [1] 2.45595e-182
## 
## $df
## [1] 10

KMO: Test alternativo al test de barlett. Measure of Sampling Adequacy

This is just a function of the squared elements of the ‘image’ matrix compared to the squares of the original correlations. The overall MSA as well as estimates for each item are found.The index is known as the Kaiser-Meyer-Olkin (KMO) index.Kaiser-Meyer-Olkin (KMO). KMO Test is a measure of how suited your data is for Factor Analysis.

KMO(morfo)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = morfo)
## Overall MSA =  0.78
## MSA for each item = 
##     Altura_Cabeza       Altura_Cola     Longitud_Cola Perimetro_Musculo 
##              0.90              0.94              0.70              0.69 
##   Longitud_Cabeza 
##              0.78
Componentes:

componentes factorMineR

#PCA usando factorMineR y factorextra
elpca<-PCA(morfo, scale.unit = TRUE,ncp=5,graph = FALSE)
summary(elpca)
## 
## Call:
## PCA(X = morfo, scale.unit = TRUE, ncp = 5, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5
## Variance               4.229   0.312   0.255   0.198   0.005
## % of var.             84.586   6.247   5.107   3.962   0.098
## Cumulative % of var.  84.586  90.833  95.940  99.902 100.000
## 
## Individuals (the 10 first)
##                       Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## 1                 |  2.206 |  1.971  0.820  0.798 | -0.757  1.636  0.118 |
## 2                 |  1.860 |  0.776  0.127  0.174 |  1.574  7.079  0.716 |
## 3                 |  2.564 |  2.403  1.220  0.879 |  0.106  0.032  0.002 |
## 4                 |  2.284 |  2.060  0.896  0.814 | -0.507  0.736  0.049 |
## 5                 |  0.933 |  0.810  0.139  0.754 | -0.293  0.245  0.099 |
## 6                 |  1.753 |  1.413  0.421  0.649 |  0.679  1.318  0.150 |
## 7                 |  1.174 | -1.002  0.212  0.728 |  0.529  0.799  0.203 |
## 8                 |  1.762 | -1.526  0.492  0.750 |  0.558  0.891  0.100 |
## 9                 |  2.020 | -1.682  0.597  0.693 | -0.517  0.764  0.065 |
## 10                |  0.694 |  0.090  0.002  0.017 |  0.218  0.136  0.098 |
##                    Dim.3    ctr   cos2  
## 1                 -0.477  0.796  0.047 |
## 2                  0.287  0.288  0.024 |
## 3                  0.883  2.724  0.118 |
## 4                 -0.549  1.054  0.058 |
## 5                 -0.050  0.009  0.003 |
## 6                  0.757  2.005  0.187 |
## 7                  0.018  0.001  0.000 |
## 8                 -0.630  1.386  0.128 |
## 9                 -0.580  1.176  0.082 |
## 10                 0.538  1.013  0.602 |
## 
## Variables
##                      Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Altura_Cabeza     |  0.888 18.640  0.788 |  0.356 40.471  0.126 | -0.138  7.483
## Altura_Cola       |  0.907 19.451  0.823 |  0.066  1.399  0.004 | -0.270 28.445
## Longitud_Cola     |  0.948 21.241  0.898 | -0.296 28.052  0.088 | -0.017  0.110
## Perimetro_Musculo |  0.963 21.909  0.927 | -0.255 20.896  0.065 |  0.025  0.239
## Longitud_Cabeza   |  0.891 18.759  0.793 |  0.169  9.182  0.029 |  0.403 63.722
##                     cos2  
## Altura_Cabeza      0.019 |
## Altura_Cola        0.073 |
## Longitud_Cola      0.000 |
## Perimetro_Musculo  0.001 |
## Longitud_Cabeza    0.163 |
#diagrama de codos
fviz_eig(elpca, addlabels = TRUE, ylim = c(0, 90))

#proporcion de varianza explicada:
eig.val <- get_eigenvalue(elpca)
variables <- get_pca_var(elpca)
# Contribucion de las variables en los componentes
corrplot(variables$contrib, is.corr=FALSE)

#el biplot de las cargas
fviz_pca_var (elpca, col.var = "contrib",
             gradient.cols=c("#481567FF","#287D8EFF","#3CBB75FF","#B8DE29FF"))

#Para exportar los componentes:
elpca$ind$coord
##            Dim.1        Dim.2        Dim.3        Dim.4         Dim.5
## 1    1.971417695 -0.756638295 -0.477155484  0.425086834 -0.0279801995
## 2    0.776411577  1.573712458  0.287113650 -0.543131292  0.0600948479
## 3    2.403476531  0.106129259  0.882666316 -0.014853194  0.0911779889
## 4    2.060087543 -0.507492988 -0.549039404  0.641435481  0.0101047462
## 5    0.809968198 -0.292772133 -0.049787470  0.349663487 -0.0554580796
## 6    1.412549803  0.679062245  0.757310497  0.197820815  0.0583394330
## 7   -1.001847838  0.528747467  0.017734936  0.295198745  0.0903691309
## 8   -1.526416716  0.558376344 -0.629520809 -0.254413120 -0.0567935325
## 9   -1.681836349 -0.516914621 -0.579862481  0.799464104  0.0919674147
## 10   0.089849974  0.217772823  0.538310797  0.363413391 -0.0662889709
## 11  -1.660479964  0.289539215  0.012022645  0.145250817 -0.0223885519
## 12  -2.660442862 -0.006239710 -0.426766861  0.540413636  0.0958374070
## 13   1.373499949 -0.315088373 -0.314195325  0.180963484  0.0631165831
## 14   1.159377704  0.165097563  0.823658320  0.128768371  0.0140157663
## 15  -2.136215549 -0.294840766  0.831254899  0.201850397  0.0021528130
## 16  -2.567155075  0.219855836 -0.321042316  0.345161847  0.0241900031
## 17  -1.794900692 -0.114111457  0.259753145 -0.623526054  0.0403628242
## 18  -1.356177464  0.377478335  0.501092695  0.009155788 -0.1721827107
## 19  -1.296125970  0.677467421  0.131762731 -0.230745496 -0.0342607902
## 20  -1.865307155 -0.146174965 -0.115963634  0.061750280 -0.0075188623
## 21  -1.006750744  0.122629602 -0.843374396  0.209026101 -0.0367531138
## 22  -3.136540620 -0.410808011 -0.195444725  0.227814977 -0.0557286059
## 23  -1.115751821  0.061017430 -0.160573055  0.382321725 -0.0192236391
## 24  -1.315858103  0.009746118 -0.617147244  1.004327922 -0.0329013975
## 25   0.246590956 -0.388084768  0.032405802 -0.022029690 -0.0917878521
## 26  -0.682192711 -0.206541787 -0.483463675 -0.104169641 -0.0521253573
## 27   0.333303945  0.021094277 -0.656586622 -0.732984488  0.0361807617
## 28  -0.825794755  0.459082631 -0.056026789 -0.585748757  0.0274731719
## 29   0.246961204 -0.251435221 -0.425781169 -0.017339584 -0.0213978448
## 30  -0.340889610 -0.184954331 -0.172507073  0.381872611 -0.0247863233
## 31  -0.606315697 -0.481907969  0.152866869 -0.775756978  0.0017491764
## 32  -2.814582666 -0.868631220  0.509526334  0.070025468  0.0534490304
## 33  -0.692635706 -0.304125670  0.307897536 -0.253250426  0.0448202141
## 34  -1.197296806 -0.690714672  0.011562363 -0.513302560  0.0217996708
## 35   6.575672257 -0.523907217  0.103080049  0.126498394 -0.0360533837
## 36   5.990249900 -0.188123888  0.344220080 -0.265288518 -0.0276080823
## 37   5.920061399  0.633904787 -0.213945615 -0.693376762  0.0413358859
## 38   4.888903774 -1.719155284  0.499414097  0.430574644 -0.0182586406
## 39   5.298095256 -0.587405212  0.246542108  0.775188733 -0.0399952552
## 40   0.634282010  0.694877732 -0.323334754  0.390829993  0.0246253329
## 41   0.194177348  0.522868652  1.070601318  0.262214596 -0.0292361990
## 42   0.263782478  0.049792034  0.061453047  0.462375214  0.0039599784
## 43   3.843344905 -0.278135989  0.169225345  0.627385755 -0.1525049509
## 44  -0.067951299 -0.148759730  0.092513457  0.079162186 -0.1703096556
## 45   3.664574327 -0.146571570 -0.197617842  0.449982936  0.1046614101
## 46   0.070010178  0.503117300 -0.204119626 -0.248884896 -0.0820017391
## 47  -0.454933471  0.068024046  0.317264970 -0.062388412 -0.0940471766
## 48   1.189943732 -0.003987709  0.645943141  0.146891759  0.0742235851
## 49  -1.672084859  0.495061991  0.038307066 -0.080526709 -0.0206963304
## 50  -1.889611725 -0.014874418  0.095516521  0.058419825  0.0103230571
## 51  -0.738472208  0.809721340  0.276743998  0.057812536  0.0137373916
## 52  -1.561626703 -0.071289969  0.542856054 -0.479753229  0.1081942157
## 53   1.055908658  0.860868567 -0.079701476 -0.273030149 -0.0409074256
## 54  -0.353881742  0.070177407  0.634408810 -1.057765233  0.0149369271
## 55  -1.572948042 -0.252398279  0.058324734  0.283004646  0.0480392282
## 56  -2.507591602 -0.269524549 -0.084943224  0.070521440  0.0393038902
## 57   0.974636096  0.015143203 -0.007664674  0.486292822  0.0423990716
## 58   0.250874705  0.789104736 -0.691534051  1.162277335  0.0404221838
## 59   0.363605833  0.922335370 -0.294000911  0.050162213  0.1065021167
## 60  -0.559301927 -0.138729736 -0.110464954 -0.058472376  0.0351539779
## 61  -0.909764310 -0.031322655 -0.196571479  0.251313442  0.0640675960
## 62   4.460114354  0.681696256 -0.943110173 -0.584174814  0.0872465948
## 63  -1.651801171 -0.130263891 -0.137721960 -0.036433030  0.0003102448
## 64   2.795022463  1.099105533 -1.490456706 -0.524438451 -0.0182565312
## 65   0.125804347 -0.296874033 -0.853862039 -0.230866257  0.0487989848
## 66   0.000646991  0.502282247 -0.756089477 -0.393528050 -0.1530378720
## 67   2.013197183 -0.808519601 -0.280846155  0.343665215  0.0916229142
## 68  -0.714605878 -0.133140808  0.139807336 -0.096776805 -0.0997726168
## 69   0.204159702 -0.675505478  0.135968182 -0.171570685  0.0056693906
## 70   1.206889268 -0.126527164 -0.285848094  0.471607380 -0.1339756881
## 71   0.370627540 -0.241902155  0.229484909  0.362824779  0.0013195684
## 72   1.006874465 -2.020720815 -0.692974067 -1.551159896 -0.0033772908
## 73  -2.299001151 -0.611615759  0.064683156 -0.142865898 -0.0527659421
## 74  -0.489823975 -0.145575326  0.231158218  0.466605956 -0.0330037038
## 75  -0.586355906 -0.469806835 -0.035094201 -0.044408968  0.1148953184
## 76   0.539843874  0.068338620 -0.204569968  0.481107160  0.0765827700
## 77   2.709508445  0.184390724  1.069613535 -0.111904465 -0.0057425414
## 78  -0.598517858 -0.684886992  0.736401669  0.122906836 -0.0301604617
## 79  -1.129425854 -0.037472693  0.192071135 -0.036860904 -0.0332344697
## 80   0.917674127 -0.358833541  0.263462207 -0.170367757 -0.0213321057
## 81   0.983350472 -0.249452885 -0.101729528 -0.249546171  0.0627365277
## 82  -0.167196174  0.217620633  0.554271618  0.071258067  0.1132321538
## 83   2.800364565 -0.224901786  0.293241647 -0.012314798  0.0550320120
## 84  -0.273586196 -0.964025141  0.078468337  0.066111289 -0.2048111185
## 85  -1.705268590  0.356964952 -0.004006259  0.053282812 -0.0017568897
## 86   2.374316024  0.461937396 -0.218866603 -0.562354480  0.1171307333
## 87  -1.343067854 -0.458147551  0.189699333 -0.013221570  0.0025875708
## 88  -1.144739243  0.967704595 -0.008880218 -0.113907232 -0.1012491473
## 89  -1.877461334 -0.210902786  0.017165947  0.110546878 -0.0082713640
## 90  -2.593741321 -0.229760016 -0.122828693  0.257701530 -0.0861822280
## 91   1.956269818  1.116730685  0.960537377 -0.381279830 -0.0553624675
## 92  -3.203217727  0.421148561  0.413142302 -0.081277949  0.0183839418
## 93  -1.395996061  0.330002568 -0.527620756  0.429268955  0.1234390749
## 94   1.411369295 -0.249883429 -0.502697795 -1.287431112 -0.0040423337
## 95  -1.672707825 -0.276646955 -0.416222053  0.187411870 -0.0327418462
## 96  -1.917631974 -0.712841291  0.093844433 -0.309530909  0.1668272506
## 97   1.629779123 -0.045574773 -0.080762730 -0.457076843 -0.0543352923
## 98   4.155757724  0.411220513  0.701571229 -0.138739472 -0.0014964822
## 99  -0.507373643  0.072407598  0.363520497 -0.122410253  0.0410483137
## 100 -1.278327292 -0.573151884 -0.004700194 -0.114299724 -0.0547009274
## 101  0.984420372  0.732797023  0.940067209  0.667525080 -0.0046767981
## 102 -3.307648339 -0.134111957  0.249885369  0.459839642  0.0087456766
## 103 -0.487073553 -0.770524486 -0.062915883 -0.115621706 -0.0171155008
## 104 -1.305656706  0.426465859  0.002025776 -0.678778164 -0.0755041334
## 105 -1.675392502  0.252760324 -0.368220754 -0.394825081 -0.0653393514
## 106 -0.994943073  1.467362433  1.065372970  0.241506102 -0.0337639255
## 107 -0.602836793  0.481643298 -0.936321113 -0.258329955 -0.1340707283
## 108 -2.288242348  0.025304813  0.196100489 -0.318132177 -0.0048696490
## 109 -3.414592500  0.121109282 -0.066535010 -0.751151590  0.0555402422
## 110  0.221775510  0.051907734 -0.697796429 -0.459406673  0.0119945129
## 111  1.400806902  0.644047640 -1.803393057  0.966969009  0.0581945396
## 112 -2.134274901 -0.643522282  0.647287842 -0.086104076  0.1537189091
# Export into a TXT file
write.infile(elpca$ind$coord, "pca.txt", sep = "\t")

# Export into a CSV file
write.infile(elpca$ind$coord, "pca.csv", sep = ";")
Rotando los componentes

Abdi & Williams, 2010, Principal component analysis: After the number of components has been determined, and in order to facilitate the interpretation,the analysis often involves a rotation of the components that were retained. principal: Standardized loadings (pattern matrix) based upon correlation matrix

morforot <- principal(morfo,1,rotate="varimax")
morforot1<-principal(morfo, rotate="varimax", nfactors=1, scores=TRUE)
cargas_rot<-morforot1$loadings

Para extraer los componentes ya rotados

componentes_rot<-morforot1$scores[1:112,]

write.csv(componentes_rot,file ="morfo_rot1.csv")

Figuras con el PC para ver el efecto del tratamiento

datos_morfo_med <- read.table("morfologia_mediana_expcasa.txt",h=T)
head(datos_morfo_med)
##   Id_Random_Foto Tratamiento Postura Tra_Pos Caja Individuo Longitud_Total
## 1             23          AB      P1    ABP1    1         1         24.217
## 2            112          AB      P1    ABP1    1         2         21.001
## 3             27          AB      P1    ABP1    1         3         25.013
## 4             53          AB      P1    ABP1    1         4         23.821
## 5             28          AB      P1    ABP1    1         5         22.316
## 6             51          AB      P2    ABP2    2         1         23.070
##   Altura_Cabeza Altura_Cola Longitud_Cola Perimetro_Musculo Longitud_Cabeza
## 1         2.145       1.691        17.326            37.880           6.600
## 2         2.536       1.336        14.029            30.549           7.144
## 3         2.328       1.495        17.072            37.172           8.039
## 4         2.172       1.760        17.075            37.207           6.721
## 5         2.062       1.496        15.598            34.448           6.555
## 6         2.283       1.455        15.270            33.554           7.736
##     PC1psych     Pcrot    PC1     PC2   Dim1    Dim2
## 1 -2.2285871 0.9543284 2.2386 -0.6195 1.9714 -0.7566
## 2 -0.6921240 0.3758471 0.6952  1.5081 0.7764  1.5737
## 3 -2.7291243 1.1634804 2.7414 -0.0910 2.4035  0.1061
## 4 -2.2501256 0.9972519 2.2602 -0.3228 2.0601 -0.5075
## 5 -0.9182887 0.3920913 0.9224 -0.2631 0.8100 -0.2928
## 6 -1.5618954 0.6837903 1.5689  0.5065 1.4125  0.6791
#para volver una variable categorica:
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)
#para cambiar el orden de los niveles tanto en las figuras como en los analisis y que aparezca primero el tratamiento Agua Estable (AE)

datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AE")

ggplot(datos_morfo_med, aes(x=Postura, y=Dim1, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "PC1 (Tamaño Corporal)")

Figuras con el PCrotado para ver el efecto del tratamiento

ggplot(datos_morfo_med, aes(x=Postura, y=Pcrot, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "PC Rotado(Tamaño Corporal)")

El modelo mixto para evaluar el efecto del tratamiento de desecación:

1.Modelo mixto con PC1

mixtopc1<-lmer(Dim1~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixtopc1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dim1 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 438.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1087 -0.7158 -0.1341  0.4910  2.7505 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.707    0.8408  
##  Postura  (Intercept) 1.039    1.0192  
##  Residual             2.359    1.5359  
## Number of obs: 112, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)    -0.4593     0.5570 12.1157  -0.825   0.4255  
## TratamientoAB  -0.1445     0.5779 12.6247  -0.250   0.8066  
## TratamientoAD   1.4615     0.5689 11.9002   2.569   0.0247 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.503       
## TratamintAD -0.511  0.492

Aqui otra forma de presentar los resultados:

Como crear una tabla con los resultados de los análisis de modelos mixtos

*Daniel Lüdecke

The marginal R-squared considers only the variance of the fixed effects, while the conditional R-squared takes both the fixed and random effects into account.

The p-value is a simple approximation, based on the t-statistics and using the normal distribution function. A more precise p-value can be computed using p.val = "kr". In this case, which only applies to linear mixed models, the computation of p-values is based on conditional F-tests with Kenward-Roger approximation for the degrees of freedom (using the using the pbkrtest-package). Note that here the computation is more time consuming and thus not used as default. You can also display the approximated degrees of freedom with show.df.

Para calcular el R*2 del modelo

R*2 en modelos mixtos

R* en modelos mixtos 2

Nakagawa S, Johnson P, Schielzeth H (2017) The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisted and expanded. J. R. Soc. Interface 14. doi: 10.1098/rsif.2017.0213.

library(sjPlot)
library(lme4)

tab_model(mixtopc1)
  Dim1
Predictors Estimates CI p
(Intercept) -0.46 -1.55 – 0.63 0.410
Tratamiento [AB] -0.14 -1.28 – 0.99 0.803
Tratamiento [AD] 1.46 0.35 – 2.58 0.010
Random Effects
σ2 2.36
τ00 Caja 0.71
τ00 Postura 1.04
ICC 0.43
N Caja 21
N Postura 7
Observations 112
Marginal R2 / Conditional R2 0.116 / 0.492 
tab_model(mixtopc1, p.val = "kr", show.df = TRUE)
  Dim1
Predictors Estimates CI p df
(Intercept) -0.46 -1.67 – 0.75 0.426 12.17
Tratamiento [AB] -0.14 -1.40 – 1.11 0.807 12.21
Tratamiento [AD] 1.46 0.22 – 2.71 0.025 11.51
Random Effects
σ2 2.36
τ00 Caja 0.71
τ00 Postura 1.04
ICC 0.43
N Caja 21
N Postura 7
Observations 112
Marginal R2 / Conditional R2 0.116 / 0.492

calculo de eta-squared

library(effectsize)
#proporcion de varianza atribuida a una variable

eta_squared(mixtopc1)
## Parameter   | Eta2 (partial) |       90% CI
## -------------------------------------------
## Tratamiento |           0.44 | [0.04, 0.66]
anova(mixtopc1)
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Tratamiento  22.71  11.355     2 12.362  4.8133 0.02845 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Prueba de homogeneidad del modelo mixto pc1
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(mixtopc1)

trat<-lm(Dim1~Tratamiento,data=datos_morfo_med)
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(trat)

Breusch Pagan Test

Breusch Pagan Test was introduced by Trevor Breusch and Adrian Pagan in 1979. It is used to test for heteroskedasticity in a linear regression model and assumes that the error terms are normally distributed. It tests whether the variance of the errors from a regression is dependent on the values of the independent variables. It is a χ2 test.

You can perform the test using the fitted values of the model, the predictors in the model and a subset of the independent variables. It includes options to perform multiple tests and p value adjustments. The options for p value adjustments include Bonferroni, Sidak and Holm’s method.

library(olsrr) #para el analisis de homog. var de breusch_pagan. No funciona si el modelo no es lm

datos_morfo_med <- read.table("morfologia_mediana_expcasa.txt",h=T)


#hice una columna combinando tratamiento y postura. En esa variable Tra_Pos quedaron 21 categorias.

tra_pos<-lm(Dim1~Tra_Pos, data=datos_morfo_med)
plot(fitted(mixtopc1),residuals(mixtopc1))

ols_test_breusch_pagan(tra_pos)
## 
##  Breusch Pagan Test for Heteroskedasticity
##  -----------------------------------------
##  Ho: the variance is constant            
##  Ha: the variance is not constant        
## 
##               Data               
##  --------------------------------
##  Response : Dim1 
##  Variables: fitted values of Dim1 
## 
##         Test Summary         
##  ----------------------------
##  DF            =    1 
##  Chi2          =    0.2599549 
##  Prob > Chi2   =    0.6101511
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(tra_pos)

Otra forma de evaluar homogeneidad:

Otra forma de evaluar la homogeneidad propuesta por Michael Palmeri: evaluando homocedasticidad en modelos mixtos

library(lmerTest)
library(lme4)
#extracts the residuals and places them in a new column in our original data table
datos_morfo_med$mixtopc1.Res<-residuals(mixtopc1)
#creates a new column with the absolute value of the residuals
datos_morfo_med$Abs.mixtopc1.Res<-abs(datos_morfo_med$mixtopc1.Res)
#squares the absolute values of the residuals to provide the more robust estimate

datos_morfo_med$mixtopc1.Res2 <- datos_morfo_med$Abs.mixtopc1.Res^2 
#ANOVA of the squared residuals
Levene.mixtopc1 <- lm(mixtopc1.Res2 ~Caja+Postura, data=datos_morfo_med) 

anova(Levene.mixtopc1) #displays the results
## Analysis of Variance Table
## 
## Response: mixtopc1.Res2
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Caja        1  47.31  47.312  5.8087 0.0177 *
## Postura     6  71.72  11.953  1.4675 0.1966  
## Residuals 104 847.09   8.145                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para revisar los supuestos en modelos mixtos:

evaluando los supuestos en modelos mixtos

Francisco Juretig. “R Statistics Cookbook.”: Mixed models can be impacted greatly by outliers. Even a minor contamination causes major estimation problems. Solución: evaluar que tan robusto es el modelo.

#para volver una variable categorica:
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)

Re-estimando el modelo con "robust errors"

Based on central model approach: That is, we assume the model to be true, but a part of the data to possibly be contaminated. Irrespective of this contamination, we would like to estimate the parameters that define the central model and these estimates should be only minimally influenced by the contamination. A robust estimation method should give reasonable results in the presence of such contaminated error distributions, while still maintain- ing efficiency in case there is no contamination.

library(robustlmm)

mixtopc1_robust <- rlmer(Dim1~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixtopc1_robust)
## Robust linear mixed model fit by DAStau 
## Formula: Dim1 ~ Tratamiento + (1 | Caja) + (1 | Postura) 
##    Data: datos_morfo_med 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.58392 -0.60911 -0.08615  0.59403  3.12392 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.0000   0.0000  
##  Postura  (Intercept) 0.6729   0.8203  
##  Residual             2.7695   1.6642  
## Number of obs: 112, groups: Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)   -0.80174    0.43293  -1.852
## TratamientoAD  1.48947    0.40161   3.709
## TratamientoAE  0.02351    0.40161   0.059
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAD
## TratamintAD -0.497       
## TratamintAE -0.497  0.536
## 
## Robustness weights for the residuals: 
##  88 weights are ~= 1. The remaining 24 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.431   0.725   0.843   0.808   0.904   0.998 
## 
## Robustness weights for the random effects: 
##  27 weights are ~= 1. The remaining one are
##    22 
## 0.509 
## 
## Rho functions used for fitting:
##   Residuals:
##     eff: smoothed Huber (k = 1.345, s = 10) 
##     sig: smoothed Huber, Proposal II (k = 1.345, s = 10) 
##   Random Effects, variance component 1 (Caja):
##     eff: smoothed Huber (k = 1.345, s = 10) 
##     vcp: smoothed Huber, Proposal II (k = 1.345, s = 10) 
##   Random Effects, variance component 2 (Postura):
##     eff: smoothed Huber (k = 1.345, s = 10) 
##     vcp: smoothed Huber, Proposal II (k = 1.345, s = 10)
library(sjPlot)
library(lme4)
tab_model(mixtopc1_robust)
  Dim1
Predictors Estimates CI p
(Intercept) -0.80 -1.65 – 0.05 0.064
Tratamiento [AD] 1.49 0.70 – 2.28 <0.001
Tratamiento [AE] 0.02 -0.76 – 0.81 0.953
Random Effects
σ2 2.77
τ00 Caja 0.00
τ00 Postura 0.67
ICC 0.20
N Caja 21
N Postura 7
Observations 112
Marginal R2 / Conditional R2 0.127 / 0.297 
plot(mixtopc1_robust)

Evaluando los residuos del modelo

Excerpt From: Francisco Juretig. “R Statistics Cookbook.” Apple Books.

“In the fitted-residual relationship; ideally, there should be no structure there. We then plotted boxplots for the residuals in each group.

Excerpt From: Francisco Juretig. “R Statistics Cookbook.” Apple Books.

Efecto aleatorio:La postura
plot(mixtopc1, resid(., scaled=TRUE) ~ fitted(.) |Postura, abline = 0)

#The following output shows the scaled residuals for each group:
plot(mixtopc1, Postura ~ resid(., scaled=TRUE))

“Now, let's plot the fitted versus the actuals. There should be a positive linear relationship here. If this weren't the case, it would imply that we were missing some structure in the model:”

plot (mixtopc1, Dim1 ~ fitted(.) | Postura, abline = c (0,1))

Efecto aleatorio: La caja
plot(mixtopc1, resid(., scaled=TRUE) ~ fitted(.) |Caja,abline = 0)

#The following output shows the scaled residuals for each group:
plot(mixtopc1, Caja ~ resid(., scaled=TRUE)) 

plot (mixtopc1, Dim1 ~ fitted(.) | Caja, abline = c (0,1))

3.La magnitud del efecto:

Effect sizes, are metrics that represent the amount of differences between two sample means.

Cohen’s d and Hedges’ g are interpreted in a similar way. Cohen suggested using the following rule of thumb for interpreting results:

Small effect (cannot be discerned by the naked eye) = 0.2 Medium Effect = 0.5 Large Effect (can be seen by the naked eye) = 0.8

Cohen, J. (1977). Statistical power analysis for the behavioral sciences. Routledge.

library(emmeans)
trat.emm.s <- emmeans(mixtopc1, "Tratamiento")
pairs(trat.emm.s)
##  contrast estimate    SE   df t.ratio p.value
##  AE - AB     0.144 0.578 12.2  0.250  0.9663 
##  AE - AD    -1.461 0.569 11.5 -2.568  0.0609 
##  AB - AD    -1.606 0.578 12.2 -2.777  0.0407 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
eff_size(trat.emm.s, sigma = sigma(mixtopc1), edf = 12)
##  contrast effect.size    SE   df lower.CL upper.CL
##  AE - AB       0.0941 0.377 12.2   -0.726    0.914
##  AE - AD      -0.9515 0.377 12.2   -1.773   -0.130
##  AB - AD      -1.0456 0.380 12.2   -1.873   -0.218
## 
## sigma used for effect sizes: 1.536 
## Degrees-of-freedom method: inherited from kenward-roger when re-gridding 
## Confidence level used: 0.95

Effect size in mixed models

Effect size in mixed models1

Plotting Random Effects of Mixed Models

Visualizing (generalized) linear mixed effects models

AB y AD
library(sjPlot)
library(lme4)

#para volver una variable categorica:
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)

datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AE")
mixtopc1<-lmer(Dim1~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
plot_model(mixtopc1)

AE y AD
#para volver una variable categorica:
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)

datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AB")
mixtopc_AB<-lmer(Dim1~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixtopc_AB)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dim1 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 438.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1087 -0.7158 -0.1341  0.4910  2.7505 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.707    0.8408  
##  Postura  (Intercept) 1.039    1.0192  
##  Residual             2.359    1.5359  
## Number of obs: 112, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)    -0.6038     0.5661 12.8986  -1.067    0.306  
## TratamientoAE   0.1445     0.5779 12.6247   0.250    0.807  
## TratamientoAD   1.6060     0.5779 12.6247   2.779    0.016 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAE
## TratamintAE -0.526       
## TratamintAD -0.526  0.515
plot_model(mixtopc_AB)

AE vs AD
#para volver una variable categorica:
datos_morfo_med$Tratamiento <- as.factor(datos_morfo_med$Tratamiento)

datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AD")
mixtopc_AD<-lmer(Dim1~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixtopc_AD)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dim1 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 438.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1087 -0.7158 -0.1341  0.4910  2.7505 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.707    0.8408  
##  Postura  (Intercept) 1.039    1.0192  
##  Residual             2.359    1.5359  
## Number of obs: 112, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)     1.0022     0.5570 12.1157   1.799   0.0969 .
## TratamientoAB  -1.6060     0.5779 12.6247  -2.779   0.0160 *
## TratamientoAE  -1.4615     0.5689 11.9002  -2.569   0.0247 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.503       
## TratamintAE -0.511  0.492
plot_model(mixtopc_AD)

Probando otras opciones para ver magnitud del efecto:

library(effectsize)
effectsize(mixtopc1)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |              -0.22 | [-0.75, 0.31]
## TratamientoAB |              -0.07 | [-0.62, 0.48]
## TratamientoAD |               0.71 | [ 0.17, 1.25]
## 
## # Standardization method: refit
effectsize(mixtopc_AB)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |              -0.29 | [-0.83, 0.24]
## TratamientoAE |               0.07 | [-0.48, 0.62]
## TratamientoAD |               0.78 | [ 0.23, 1.33]
## 
## # Standardization method: refit
effectsize(mixtopc_AD)
## Parameter     | Coefficient (std.) |         95% CI
## ---------------------------------------------------
## (Intercept)   |               0.49 | [-0.04,  1.01]
## TratamientoAB |              -0.78 | [-1.33, -0.23]
## TratamientoAE |              -0.71 | [-1.25, -0.17]
## 
## # Standardization method: refit

4.El modelo mixto con PC rotado

mixto1pc1<-lmer(Pcrot~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixto1pc1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pcrot ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 280.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1087 -0.7158 -0.1341  0.4910  2.7505 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.1657   0.4070  
##  Postura  (Intercept) 0.2434   0.4934  
##  Residual             0.5528   0.7435  
## Number of obs: 112, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)     0.4851     0.2696 12.1157   1.799   0.0969 .
## TratamientoAB  -0.7774     0.2798 12.6247  -2.779   0.0160 *
## TratamientoAE  -0.7075     0.2754 11.9002  -2.569   0.0247 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.503       
## TratamintAE -0.511  0.492
library(sjPlot)
library(lme4)
tab_model(mixto1pc1)
  Pcrot
Predictors Estimates CI p
(Intercept) 0.49 -0.04 – 1.01 0.072
Tratamiento [AB] -0.78 -1.33 – -0.23 0.005
Tratamiento [AE] -0.71 -1.25 – -0.17 0.010
Random Effects
σ2 0.55
τ00 Caja 0.17
τ00 Postura 0.24
ICC 0.43
N Caja 21
N Postura 7
Observations 112
Marginal R2 / Conditional R2 0.116 / 0.492 
anova(mixto1pc1)
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Tratamiento 5.3219  2.6609     2 12.362  4.8133 0.02845 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(mixto1pc1)

5. El modelo mixto con PC2

mixto1pc2<-lmer(Dim2~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)

summary(mixto1pc2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dim2 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 192.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1792 -0.5581 -0.0878  0.5469  3.1387 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.008841 0.09402 
##  Postura  (Intercept) 0.027644 0.16626 
##  Residual             0.288700 0.53731 
## Number of obs: 112, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)   -0.02753    0.11240 14.18087  -0.245    0.810
## TratamientoAB  0.01193    0.13609 11.73078   0.088    0.932
## TratamientoAE  0.06033    0.13176 10.42076   0.458    0.656
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.568       
## TratamintAE -0.586  0.484
library(sjPlot)
library(lme4)
tab_model(mixto1pc2)
  Dim2
Predictors Estimates CI p
(Intercept) -0.03 -0.25 – 0.19 0.806
Tratamiento [AB] 0.01 -0.25 – 0.28 0.930
Tratamiento [AE] 0.06 -0.20 – 0.32 0.647
Random Effects
σ2 0.29
τ00 Caja 0.01
τ00 Postura 0.03
ICC 0.11
N Caja 21
N Postura 7
Observations 112
Marginal R2 / Conditional R2 0.002 / 0.114 
anova(mixto1pc2)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## Tratamiento 0.067306 0.033653     2 11.217  0.1166  0.891
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(mixto1pc2)

datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AE")
ggplot(datos_morfo_med, aes(x=Postura, y=Dim2, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "Postura",y = "PC2 (Altura cabeza)")

Resultados Desempeño Locomotor

Para este experimento se utilizaron 7 posturas diferentes. Cada postura se dividió en 3 tratamientos (Agua Baja AB,Agua Estable AE ,Agua Desecación AD), para un total de 21 cajas.5 individuos provenientes de cada caja fueron utilizados para las pruebas de desempeño locomotor. Estas pruebas se realizaron entre las 9 y las 16 horas del día. Se evaluó el desempeño locomotor a 25, 30, 33, 37 y 40 grados centígrados. El desempeño de cada renacuajo fue evaluado a una única temperatura y para esto el animal fue introducido en recipiente con agua durante 3 minutos y se filmó su respuesta al ser estimulado con un estilete. A partir de los videos y utlizando el programa Tracker se determinó para cada renacuajo los valores máximos de velocidad, aceleración y distancia.

programa para análisis de videos

knitr::include_graphics("locom.png")

Los datos:

desloc<- read.table("desloc_def.txt",h=T)
head(desloc)
##   tratamiento temp replica       vel     acel      long     vel_a      mag_p
## 1          AB 25.5      R1 0.3994130 4.304755 1.0157336 113.82766 0.14997262
## 2          AB 25.5      R2 0.3856361 4.718915 0.8177761 232.17318 0.14841826
## 3          AB 25.5      R3 0.3862476 3.794305 0.1733697 181.24012 0.13290586
## 4          AB 25.5      R4 0.2625172 2.099849 0.1875891 222.44808 0.14688139
## 5          AB 25.5      R5 0.2067526 2.053904 0.1393907  13.37304 0.09035338
## 6          AB 25.5      R6 0.3568749 2.757037 0.2041926 356.23589 0.15353251

Librerias para las figuras:

library(ggplot2)
library(yarrr)
library(ggstance)
library(ggformula)
library(viridisLite)
library(viridis)

Figuras exploratorias

Velocidad
#velocidad

ggplot(data = desloc,aes(x = temp, y =vel, group=replica))+
  facet_grid(.~tratamiento)+theme_bw()+geom_spline(aes(x=temp, y=vel,color=replica))+
  coord_cartesian(xlim = c(25,40),ylim = c(0,0.6))+ xlab("Temperatura(°C)")+ylab("velocidad(m/s)")+ 
  theme(legend.position = "top", 
        legend.title=element_blank())

desloc$tratamiento <- as.factor(desloc$tratamiento)

desloc$tratamiento <-relevel(desloc$tratamiento, ref="AE")

ggplot(data = desloc,aes(x = temp, y =vel, group=tratamiento))+
  facet_grid(.~replica)+geom_spline(aes(x=temp, y=vel,color=tratamiento))+theme_bw()+
  coord_cartesian(xlim = c(25,40),ylim = c(0,0.6))+ xlab("Temperatura(°C)")+ylab("velocidad(m/s)")+ 
  theme(legend.position = "top", 
        legend.title=element_blank())+scale_colour_brewer(palette = "Oranges")

```

ggplot(desloc, aes(x=temp, y=vel,color=tratamiento)) +
  geom_point() + geom_smooth(se =FALSE) +facet_grid(.~replica) +scale_color_brewer(palette="Oranges")+ 
  theme_bw()+labs(x = "Temperatura(°C)", y = "Velocidad(m/s)")+
  theme(plot.title = element_text(size = 14, family = "Tahoma", face = "bold"), text = element_text(size = 14, family = "Tahoma"),axis.title = element_text(face="bold"),axis.text.x=element_text(size = 12))

ggplot(desloc, aes(x=temp, y=vel,color=tratamiento)) +
  geom_point() + geom_smooth(method='lm', formula= y~poly(x,3), se=FALSE) +facet_grid(.~replica) +scale_color_brewer(palette="Oranges")+ 
  theme_bw()+labs(x = "Temperatura(°C)", y = "Velocidad(m/s)")+
  theme(plot.title = element_text(size = 14, family = "Tahoma", face = "bold"), text = element_text(size = 14, family = "Tahoma"),
        axis.title = element_text(face="bold"),axis.text.x=element_text(size = 12))

ggplot(desloc, aes(x=temp, y=vel,color=tratamiento)) +
  geom_point() + geom_smooth(aes(fill=tratamiento)) + scale_color_brewer(palette="Oranges")+ 
  theme_bw() +scale_fill_brewer(palette="Oranges")+labs(x = "Temperatura(°C)", y = "Velocidad(m/s)")+
  theme(plot.title = element_text(size = 14, family = "Tahoma", face = "bold"), text = element_text(size = 14, family = "Tahoma"),
        axis.title = element_text(face="bold"),axis.text.x=element_text(size = 12))

Aceleración
ggplot(desloc, aes(x=temp, y=acel,color=tratamiento)) +
  geom_point() + geom_smooth(aes(fill=tratamiento)) + scale_color_brewer(palette="Oranges")+ 
  theme_bw() +scale_fill_brewer(palette="Oranges")+labs(x = "Temperatura(°C)", y = "Aceleración(m/s^2)")+
  theme(plot.title = element_text(size = 14, family = "Tahoma", face = "bold"), text = element_text(size = 14, family = "Tahoma"),
        axis.title = element_text(face="bold"),axis.text.x=element_text(size = 12))

Máxima distancia
ggplot(desloc, aes(x=temp, y=long,color=tratamiento)) +
  geom_point() + geom_smooth(aes(fill=tratamiento)) + scale_color_brewer(palette="Oranges")+ 
  theme_bw() +scale_fill_brewer(palette="Oranges")+labs(x = "Temperatura(°C)", y = "Máxima Distancia(m)")+
  theme(plot.title = element_text(size = 14, family = "Tahoma", face = "bold"), text = element_text(size = 14, family = "Tahoma"),
        axis.title = element_text(face="bold"),axis.text.x=element_text(size = 12))

Modelo MixtoGAM

effect size

Paquetes para hacer el análisis

library(lme4)
library(mgcv)
library(gamm4)

Los modelos que probé:

mod.gam = gam(vel ~ s(temperatura)+s(tratamiento), data=desloc)

Error in smooth.construct.tp.smooth.spec(object, dk\(data, dk\)knots) : A term has fewer unique covariate combinations than specified maximum degrees of freedom

mixtogam<-gamm(vel~s(tratamiento)+s(temp)+s(replica),data=desloc)

Error in names(dat) <- object$term : 'names' attribute [1] must be the same length as the vector [0]

mixtogam<-gamm(vels(tratamiento)+s(temp),random=list(replica=1), data=desloc)

Error in smooth.construct.tp.smooth.spec(object, dk\(data, dk\)knots) : A term has fewer unique covariate combinations than specified maximum degrees of freedom