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_rep1.txt",h=T)
head(datos_morfo)
##   Id_Random_Foto Tratamiento Postura Caja Individuo Medicion Altura_Cabeza
## 1             23          AB      P1    1         1        1         2.167
## 2             23          AB      P1    1         1        2         2.129
## 3             23          AB      P1    1         1        3         2.145
## 4            112          AB      P1    1         2        1         2.437
## 5            112          AB      P1    1         2        2         2.536
## 6            112          AB      P1    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)
library(sjPlot) ###para hacer las tablas

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: 451.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.0576 -0.1748 -0.0067  0.1780  5.0469 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 0.89042  0.9436  
##  Residual                   0.03106  0.1762  
## Number of obs: 525, groups:  Id_Random_Foto, 175
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.11690    0.07174 173.99999   85.26   <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.966
tab_model(mixtolong_ca)
  Longitud_Cabeza
Predictors Estimates CI p
(Intercept) 6.12 5.98 – 6.26 <0.001
Random Effects
σ2 0.03
τ00 Id_Random_Foto 0.89
ICC 0.97
N Id_Random_Foto 175
Observations 525
Marginal R2 / Conditional R2 0.000 / 0.966

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: -437.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4662 -0.3335 -0.0371  0.3043  4.2060 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 0.087903 0.29648 
##  Residual                   0.007694 0.08772 
## Number of obs: 525, groups:  Id_Random_Foto, 175
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.01390    0.02274 173.99998   88.58   <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.920
tab_model(mixtoalt_cabeza)
  Altura_Cabeza
Predictors Estimates CI p
(Intercept) 2.01 1.97 – 2.06 <0.001
Random Effects
σ2 0.01
τ00 Id_Random_Foto 0.09
ICC 0.92
N Id_Random_Foto 175
Observations 525
Marginal R2 / Conditional R2 0.000 / 0.920

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: 1621.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.0420 -0.1362 -0.0116  0.1458  5.3295 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 8.9103   2.9850  
##  Residual                   0.2803   0.5295  
## Number of obs: 525, groups:  Id_Random_Foto, 175
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  20.4855     0.2268 174.0000   90.31   <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.969
tab_model(mixtolong_total)
  Longitud_Total
Predictors Estimates CI p
(Intercept) 20.49 20.04 – 20.93 <0.001
Random Effects
σ2 0.28
τ00 Id_Random_Foto 8.91
ICC 0.97
N Id_Random_Foto 175
Observations 525
Marginal R2 / Conditional R2 0.000 / 0.969

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: -910.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0375 -0.3618 -0.0058  0.3258  4.0340 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 0.058557 0.24199 
##  Residual                   0.002456 0.04955 
## Number of obs: 525, groups:  Id_Random_Foto, 175
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.30681    0.01842 173.99999   70.95   <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.960
tab_model(mixtoalt_cola)
  Altura_Cola
Predictors Estimates CI p
(Intercept) 1.31 1.27 – 1.34 <0.001
Random Effects
σ2 0.00
τ00 Id_Random_Foto 0.06
ICC 0.96
N Id_Random_Foto 175
Observations 525
Marginal R2 / Conditional R2 0.000 / 0.960

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: 1286.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3658 -0.2034 -0.0094  0.1673  5.3553 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 4.6936   2.1665  
##  Residual                   0.1478   0.3845  
## Number of obs: 525, groups:  Id_Random_Foto, 175
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  14.3673     0.1646 174.0000   87.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(mixtolong_cola)
## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.969
##   Conditional ICC: 0.969
tab_model(mixtolong_cola)
  Longitud_Cola
Predictors Estimates CI p
(Intercept) 14.37 14.04 – 14.69 <0.001
Random Effects
σ2 0.15
τ00 Id_Random_Foto 4.69
ICC 0.97
N Id_Random_Foto 175
Observations 525
Marginal R2 / Conditional R2 0.000 / 0.969

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: 2172.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0711 -0.3258  0.0001  0.3389  4.1643 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Id_Random_Foto (Intercept) 23.4744  4.8450  
##  Residual                    0.8357  0.9142  
## Number of obs: 525, groups:  Id_Random_Foto, 175
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  30.7450     0.3684 174.0000   83.45   <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.966
tab_model(mixtoperim_mus)
  Perim_Musculo1
Predictors Estimates CI p
(Intercept) 30.75 30.02 – 31.47 <0.001
Random Effects
σ2 0.84
τ00 Id_Random_Foto 23.47
ICC 0.97
N Id_Random_Foto 175
Observations 525
Marginal R2 / Conditional R2 0.000 / 0.966

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 esta no es afectada por los valores extremos y se recomienda mas cuando los valores no tienen distribuución normal.

datos_morfo_med <- read.table("morfo_med_completa_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=grupo_datos, y=Longitud_Total, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "PoGrupo datos",y = "Longitud Total (mm)")

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)")

Peso

ggplot(datos_morfo_med, aes(x=Postura, y=Peso, fill=Tratamiento)) +
  geom_boxplot()+scale_fill_brewer(palette = "Oranges") + theme_bw()+labs(x = "postura",y = "Peso (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.960 0.0461 
## 2 AB          Altura_Cabeza     0.989 0.904  
## 3 AD          Altura_Cabeza     0.937 0.00349
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.967 0.108     
## 2 AB          Longitud_Total     0.986 0.772     
## 3 AD          Longitud_Total     0.858 0.00000466
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.977 0.325  
## 2 AB          Altura_Cola     0.963 0.0977 
## 3 AD          Altura_Cola     0.943 0.00699
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.980 0.442    
## 2 AB          Longitud_Cola     0.984 0.698    
## 3 AD          Longitud_Cola     0.880 0.0000223
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.970 0.139    
## 2 AB          Perimetro_Musculo     0.988 0.852    
## 3 AD          Perimetro_Musculo     0.888 0.0000448
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.000450
## 2 AB          Longitud_Cabeza     0.943 0.0120  
## 3 AD          Longitud_Cabeza     0.913 0.000358

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.5161 0.9567
##       154
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  3.0294 0.05093 .
##       172                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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.6739  0.511
##       172
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.0103 0.3662
##       172
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.0209 0.01965 *
##       172                  
## ---
## 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.8058 0.06323 .
##       172                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Peso
# Levene's test with multiple independent variables
leveneTest(Peso~Tratamiento, data = datos_morfo_med)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group   2  5.4672 0.004989 **
##       172                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis de componentes pricipales

Correlación entre variables

morfo<-datos_morfo_med[,8:13] #primero los individuos activos y luego las variables
head(morfo)
##    Peso Altura_Cabeza Altura_Cola Longitud_Cola Perimetro_Musculo
## 1 0.036         2.145       1.691        17.326            37.880
## 2 0.027         2.536       1.336        14.029            30.549
## 3 0.043         2.328       1.495        17.072            37.172
## 4 0.025         2.172       1.760        17.075            37.207
## 5 0.036         2.062       1.496        15.598            34.448
## 6 0.042         2.283       1.455        15.270            33.554
##   Longitud_Cabeza
## 1           6.600
## 2           7.144
## 3           8.039
## 4           6.721
## 5           6.555
## 6           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=175)
## $chisq
## [1] 1396.628
## 
## $p.value
## [1] 8.798109e-289
## 
## $df
## [1] 15

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.86
## MSA for each item = 
##              Peso     Altura_Cabeza       Altura_Cola     Longitud_Cola 
##              0.94              0.91              0.94              0.78 
## Perimetro_Musculo   Longitud_Cabeza 
##              0.78              0.91
Componentes:

componentes factorMineR

#PCA usando factorMineR y factorextra
elpca<-PCA(morfo, scale.unit = TRUE,ncp=6,graph = FALSE)
summary(elpca)
## 
## Call:
## PCA(X = morfo, scale.unit = TRUE, ncp = 6, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6
## Variance               4.909   0.386   0.272   0.246   0.174   0.013
## % of var.             81.811   6.439   4.534   4.103   2.903   0.210
## Cumulative % of var.  81.811  88.250  92.785  96.887  99.790 100.000
## 
## Individuals (the 10 first)
##                       Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## 1                 |  2.806 |  2.609  0.792  0.865 | -0.569  0.479  0.041 |
## 2                 |  2.077 |  1.182  0.163  0.324 |  1.465  3.177  0.498 |
## 3                 |  3.382 |  3.247  1.228  0.922 | -0.131  0.025  0.001 |
## 4                 |  2.727 |  2.339  0.637  0.736 |  0.066  0.007  0.001 |
## 5                 |  1.644 |  1.527  0.271  0.863 | -0.464  0.318  0.079 |
## 6                 |  2.595 |  2.295  0.613  0.782 |  0.188  0.052  0.005 |
## 7                 |  0.911 | -0.720  0.060  0.624 |  0.478  0.338  0.275 |
## 8                 |  1.782 | -1.466  0.250  0.676 |  0.754  0.842  0.179 |
## 9                 |  1.924 | -1.545  0.278  0.645 | -0.297  0.130  0.024 |
## 10                |  1.181 |  0.144  0.002  0.015 |  0.568  0.478  0.232 |
##                    Dim.3    ctr   cos2  
## 1                  0.800  1.343  0.081 |
## 2                 -0.722  1.096  0.121 |
## 3                 -0.866  1.576  0.066 |
## 4                  1.022  2.194  0.140 |
## 5                  0.206  0.089  0.016 |
## 6                 -0.809  1.374  0.097 |
## 7                  0.011  0.000  0.000 |
## 8                  0.429  0.387  0.058 |
## 9                  0.975  1.997  0.257 |
## 10                -0.283  0.169  0.058 |
## 
## Variables
##                      Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Peso              |  0.884 15.928  0.782 | -0.282 20.527  0.079 | -0.058  1.238
## Altura_Cabeza     |  0.862 15.135  0.743 |  0.440 50.149  0.194 | -0.029  0.317
## Altura_Cola       |  0.884 15.922  0.782 |  0.125  4.052  0.016 |  0.396 57.702
## Longitud_Cola     |  0.945 18.209  0.894 | -0.213 11.696  0.045 |  0.014  0.076
## Perimetro_Musculo |  0.954 18.533  0.910 | -0.178  8.158  0.032 |  0.010  0.040
## Longitud_Cabeza   |  0.894 16.273  0.799 |  0.145  5.418  0.021 | -0.332 40.629
##                     cos2  
## Peso               0.003 |
## Altura_Cabeza      0.001 |
## Altura_Cola        0.157 |
## Longitud_Cola      0.000 |
## Perimetro_Musculo  0.000 |
## Longitud_Cabeza    0.111 |
#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    2.608957825 -0.56907394  0.799633373 -0.173122420 -0.260700468
## 2    1.182030543  1.46547465 -0.722476721  0.209701931  0.449748936
## 3    3.247384178 -0.13061518 -0.866141097  0.306373442 -0.182925545
## 4    2.339435404  0.06643850  1.022062907 -0.507187287 -0.811798924
## 5    1.527019762 -0.46351706  0.205692753  0.309522794 -0.028762565
## 6    2.294529136  0.18807370 -0.808750864  0.879268541  0.009825349
## 7   -0.719888086  0.47770433  0.010794894  0.255833178 -0.133785523
## 8   -1.465861170  0.75448868  0.429232614 -0.338803144  0.389143732
## 9   -1.545003466 -0.29661347  0.975056319 -0.095313121 -0.514270134
## 10   0.144322268  0.56819434 -0.283280337 -0.334630877 -0.917273332
## 11  -1.439779454  0.25438553 -0.019625750  0.062105316 -0.027237993
## 12  -2.432430290 -0.10885656  0.561197950  0.284597622 -0.027559422
## 13   1.556037887  0.23115594  0.562267727 -0.708834391 -0.540334987
## 14   1.238716103  0.61898811 -0.603652792 -0.588162186 -1.008886565
## 15  -2.036430098 -0.25544250 -0.669733545 -0.243225424 -0.606092799
## 16  -2.448223566  0.21423313  0.377101951  0.100737197 -0.032691959
## 17  -1.619886712 -0.06866157 -0.512092353 -0.581298277  0.350207912
## 18  -1.121778566  0.31026871 -0.542465176  0.024261485 -0.148676483
## 19  -1.011072444  0.55843478 -0.357778109  0.081796451  0.299863838
## 20  -1.538342103 -0.25345553  0.098724181  0.005490799  0.181575681
## 21  -0.576787106  0.04701803  0.824139543  0.201775063  0.426737347
## 22  -2.984510113 -0.45225356  0.277199705 -0.112049929  0.026101488
## 23  -0.697016985 -0.07788810  0.254822234  0.305640093 -0.013580848
## 24  -0.968143215 -0.03483548  0.973681927  0.521409700 -0.358224411
## 25   0.760443740 -0.39400141  0.015139024 -0.166551934  0.042381902
## 26  -0.377582042 -0.07742470  0.453489254 -0.342264676  0.247949729
## 27   0.779035557  0.12557039  0.320460142 -0.524285066  0.861459637
## 28  -0.454531291  0.36333573 -0.296942235 -0.157545186  0.649877911
## 29   0.625272897 -0.07816309  0.466733871 -0.332756621  0.093248853
## 30  -0.136677577  0.06340631  0.404408957 -0.250327987 -0.437506959
## 31  -0.124232471 -0.55472180 -0.440732411 -0.540497029  0.646374432
## 32  -2.581664348 -0.95751979 -0.396239649 -0.287739360 -0.192354071
## 33  -0.381422553 -0.24463761 -0.353227938 -0.411352540  0.007230492
## 34  -0.996880660 -0.48600243 -0.110725444 -0.843966825  0.212962820
## 35   8.511332636 -1.25884034 -0.055975683  1.122988880  0.643884344
## 36   7.628778808 -0.73191341 -0.454074063  0.700036899  0.621177871
## 37   6.662731811  1.04115972  0.003708405 -0.560009273  0.235884456
## 38   6.849209358 -2.58697646 -0.196170982  1.055622978  0.324497749
## 39   6.849723876 -1.07907875  0.138969754  1.056293577 -0.205650869
## 40   1.101949077  0.66685254  0.379063675  0.410109440 -0.054062174
## 41   0.670849660  0.31410488 -0.997728307  0.412099190 -0.513568404
## 42   0.761393696 -0.00943285  0.123774933  0.277373126 -0.293149502
## 43   4.079614040  0.48765422  0.397766613 -0.705043891 -1.404945625
## 44   0.404975593 -0.19399853 -0.040101501 -0.018644383 -0.020432048
## 45   4.328942007  0.12910576  0.514005819 -0.069182020 -0.532014617
## 46   0.390980091  0.59040507  0.041416709 -0.171373204  0.262538563
## 47  -0.088024240  0.03898466 -0.342037336 -0.090491655 -0.044946135
## 48   1.298390209  0.46268630 -0.407138251 -0.621878075 -0.944206472
## 49  -1.359142418  0.32528680 -0.201126964  0.175512956  0.286751985
## 50  -1.786656058  0.08923760 -0.065195353 -0.266014872 -0.160116215
## 51  -0.412723017  0.68342560 -0.380532800  0.292002264 -0.001124135
## 52  -1.523597380  0.12566817 -0.680225743 -0.709643592 -0.084993418
## 53   1.085829945  1.35159137 -0.018631692 -0.606364989 -0.227802403
## 54  -0.281825049  0.37586706 -0.984724539 -1.087750490  0.199828895
## 55  -1.288151599 -0.27685946  0.068736920 -0.016566083 -0.182798968
## 56  -2.170587826 -0.46229362  0.057736478  0.067049199  0.251968255
## 57   1.099431829  0.48885133  0.342931481 -0.405065441 -0.861607454
## 58   0.894521695  0.51751233  0.974612247  1.187199778 -0.193255373
## 59   1.068068687  0.55846870  0.086701085  0.733228851  0.568881104
## 60  -0.125014683 -0.18676304  0.073491789 -0.081513881  0.200687617
## 61  -0.410684858 -0.21368657  0.233667430  0.285867836  0.152722919
## 62   4.868094649  1.32986572  0.757764470 -0.817187779  0.291877418
## 63  -1.309440204 -0.22395963  0.080887553 -0.050004357  0.254979031
## 64   3.273854648  1.42887784  1.138050708 -0.275679086  0.908254406
## 65   0.369523374  0.06477640  0.820630228 -0.681397600  0.287938246
## 66   0.379184566  0.57117621  0.488643770 -0.185194834  0.699307439
## 67   2.648754178 -0.62302222  0.585383565 -0.236060067 -0.304495785
## 68  -0.301923021 -0.20307180 -0.179585685 -0.102062752  0.122131220
## 69   0.554223090 -0.46727149 -0.059131727 -0.616426073 -0.135785047
## 70   1.865995424 -0.17973231  0.492991662  0.289714718 -0.122990549
## 71   0.957598083 -0.36050543 -0.050119025  0.220892066 -0.246132811
## 72   1.573534187 -1.64848644  0.344588054 -1.818267570  1.201550393
## 73  -2.047250969 -0.64063012 -0.082372980 -0.362992406  0.155680410
## 74  -0.039400241 -0.24401989 -0.026267235  0.236448063 -0.344024095
## 75  -0.344856584 -0.26245051  0.112353589 -0.509695986 -0.136685816
## 76   1.082566861  0.01928532  0.380794697  0.313513814 -0.176842274
## 77   2.781021939  0.83867063 -0.865321264 -0.938599438 -1.179770948
## 78  -0.146609576 -0.76586463 -0.570682754 -0.152783781 -0.374453004
## 79  -0.719037371 -0.18194237 -0.241322099  0.036705315  0.130186213
## 80   1.255942705 -0.09679998 -0.193539407 -0.588715430 -0.263739623
## 81   1.474282941 -0.13008913  0.060666169 -0.374215558  0.153087689
## 82   0.418796772 -0.06700617 -0.589119159  0.366439780 -0.004002893
## 83   3.185439819  0.21945023 -0.100365965 -0.639343096 -0.613971242
## 84   0.340890801 -1.07933372  0.054082757 -0.149579253  0.063814838
## 85  -1.451849118  0.27195921 -0.066607339  0.098568331  0.109373185
## 86   3.156710981  0.35714033 -0.093654053  0.016192496  0.710221435
## 87  -0.918569396 -0.58762524 -0.174907128 -0.086866995  0.090712590
## 88  -0.769594347  0.73710438 -0.237656490  0.384608843  0.410852362
## 89  -1.513911135 -0.36595397 -0.007408362  0.068337792  0.121816272
## 90  -2.409793127 -0.26620335  0.209315636 -0.040814165 -0.042290314
## 91   2.068333855  1.50168932 -1.076731949 -0.500048891 -0.594573200
## 92  -2.916098329  0.04306536 -0.615029616  0.347104678  0.317367078
## 93  -1.254427744  0.46369503  0.639260818  0.053319447 -0.145619730
## 94   1.672143047  0.24705319  0.087864960 -1.389173732  0.777468116
## 95  -1.521859895 -0.11136528  0.519157157 -0.300344110 -0.054191172
## 96  -1.759512143 -0.56776293 -0.126496598 -0.691550004  0.073245327
## 97   1.973670808  0.30195881 -0.005276784 -0.720349524  0.058038677
## 98   4.715672219  0.71479836 -0.629335268 -0.329946620 -0.580006971
## 99  -0.039332941 -0.08252692 -0.447317173  0.048861729  0.112385930
## 100 -0.897120236 -0.60728706  0.003277711 -0.283472668  0.171373190
## 101  1.846028907  0.17185779 -0.804979913  1.202248872 -0.379686925
## 102 -2.966970184 -0.50075695 -0.156961960  0.483263593 -0.040788925
## 103 -0.158358615 -0.60143464  0.142857076 -0.571888954 -0.025913369
## 104 -1.167466509  0.52938040 -0.341726200 -0.512884899  0.479028177
## 105 -1.443669286  0.25482104  0.116279122 -0.256120169  0.566018453
## 106 -0.668200383  1.12914353 -1.173744541  0.771449392 -0.332236525
## 107 -0.313426346  0.58663993  0.716194519 -0.177745309  0.654065691
## 108 -2.072148321 -0.06330471 -0.378537587 -0.218867842  0.295934540
## 109 -3.239591314 -0.05701732 -0.373868727 -0.325573045  0.885169015
## 110  0.650244063  0.15538648  0.471261760 -0.372320571  0.670597423
## 111  1.589576080  1.18106702  2.157109671  0.086256404 -0.325307170
## 112 -1.947736130 -0.62875728 -0.588208565 -0.431881231 -0.259837361
## 113  0.282629591 -0.44851828 -0.506646982 -0.601496025  0.490659618
## 114  0.376503356  0.16432282 -0.848266677  0.351159596 -0.282800953
## 115  1.657778519  0.30122114 -0.039404916 -0.119036978 -0.258313548
## 116 -3.916057698  0.77777907 -0.134624858  0.786024082  0.540468476
## 117 -1.621753367  0.06732067  0.372825603  0.213837657  0.191520609
## 118  1.168392265  0.97792154 -0.593321993 -0.288152863 -0.641367746
## 119 -3.157987959  0.60450033 -0.132615390  0.158133845  0.548265193
## 120 -2.532036301 -0.02797345  0.379737235  0.060812231  0.035573395
## 121 -1.997280348 -0.59962676  1.364958497  0.292201499 -0.492450812
## 122 -1.697566829  0.50296562  0.441103732  0.594191969  0.021545340
## 123 -3.390674756 -0.31868616  0.199795462  0.551607453 -0.190899958
## 124 -0.829588794  0.64282987 -0.128302222 -0.059691304  0.642228991
## 125 -2.012414820 -0.01067511 -0.099484720  0.078828527  0.494803602
## 126 -3.526142950  0.01989272  0.180001335  0.438901210  0.361510873
## 127 -1.595772927  0.23981994 -0.485243424  0.398631755 -0.094738120
## 128 -2.057785553 -0.55409084  0.312151623  0.116399246 -0.251319611
## 129 -0.741870585 -0.30953764 -0.139464331 -0.360055097  0.065691806
## 130 -1.397839322 -0.86241187 -0.728475598  0.038821525 -0.207871477
## 131 -1.645064413 -0.66459143 -0.455952686 -0.053785294 -0.147691015
## 132 -3.081659765 -0.80536922 -0.224938200  0.506978225 -0.262033315
## 133  6.647103810 -1.66809317 -0.954434841 -0.063066369  0.590568433
## 134  1.706695531 -0.71049289  0.331509578 -0.153787239 -0.108952131
## 135  5.711084144 -0.53753299 -0.090197062  0.866734356  0.755291399
## 136  0.059145062  0.18088764  0.072432302  0.628382823  0.343726584
## 137  0.769397914 -0.37369175  0.834055377  0.733278271 -0.531449217
## 138 -1.613509522 -0.05641079  0.173068335  0.537126903 -0.210671610
## 139 -0.441582003  0.78290356 -0.450522940 -0.145592214 -0.334429135
## 140 -0.184073356  0.16252655  0.026118575  0.308262583 -0.081974552
## 141 -0.832756411  0.25988047 -0.547896595  0.908439351 -0.345416087
## 142 -1.917749070 -0.38899433 -0.571496640 -0.140599452  0.240614565
## 143 -2.684807542 -0.32588193 -0.023674544  0.229469976  0.079859355
## 144 -2.290215742 -0.14609278 -0.511594908  0.313982959 -0.285618203
## 145  1.366947491  0.55760995 -0.364292618  0.709449286  0.104422332
## 146 -1.397457466  0.55648357 -0.660638265  1.034031416  0.036661155
## 147  0.077616571  1.00787361 -0.917314463  0.635001836  0.030462297
## 148 -0.006675255 -0.84264008  0.352061001 -0.373348175 -0.391514624
## 149 -0.251207181 -0.61690446 -0.192179750 -0.065924769  0.005600501
## 150 -1.740127594 -0.47591162 -0.264813369 -0.511047913  0.216057590
## 151  0.662555811 -0.12516683  1.021091522 -0.123041876 -0.103194458
## 152  0.829174933 -0.52526825  0.468958093 -0.138458444 -0.205852488
## 153 -0.396487321 -0.96859486  0.185852576 -0.083860368 -0.585677762
## 154 -2.257631618 -0.62535233  0.507353248 -0.054692907 -0.559080272
## 155 -0.140374138  0.40110735  0.522580293  0.483320434 -0.336940104
## 156  0.968371625  0.81961014  0.563467047  0.005683113  0.417549090
## 157 -2.865054801  0.01857965  0.310802733  0.549942189 -0.114460662
## 158  0.908610658  0.05641274  0.173943867  0.738845378 -0.235425591
## 159  1.717033935  0.97264194  0.572861939  0.146741955  0.610680844
## 160 -0.491871429  0.10159338  0.879042246  0.421324167 -0.176404997
## 161 -1.655524658 -0.18236271 -0.203847081  0.133409761  0.286963647
## 162 -1.184071278  0.08762429  0.387603658  0.184821016 -0.101540075
## 163  0.783118354  0.51205771  0.158890539  0.338805726  0.278421467
## 164 -2.621112437  0.35075678  0.671421968  0.785698105 -0.383888096
## 165 -0.391305107  1.68248757 -0.505074386 -0.568698174  0.628331054
## 166 -2.602848679 -0.13535076  0.378025419  0.313774087  0.009659203
## 167 -0.676937642  0.44338979 -0.134639391  0.513501315  0.297946138
## 168 -0.129835524 -0.03254117  0.006579340  0.548780709 -0.113529046
## 169  1.690342186  0.48594056 -0.998481686  0.686861789  0.320373729
## 170 -1.260641689 -0.33347487  0.141642107  0.257141979 -0.077482078
## 171 -2.074595362 -0.81058528  0.356282575 -0.111021082 -0.509337963
## 172 -1.872653651 -0.10105046  0.157379245  0.239488110 -0.009126399
## 173 -1.257575827 -0.56336315 -0.014544293 -0.401175912 -0.001460952
## 174 -0.647607803 -1.12099418 -0.452632834 -0.842281648 -0.221276643
## 175 -1.544002667 -1.02589574 -0.912614534 -0.270328302 -0.556573620
##            Dim.6
## 1   -0.059995846
## 2   -0.026620734
## 3   -0.010827964
## 4    0.001319216
## 5   -0.131421247
## 6   -0.071249843
## 7    0.008421899
## 8   -0.093242000
## 9    0.047846708
## 10  -0.127637863
## 11  -0.097957622
## 12   0.020846477
## 13   0.050622066
## 14  -0.041017217
## 15  -0.085435141
## 16  -0.045010347
## 17  -0.017588165
## 18  -0.260463583
## 19  -0.116066418
## 20  -0.078658337
## 21  -0.089921907
## 22  -0.121926485
## 23  -0.098417275
## 24  -0.104023810
## 25  -0.150700970
## 26  -0.089681518
## 27   0.011914329
## 28  -0.037820003
## 29  -0.056830805
## 30  -0.070900817
## 31  -0.052342618
## 32  -0.027116174
## 33  -0.016040254
## 34  -0.011753732
## 35  -0.140023087
## 36  -0.121947411
## 37   0.029582031
## 38  -0.137539049
## 39  -0.143726666
## 40  -0.046665067
## 41  -0.147778535
## 42  -0.074370583
## 43  -0.164954562
## 44  -0.238568423
## 45   0.066188161
## 46  -0.134334528
## 47  -0.168330916
## 48   0.026323230
## 49  -0.105097432
## 50  -0.052725531
## 51  -0.078726946
## 52   0.050487055
## 53  -0.069258063
## 54  -0.025410657
## 55  -0.024426966
## 56  -0.039119953
## 57   0.006046927
## 58  -0.054969138
## 59   0.012480262
## 60  -0.027470289
## 61  -0.013910858
## 62   0.107404114
## 63  -0.067251536
## 64  -0.015036529
## 65   0.042687312
## 66  -0.188701579
## 67   0.055611217
## 68  -0.169655655
## 69  -0.034853132
## 70  -0.197934368
## 71  -0.080422766
## 72   0.033783984
## 73  -0.114322681
## 74  -0.117780049
## 75   0.072053075
## 76   0.005753144
## 77  -0.045539858
## 78  -0.114861544
## 79  -0.113409715
## 80  -0.063658791
## 81   0.019704999
## 82   0.008829135
## 83   0.023773980
## 84  -0.268130446
## 85  -0.079830914
## 86   0.060553880
## 87  -0.072507380
## 88  -0.192323751
## 89  -0.086695891
## 90  -0.155494925
## 91  -0.119019162
## 92  -0.093668429
## 93   0.068789621
## 94   0.012014650
## 95  -0.076034651
## 96   0.120937030
## 97  -0.077872588
## 98  -0.056483462
## 99  -0.043547969
## 100 -0.113822046
## 101 -0.154165317
## 102 -0.101680368
## 103 -0.056162110
## 104 -0.125989965
## 105 -0.117118904
## 106 -0.174595608
## 107 -0.165691390
## 108 -0.078889019
## 109 -0.013409476
## 110 -0.016744195
## 111  0.059678826
## 112  0.077917224
## 113  0.161270182
## 114  0.168716777
## 115  0.194297198
## 116  0.017483566
## 117  0.089887921
## 118  0.278145133
## 119 -0.011210522
## 120  0.058009319
## 121  0.004241535
## 122 -0.027345982
## 123  0.011065588
## 124  0.164739170
## 125  0.062599757
## 126 -0.024653694
## 127  0.076125904
## 128  0.059576486
## 129  0.105292781
## 130  0.111460664
## 131  0.104495417
## 132  0.046179761
## 133  0.280197739
## 134  0.140831768
## 135  0.218146176
## 136  0.113435397
## 137  0.101678057
## 138  0.069421637
## 139  0.177456105
## 140  0.117840042
## 141  0.091978101
## 142  0.097844878
## 143  0.024690172
## 144  0.066620415
## 145  0.231763208
## 146  0.115525100
## 147  0.143912508
## 148  0.118028192
## 149  0.132394224
## 150  0.101912947
## 151  0.127569248
## 152  0.231621475
## 153  0.169871592
## 154  0.083105545
## 155  0.092160143
## 156  0.120717282
## 157 -0.001402911
## 158  0.167114929
## 159  0.164617620
## 160  0.005124170
## 161  0.028587186
## 162  0.120754829
## 163  0.141764277
## 164  0.000621668
## 165  0.137273069
## 166  0.041949940
## 167  0.357283845
## 168  0.082486988
## 169  0.163797998
## 170  0.073410107
## 171  0.087250443
## 172  0.042624618
## 173  0.119282023
## 174  0.164279578
## 175  0.145839691
# Export into a TXT file
write.infile(elpca$ind$coord, "pca_completo.txt", sep = "\t")

# Export into a CSV file
write.infile(elpca$ind$coord, "pca_completo.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<-morforot$scores[1:175,]

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

Figuras con el PC para ver el efecto del tratamiento

datos_morfo_med <- read.table("morfo_med_completa_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
##    Peso Altura_Cabeza Altura_Cola Longitud_Cola Perimetro_Musculo
## 1 0.036         2.145       1.691        17.326            37.880
## 2 0.027         2.536       1.336        14.029            30.549
## 3 0.043         2.328       1.495        17.072            37.172
## 4 0.025         2.172       1.760        17.075            37.207
## 5 0.036         2.062       1.496        15.598            34.448
## 6 0.042         2.283       1.455        15.270            33.554
##   Longitud_Cabeza grupo_datos  Dim.1   Dim.2    PC1rot
## 1           6.600      desloc 2.6090 -0.5691 1.1741947
## 2           7.144      desloc 1.1820  1.4655 0.5319879
## 3           8.039      desloc 3.2474 -0.1306 1.4615267
## 4           6.721      desloc 2.3394  0.0664 1.0528927
## 5           6.555      desloc 1.5270 -0.4635 0.6872547
## 6           7.736      desloc 2.2945  0.1881 1.0326821
#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=Dim.1, 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=PC1rot, 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:

Modelo mixto con PC1

mixtopc1<-lmer(Dim.1~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: Dim.1 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 686.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5505 -0.6150 -0.1473  0.5354  2.7810 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 1.173    1.083   
##  Postura  (Intercept) 1.251    1.119   
##  Residual             2.381    1.543   
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)    -0.3371     0.6213 12.5049  -0.543   0.5970  
## TratamientoAB  -0.4517     0.6481 12.2773  -0.697   0.4988  
## TratamientoAD   1.3648     0.6434 11.9367   2.121   0.0555 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.515       
## TratamintAD -0.519  0.497

Aqui otra forma de presentar los resultados, calculando el R2 del modelo y haciendo un ajuste en la significancia:

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)
  Dim.1
Predictors Estimates CI p
(Intercept) -0.34 -1.55 – 0.88 0.587
Tratamiento [AB] -0.45 -1.72 – 0.82 0.486
Tratamiento [AD] 1.36 0.10 – 2.63 0.034
Random Effects
σ2 2.38
τ00 Caja 1.17
τ00 Postura 1.25
ICC 0.50
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.111 / 0.560

Haciendo el ajuste en la significancia:

tab_model(mixtopc1, p.val = "kr", show.df = TRUE)
  Dim.1
Predictors Estimates CI p df
(Intercept) -0.34 -1.68 – 1.01 0.597 12.53
Tratamiento [AB] -0.45 -1.86 – 0.96 0.499 12.12
Tratamiento [AD] 1.36 -0.04 – 2.77 0.056 11.78
Random Effects
σ2 2.38
τ00 Caja 1.17
τ00 Postura 1.25
ICC 0.50
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.111 / 0.560
1.Prueba de homogeneidad del modelo mixto pc1

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

par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(mixtopc1)

Prueba de homogeneidad sin considerar los factores aleatorios (lineal model)
trat<-lm(Dim.1~Tratamiento,data=datos_morfo_med)
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(trat)

Prueba de homogeneidad Breusch Pagan

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("morfo_med_completa_expcasa.txt",h=T)


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

tra_pos<-lm(Dim.1~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 : Dim.1 
##  Variables: fitted values of Dim.1 
## 
##         Test Summary          
##  -----------------------------
##  DF            =    1 
##  Chi2          =    4.094103 
##  Prob > Chi2   =    0.04303307
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    4.97  4.9704  0.4512 0.5027
## Postura     6   56.98  9.4973  0.8622 0.5240
## Residuals 167 1839.52 11.0151
2.Evaluando los supuestos de los 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(Dim.1~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixtopc1_robust)
## Robust linear mixed model fit by DAStau 
## Formula: Dim.1 ~ Tratamiento + (1 | Caja) + (1 | Postura) 
##    Data: datos_morfo_med 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9617 -0.6533 -0.0823  0.5432  3.2641 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.3360   0.5796  
##  Postura  (Intercept) 0.7746   0.8801  
##  Residual             2.0389   1.4279  
## Number of obs: 175, groups: Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    -0.9195     0.4554  -2.019
## TratamientoAD   1.4072     0.4203   3.348
## TratamientoAE   0.3978     0.4211   0.945
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAD
## TratamintAD -0.475       
## TratamintAE -0.474  0.514
## 
## Robustness weights for the residuals: 
##  140 weights are ~= 1. The remaining 35 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.412   0.615   0.744   0.750   0.925   0.999 
## 
## Robustness weights for the random effects: 
##  26 weights are ~= 1. The remaining 2 ones are
##     8    22 
## 0.217 0.637 
## 
## 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)
  Dim.1
Predictors Estimates CI p
(Intercept) -0.92 -1.81 – -0.03 0.043
Tratamiento [AD] 1.41 0.58 – 2.23 0.001
Tratamiento [AE] 0.40 -0.43 – 1.22 0.345
Random Effects
σ2 2.04
τ00 Caja 0.34
τ00 Postura 0.77
ICC 0.35
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.101 / 0.418 
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, Dim.1 ~ 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, Dim.1 ~ fitted(.) | Caja, abline = c (0,1))

3. Calculando la magnitud del efecto del modelo

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.

calculo de eta-squared

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

eta_squared(mixtopc1)
## Parameter   | Eta2 (partial) |       90% CI
## -------------------------------------------
## Tratamiento |           0.41 | [0.02, 0.64]
library(emmeans)
trat.emm.s <- emmeans(mixtopc1, "Tratamiento")
pairs(trat.emm.s)
##  contrast estimate    SE   df t.ratio p.value
##  AE - AB     0.452 0.648 12.1  0.697  0.7698 
##  AE - AD    -1.365 0.643 11.8 -2.121  0.1278 
##  AB - AD    -1.816 0.648 12.1 -2.805  0.0391 
## 
## 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.293 0.424 12.5   -0.627   1.2129
##  AE - AD       -0.884 0.427 12.5   -1.811   0.0416
##  AB - AD       -1.177 0.421 12.5   -2.090  -0.2640
## 
## sigma used for effect sizes: 1.543 
## 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(Dim.1~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(Dim.1~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: Dim.1 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 686.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5505 -0.6150 -0.1473  0.5354  2.7810 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 1.173    1.083   
##  Postura  (Intercept) 1.251    1.119   
##  Residual             2.381    1.543   
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)    -0.7888     0.6256 12.8484  -1.261   0.2298  
## TratamientoAE   0.4517     0.6481 12.2773   0.697   0.4988  
## TratamientoAD   1.8165     0.6476 12.2377   2.805   0.0156 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAE
## TratamintAE -0.525       
## TratamintAD -0.525  0.507
plot_model(mixtopc_AB)

AB y AE
#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(Dim.1~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: Dim.1 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 686.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5505 -0.6150 -0.1473  0.5354  2.7810 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 1.173    1.083   
##  Postura  (Intercept) 1.251    1.119   
##  Residual             2.381    1.543   
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)     1.0277     0.6208 12.4615   1.655   0.1228  
## TratamientoAB  -1.8165     0.6476 12.2377  -2.805   0.0156 *
## TratamientoAE  -1.3648     0.6434 11.9367  -2.121   0.0555 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.514       
## TratamintAE -0.517  0.496
plot_model(mixtopc_AD)

Probando otras opciones para ver magnitud del efecto:

library(effectsize)
effectsize(mixtopc1)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |              -0.15 | [-0.70, 0.40]
## TratamientoAB |              -0.20 | [-0.77, 0.37]
## TratamientoAD |               0.61 | [ 0.05, 1.18]
## 
## # Standardization method: refit
effectsize(mixtopc_AB)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |              -0.36 | [-0.91, 0.20]
## TratamientoAE |               0.20 | [-0.37, 0.77]
## TratamientoAD |               0.82 | [ 0.25, 1.39]
## 
## # Standardization method: refit
effectsize(mixtopc_AD)
## Parameter     | Coefficient (std.) |         95% CI
## ---------------------------------------------------
## (Intercept)   |               0.46 | [-0.09,  1.01]
## TratamientoAB |              -0.82 | [-1.39, -0.25]
## TratamientoAE |              -0.61 | [-1.18, -0.05]
## 
## # Standardization method: refit

El modelo mixto con PC rotado

mixto1pc1<-lmer(PC1rot~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: PC1rot ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 411.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5505 -0.6150 -0.1473  0.5354  2.7810 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.2375   0.4874  
##  Postura  (Intercept) 0.2534   0.5034  
##  Residual             0.4823   0.6944  
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)     0.4625     0.2794 12.4615   1.655   0.1228  
## TratamientoAB  -0.8175     0.2914 12.2377  -2.805   0.0156 *
## TratamientoAE  -0.6142     0.2896 11.9367  -2.121   0.0555 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.514       
## TratamintAE -0.517  0.496
library(sjPlot)
library(lme4)
tab_model(mixto1pc1)
  PC1rot
Predictors Estimates CI p
(Intercept) 0.46 -0.09 – 1.01 0.098
Tratamiento [AB] -0.82 -1.39 – -0.25 0.005
Tratamiento [AE] -0.61 -1.18 – -0.05 0.034
Random Effects
σ2 0.48
τ00 Caja 0.24
τ00 Postura 0.25
ICC 0.50
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.111 / 0.560 
anova(mixto1pc1)
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Tratamiento 4.1353  2.0676     2 12.146  4.2874 0.03901 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1.Prueba de homogeneidad del modelo mixto rotado pc1

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

par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(mixto1pc1)

Prueba de homogeneidad sin considerar los factores aleatorios (lineal model)

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

Prueba de homogeneidad Breusch Pagan:

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

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


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

tra_pos_pcrot<-lm(PC1rot~Tra_Pos, data=datos_morfo_med)
plot(fitted(tra_pos_pcrot),residuals(tra_pos_pcrot))

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

Otra forma de evaluar la homogeneidad propuesta por Michael Palmeri

library(lmerTest)
library(lme4)
mixto1pc1<-lmer(PC1rot~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
#extracts the residuals and places them in a new column in our original data table
datos_morfo_med$mixto1pc1.Res<-residuals(mixto1pc1)
#creates a new column with the absolute value of the residuals
datos_morfo_med$Abs.mixto1pc1.Res<-abs(datos_morfo_med$mixto1pc1.Res)
#squares the absolute values of the residuals to provide the more robust estimate

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

anova(Levene.mixto1pc1) #displays the results
## Analysis of Variance Table
## 
## Response: mixto1pc1.Res2
##            Df Sum Sq Mean Sq F value Pr(>F)
## Caja        1  0.204 0.20394  0.4513 0.5027
## Postura     6  2.338 0.38967  0.8622 0.5240
## Residuals 167 75.474 0.45194

El modelo mixto con PC2

mixto1pc2<-lmer(Dim.2~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: Dim.2 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 318.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5608 -0.6316 -0.0924  0.5285  2.9662 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.03867  0.1967  
##  Postura  (Intercept) 0.05519  0.2349  
##  Residual             0.30900  0.5559  
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)    0.01762    0.13865 13.85620   0.127    0.901
## TratamientoAD -0.10882    0.14813 12.27026  -0.735    0.476
## TratamientoAE  0.04709    0.14842 12.36862   0.317    0.756
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAD
## TratamintAD -0.552       
## TratamintAE -0.551  0.516
library(sjPlot)
library(lme4)
tab_model(mixto1pc2)
  Dim.2
Predictors Estimates CI p
(Intercept) 0.02 -0.25 – 0.29 0.899
Tratamiento [AD] -0.11 -0.40 – 0.18 0.463
Tratamiento [AE] 0.05 -0.24 – 0.34 0.751
Random Effects
σ2 0.31
τ00 Caja 0.04
τ00 Postura 0.06
ICC 0.23
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.011 / 0.241 
anova(mixto1pc2)
## Type III Analysis of Variance Table with Satterthwaite's method
##              Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Tratamiento 0.37091 0.18546     2 12.066  0.6002 0.5643

Haciendo el ajuste en la significancia:

tab_model(mixto1pc2, p.val = "kr", show.df = TRUE)
  Dim.2
Predictors Estimates CI p df
(Intercept) 0.02 -0.28 – 0.32 0.901 13.93
Tratamiento [AD] -0.11 -0.43 – 0.21 0.477 12.17
Tratamiento [AE] 0.05 -0.28 – 0.37 0.757 12.27
Random Effects
σ2 0.31
τ00 Caja 0.04
τ00 Postura 0.06
ICC 0.23
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.011 / 0.241
1.Prueba de homogeneidad del modelo mixto pc2

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

par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(mixto1pc2)

Prueba de homogeneidad sin considerar los factores aleatorios (lineal model)

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

Prueba de homogeneidad Breusch Pagan:

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

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

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

tra_pos_pc2<-lm(Dim.2~Tra_Pos, data=datos_morfo_med)
plot(fitted(mixto1pc2),residuals(mixto1pc2))

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

Otra forma de evaluar la homogeneidad propuesta por Michael Palmeri

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

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

anova(Levene.mixto1pc2) #displays the results
## Analysis of Variance Table
## 
## Response: mixto1pc2.Res2
##            Df Sum Sq Mean Sq F value Pr(>F)
## Caja        1  0.404 0.40423  1.4239 0.2345
## Postura     6  2.905 0.48411  1.7053 0.1226
## Residuals 167 47.410 0.28389
2.Evaluando los supuestos del modelo mixto con el PC2
#para volver una variable categorica:
datos_morfo_med$Tratamiento<-as.factor(datos_morfo_med$Tratamiento)

Re-estimando el modelo con "robust errors"

library(robustlmm)

mixto1pc2_robust<-rlmer(Dim.2~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixto1pc2_robust)
## Robust linear mixed model fit by DAStau 
## Formula: Dim.2 ~ Tratamiento + (1 | Caja) + (1 | Postura) 
##    Data: datos_morfo_med 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5657 -0.6231 -0.0522  0.7214  3.2657 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.00000  0.0000  
##  Postura  (Intercept) 0.07251  0.2693  
##  Residual             0.25711  0.5071  
## Number of obs: 175, groups: Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    0.01424    0.12624   0.113
## TratamientoAD -0.08155    0.09738  -0.838
## TratamientoAE  0.01940    0.09777   0.198
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAD
## TratamintAD -0.410       
## TratamintAE -0.409  0.530
## 
## Robustness weights for the residuals: 
##  146 weights are ~= 1. The remaining 29 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.295   0.519   0.688   0.689   0.870   0.979 
## 
## Robustness weights for the random effects: 
##  27 weights are ~= 1. The remaining one are
##    28 
## 0.783 
## 
## 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(mixto1pc2_robust)
  Dim.2
Predictors Estimates CI p
(Intercept) 0.01 -0.23 – 0.26 0.910
Tratamiento [AD] -0.08 -0.27 – 0.11 0.402
Tratamiento [AE] 0.02 -0.17 – 0.21 0.843
Random Effects
σ2 0.26
τ00 Caja 0.00
τ00 Postura 0.07
ICC 0.22
N Caja 21
N Postura 7
Observations 175
Marginal R2 / Conditional R2 0.006 / 0.225 
plot(mixto1pc2_robust)

Evaluando los residuos del modelo:

*Efecto aleatorio:La postura:

plot(mixto1pc2, resid(., scaled=TRUE) ~ fitted(.) |Postura, abline = 0)

#The following output shows the scaled residuals for each group:
plot(mixto1pc2, 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 (mixto1pc2, Dim.2 ~ fitted(.) | Postura, abline = c (0,1))

*Efecto aleatorio: La caja

plot(mixto1pc2, resid(., scaled=TRUE) ~ fitted(.) |Caja,abline = 0)

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

plot (mixto1pc2, Dim.2 ~ fitted(.) | Caja, abline = c (0,1))

3. Calculando la magnitud del efecto del modelo
library(effectsize)
#proporcion de varianza atribuida a una variable

eta_squared(mixto1pc2)
## Parameter   | Eta2 (partial) |       90% CI
## -------------------------------------------
## Tratamiento |           0.09 | [0.00, 0.33]
library(emmeans)
trat.emm.s <- emmeans(mixto1pc2, "Tratamiento")
pairs(trat.emm.s)
##  contrast estimate    SE   df t.ratio p.value
##  AB - AD    0.1088 0.148 12.2  0.734  0.7484 
##  AB - AE   -0.0471 0.149 12.3 -0.317  0.9463 
##  AD - AE   -0.1559 0.146 11.5 -1.069  0.5510 
## 
## 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(mixto1pc2), edf = 12)
##  contrast effect.size    SE   df lower.CL upper.CL
##  AB - AD       0.1958 0.268 12.9   -0.384    0.775
##  AB - AE      -0.0847 0.268 13.1   -0.663    0.494
##  AD - AE      -0.2805 0.263 12.9   -0.848    0.287
## 
## sigma used for effect sizes: 0.5559 
## Degrees-of-freedom method: inherited from kenward-roger when re-gridding 
## Confidence level used: 0.95
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")
mixto1pc2<-lmer(Dim.2~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
plot_model(mixto1pc2)

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")
mixto1pc2_AB<-lmer(Dim.2~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixto1pc2_AB)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dim.2 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 318.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5608 -0.6316 -0.0924  0.5285  2.9662 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.03867  0.1967  
##  Postura  (Intercept) 0.05519  0.2349  
##  Residual             0.30900  0.5559  
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)    0.01762    0.13865 13.85620   0.127    0.901
## TratamientoAE  0.04709    0.14842 12.36862   0.317    0.756
## TratamientoAD -0.10882    0.14813 12.27026  -0.735    0.476
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAE
## TratamintAE -0.551       
## TratamintAD -0.552  0.516
plot_model(mixto1pc2_AB)

AB y AE
#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")
mixto1pc2_AD<-lmer(Dim.2~Tratamiento+(1|Caja)+(1|Postura),data=datos_morfo_med)
summary(mixto1pc2_AD)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dim.2 ~ Tratamiento + (1 | Caja) + (1 | Postura)
##    Data: datos_morfo_med
## 
## REML criterion at convergence: 318.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5608 -0.6316 -0.0924  0.5285  2.9662 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Caja     (Intercept) 0.03867  0.1967  
##  Postura  (Intercept) 0.05519  0.2349  
##  Residual             0.30900  0.5559  
## Number of obs: 175, groups:  Caja, 21; Postura, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)    -0.0912     0.1359 12.8588  -0.671    0.514
## TratamientoAB   0.1088     0.1481 12.2703   0.735    0.476
## TratamientoAE   0.1559     0.1459 11.6044   1.069    0.307
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmAB
## TratamintAB -0.526       
## TratamintAE -0.534  0.490
plot_model(mixto1pc2_AD)

Probando otras opciones para ver magnitud del efecto:

library(effectsize)
effectsize(mixto1pc2)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |               0.10 | [-0.32, 0.53]
## TratamientoAB |              -0.08 | [-0.54, 0.39]
## TratamientoAD |              -0.25 | [-0.71, 0.21]
## 
## # Standardization method: refit
effectsize(mixto1pc2_AB)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |               0.03 | [-0.41, 0.46]
## TratamientoAE |               0.08 | [-0.39, 0.54]
## TratamientoAD |              -0.17 | [-0.64, 0.29]
## 
## # Standardization method: refit
effectsize(mixto1pc2_AD)
## Parameter     | Coefficient (std.) |        95% CI
## --------------------------------------------------
## (Intercept)   |              -0.15 | [-0.57, 0.28]
## TratamientoAB |               0.17 | [-0.29, 0.64]
## TratamientoAE |               0.25 | [-0.21, 0.71]
## 
## # Standardization method: refit
datos_morfo_med$Tratamiento <-relevel(datos_morfo_med$Tratamiento, ref="AE")
ggplot(datos_morfo_med, aes(x=Postura, y=Dim.2, 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