Continuar utilizando las funciones del conjunto de librerías tidyverse para la manipulación y análisis de datos, y la librería flextable para la presentación de resultados.
Usar la función str()
para explorar la estructura de un objeto en R.
Aprender a acceder a las variables contenidas en data.frames
y tibbles
.
Conocer el concepto de lista en R, y aprender a acceder a sus componentes.
Utilizar BinomCI
del paquete DescTools
para calcular intervalos de confianza para proporciones.
Utilizar t.test
para calcular intervalos de confianza para medias.
Utilizar varTest
para calcular intervalos de confianza para varianzas y desviaciones típicas.
Conocer la función do()
para integrar el cálculo de intervalos de confianza en cadenas de comandos de tidyverse
enlazados por tuberías
Aprender a construir tablas resumen con estadísticos e intervalos de confianza para varias variables.
library(readxl)
library(tidyverse)
library(flextable)
library(janitor)
library(DescTools)
library(broom)
library(EnvStats)
Como ya hemos visto en prácticas anteriores, un data.frame
o un tibble
son las estructuras que utiliza R normalmente para archivar bases de datos. En ambos casos se trata de estructuras rectangulares en forma de matriz, donde cada fila representa un objeto y cada columna una variable. Las columnas (variables) tienen nombres diferentes, y los tipos de datos pueden ser distintos para cada columna: una columna puede contener datos numéricos, otra contener fechas, una tercera contener caracteres, etc.
Comenzamos por cargar el archivo de datos que hemos utilizado en las prácticas anteriores, correspondiente a un muestreo de nidos de tortugas en la isla de Boa Vista del archipiélago de Cabo Verde:
NOTA Codificación de caracteres en un archivo (fileEncoding
): en elcomando anterior hemos utilizado la opción fileEncoding="UTF-8"
para que R interprete correctamente los caracteres típicos del español: la ñ, las tildes, diéresis, etc. Ello es debido a que el archivo de las tortugas se codificó originalmente en un sistema operativo Linux, que usa esa codificación específica de caracteres; si no se especifica “UTF8” y el archivo se lee en un sistema Windows, los caracteres citados no serán legibles. Si la lectura se hiciera en un Mac no habría problemas, ya que Mac utiliza también codificación UTF8. Si un usuario de un Mac quisiera leer archivos que contienen caracteres codificados en Windows, debería usar la opción fileEncoding="latin1"
.
La función str()
nos da información sobre el tipo de objeto que es tortugas
, y sobre cuál es su contenido:
## 'data.frame': 1277 obs. of 13 variables:
## $ Año : int 1999 1999 1999 1999 1999 1999 1999 1999 1999 1999 ...
## $ periodo : int 1 1 1 1 1 1 1 1 1 1 ...
## $ LCC : num 77.4 82 87.3 85.4 79.2 80.7 78.7 85.7 80.4 89.4 ...
## $ ACC : num 71.6 76.3 85.5 79.7 74.4 74.5 75.9 78.6 73.3 81.6 ...
## $ peso : num 55.7 54.4 66.5 67.4 56.2 64.1 56.2 60.4 55.9 66.2 ...
## $ Huevos : int 81 99 127 107 101 77 104 105 96 108 ...
## $ playa : Factor w/ 4 levels "Calheta","Ervatao",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ distancia : num 28 35 28.8 29.3 41.2 42.8 15.6 24.8 20.7 22.1 ...
## $ profNido : num 61.3 59.5 60.2 60.6 69.5 74.3 49.2 57.1 60.1 64.4 ...
## $ crias_Emerg : int 65 30 NA 9 0 NA 15 NA 54 NA ...
## $ crias_Muertas: int 3 6 NA 2 1 NA 1 NA 1 NA ...
## $ hrotos : int 13 63 NA 96 100 NA 88 NA 41 NA ...
## $ cangrejos : int 0 1 NA 1 1 NA 1 NA 0 NA ...
Como podemos observar, en la primera linea se nos indica que tortugas
es un objeto de tipo data.frame
con 1277 observaciones de 13 variables, cuyos nombres figuran a continuación. Para cada variable se indica su tipo (numérico, factor, entero) y se muestran por defecto sus 10 primeros valores.
Para extraer los valores de una variable concreta de este data.frame
, por ejemplo los valores de LCC, podríamos usar cualquiera de las siguientes expresiones equivalentes: tortugas$LCC
, tortugas[["LCC"]]
o tortugas[,"LCC"]
, como podemos ver a continuación:
## [1] 77.4 82.0 87.3 85.4 79.2 80.7 78.7 85.7 80.4 89.4 84.3 80.0 79.2 80.7
## [15] 91.9 84.3 86.1 77.6 78.5 88.5 82.2 86.3 89.3 81.7 81.9 81.8 82.2 78.7
## [1] 77.4 82.0 87.3 85.4 79.2 80.7 78.7 85.7 80.4 89.4 84.3 80.0 79.2 80.7
## [15] 91.9 84.3 86.1 77.6 78.5 88.5 82.2 86.3 89.3 81.7 81.9 81.8 82.2 78.7
Ello nos permitiría, por ejemplo, calcular rápidamente la media (o cualquier otro estadístico descriptivo) de esta variable:
## [1] 81.91457
sin usar las cadenas de tuberías típicas del ámbito del tidyverse
.
Cuando se usa Rstudio es más cómodo usar la notación del $
, ya que una vez que escribimos tortugas$
, Rstudio nos muestra las distintas variables en una caja para que podamos elegir el nombre de la que queremos extraer:
De modo alternativo, en lugar de leer el archivo de tortugas con la función read.csv2()
, podemos leerlo con la función equivalente read_csv2()
(nótese que hemos sustituido el punto por la barra baja) del paquete tidyverse
; la sintaxis para especificar la codificación de caracteres del archivo es ligeramente distinta:
Nuevamente, si analizamos la estructura del objeto tortugas
mediante la función str()
obtenemos:
## tibble [1,277 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Año : num [1:1277] 1999 1999 1999 1999 1999 ...
## $ periodo : num [1:1277] 1 1 1 1 1 1 1 1 1 1 ...
## $ LCC : num [1:1277] 77.4 82 87.3 85.4 79.2 80.7 78.7 85.7 80.4 89.4 ...
## $ ACC : num [1:1277] 71.6 76.3 85.5 79.7 74.4 74.5 75.9 78.6 73.3 81.6 ...
## $ peso : num [1:1277] 55.7 54.4 66.5 67.4 56.2 64.1 56.2 60.4 55.9 66.2 ...
## $ Huevos : num [1:1277] 81 99 127 107 101 77 104 105 96 108 ...
## $ playa : chr [1:1277] "Ervatao" "Ervatao" "Ervatao" "Ervatao" ...
## $ distancia : num [1:1277] 28 35 28.8 29.3 41.2 42.8 15.6 24.8 20.7 22.1 ...
## $ profNido : num [1:1277] 61.3 59.5 60.2 60.6 69.5 74.3 49.2 57.1 60.1 64.4 ...
## $ crias_Emerg : num [1:1277] 65 30 NA 9 0 NA 15 NA 54 NA ...
## $ crias_Muertas: num [1:1277] 3 6 NA 2 1 NA 1 NA 1 NA ...
## $ hrotos : num [1:1277] 13 63 NA 96 100 NA 88 NA 41 NA ...
## $ cangrejos : num [1:1277] 0 1 NA 1 1 NA 1 NA 0 NA ...
## - attr(*, "spec")=
## .. cols(
## .. Año = col_double(),
## .. periodo = col_double(),
## .. LCC = col_double(),
## .. ACC = col_double(),
## .. peso = col_double(),
## .. Huevos = col_double(),
## .. playa = col_character(),
## .. distancia = col_double(),
## .. profNido = col_double(),
## .. crias_Emerg = col_double(),
## .. crias_Muertas = col_double(),
## .. hrotos = col_double(),
## .. cangrejos = col_double()
## .. )
Como vemos, al leer los datos con read_csv2()
la estructura de los datos es ligeramente distinta; la primera linea del resultado anterior nos indica que tortugas
es un objeto de varias clases: Classes ‘spec_tbl_df’, ‘tbl_df’, ‘tbl’ and ‘data.frame’, con 1277 observaciones y 13 variables. El acrónimo tbl
se refiere a la palabra tibble
que suena parecido a table
(tabla). Los objetos de clase tibble
pueden ser algo más complejos que los data.frames
, si bien al nivel que los utilizaremos en este curso, podemos considerar que son iguales: tablas rectangulares de datos.
Las únicas diferencias en la estructura de datos según se lean con read.csv2()
o con read_csv2()
, son que todas las variables numéricas son del mismo tipo (no se distinguen variables enteras de no enteras), y las variables tipo carácter no se convierten automáticamente en factores (obsérvese que la variable playa
al leerla con read.csv2()
se ha leido como factor, mientras que al leerla con read_csv2()
se ha leido como carácter). Estas diferencias, en cualquier caso no son relevantes.
Con los tibbles
también podemos usar la sintaxis tortugas$LCC
o tortugas[[“LCC”]] para extraer los valores de una variable concreta.
Las listas constituyen otro tipo de estructura de datos que utiliza R para guardar información relacionada que no tiene por qué tener necesariamente una estructura rectangular como una tabla. Muchas funciones de R para el análisis estadístico devuelven sus resultados en forma de lista.
Podemos imaginar una lista como un contenedor de objetos de distinta naturaleza. En el ejemplo siguiente definimos una lista que contiene unos datos experimentales medidos en tres sujetos, el nombre del laboratorio que hizo las medidas, la fecha de su realización, y si el resultado para cada sujeto es positivo o negativo:
analiticas <- list(datos = data.frame(sujeto1=c(10,12,11.2,9.8),
sujeto2=c(5.4,6.9,7.8,4.6),
sujeto3=c(8.7,7.5,6.2,0.3)),
laboratorio="X Labs SL",
fecha="23-08-2019",
resultado=c("positivo","negativo","positivo"))
Como vemos, dentro de la lista analíticas
hemos introducido un data.frame
con datos, dos variables con un único valor (laboratorio
y fecha
), y una tercera variable (resultado
) que es un vector con 3 valores de tipo carácter. Si pedimos a R que nos enseñe en contenido de esta lista nos lo muestra de la forma siguiente:
## $datos
## sujeto1 sujeto2 sujeto3
## 1 10.0 5.4 8.7
## 2 12.0 6.9 7.5
## 3 11.2 7.8 6.2
## 4 9.8 4.6 0.3
##
## $laboratorio
## [1] "X Labs SL"
##
## $fecha
## [1] "23-08-2019"
##
## $resultado
## [1] "positivo" "negativo" "positivo"
Y si le pedimos que nos muestre la estructura, R nos indica que analiticas
es una lista y nos muestra de qué clase en cada uno de los elementos que contiene:
## List of 4
## $ datos :'data.frame': 4 obs. of 3 variables:
## ..$ sujeto1: num [1:4] 10 12 11.2 9.8
## ..$ sujeto2: num [1:4] 5.4 6.9 7.8 4.6
## ..$ sujeto3: num [1:4] 8.7 7.5 6.2 0.3
## $ laboratorio: chr "X Labs SL"
## $ fecha : chr "23-08-2019"
## $ resultado : chr [1:3] "positivo" "negativo" "positivo"
Del mismo modo que ya hemos visto con data.frames
y tibbles
podemos extraer una parte de una lista mediante el uso del símbolo $
o el doble corchete:
## [1] "X Labs SL"
## [1] "positivo" "negativo" "positivo"
Es conveniente conocer el concepto de lista en R, así como la manera de acceder a su contenido, ya que R devuelve el resultado de muchos procedimientos de inferencia estadística en forma de lista, en particular los que vamos a utilizar aquí para calcular intervalos de confianza, así como los que utilizaremos más adelante en contrastes de hipótesis y en modelos lineales.
Construir un intervalo de confianza al 95% para la proporción de nidos afectados por cangrejos.
Para construir un intervalo de confianza para una proporción utilizamos la función BinomCI
de la librería DescTools
. Para aplicar esta función necesitamos saber cuántos nidos se revisaron y en cuántos había cangrejos. Para ello podemos usar la función tabyl
aplicada a la variable cangrejos
.
Como no todos los nidos de Boa Vista se revisaron para ver si tenían cangrejos, eliminamos en primer lugar todos los valores perdidos de la variable cangrejos (que están codificados como NA) y a continuación aplicamos la función tabyl
, calculando una fila adicional con los totales:
## cangrejos n percent
## 0 305 0.3388889
## 1 595 0.6611111
## Total 900 1.0000000
Por tanto, vemos que en total se revisaron 900 nidos de los cuales tenían cangrejos 595; la proporción estimada de nidos con cangrejos es entonces:
\[\hat{p}=\frac{595}{900}=0.6611\] que corresponde al porcentaje calculado por tabyl
en la última columna. Ahora, para calcular el intervalo de confianza usamos la función BinomCI
, que recibe como argumentos el número de nidos con cangrejos, el número total de nidos revisados y el nivel de confianza:
## est lwr.ci upr.ci
## [1,] 0.6611111 0.6295608 0.691292
Por tanto, podemos asegurar con un 95% de confianza que de todos los nidos de la isla, entre un 62.95% y un 69.12% están afectados por la presencia de cangrejos.
Si queremos incluir el cálculo del intervalo de confianza en una cadena de comandos enlazados por tuberías para conseguir una mejor presentación final del resultado, podemos empezar por calcular el número de éxitos (número de nidos con cangrejos) y el número total de observaciones que necesita la función BinomCI
. Téngase en cuenta que como la variable cangrejos
vale 0 cuando no hay cangrejos y 1 cuando sí hay, la suma de valores de esta variable es igual al número total de nidos que tienen cangrejos.
## # A tibble: 1 x 2
## nE n
## <dbl> <int>
## 1 595 900
Obtenemos así solo los valores que necesitamos para BinomCI()
, con una presentación más simple que la tabla calculada por tabyl
. Ahora, para poder enlazar BinomCI()
con este resultado, debemos usar la función do()
. Esta función permite incluir en una cadena de tuberías el cálculo de cualquier otra función; se deben tener en cuenta las siguientes restricciones:
Si la función que se incluye usa variables procedentes de la “tubería” anterior, dichas variables deben ir precedidas por el símbolo “.$
”.
El resultado de la función debe ser un data.frame
o un tibble
. Podemos averiguar de qué clase es el resultado de la función BinomCI()
mirando su estructura mediante str()
:
## num [1, 1:3] 0.661 0.63 0.691
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "est" "lwr.ci" "upr.ci"
Como vemos, el resultado es una matriz numérica de una fila y una columna. Para convertir una matriz m
en un data.frame utilizamos data.frame(m)
. Por tanto, podemos obtener el resultado de BinomCI()
al final de nuestra cadena de tuberías mediante:
tortugas %>%
filter(!is.na(cangrejos)) %>%
summarize(nE=sum(cangrejos),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95)))
## est lwr.ci upr.ci
## 1 0.6611111 0.6295608 0.691292
Por último, podemos mejorar la presentación final creando una función que le dé formato al estimador y al intervalo de confianza (esta función es genérica y la podremos utilizar posteriormente con otros intervalos). IMPORTANTE: Si se va a utilizar esta función, es preciso incluirla en el documento Rmd antes de llamarla por primera vez)
# La función recibe como argumentos, el estimador del parámetro y los límites
# inferior y superior del intervalo de confianza. Por defecto se redondean los
# valores a 4 dígitos y se considera que el nivel de confianza el 0.95. El valor
# del nivel de confianza se usa solamente para poner nombre a la variable que
# designa el intervalo, como suele ser habitual en publicaciones científicas
# La función devuelve un data.frame con el estimador y el intervalo formateados
# con el número de dígitos indicado.
formatInterval <- function(est, inf, sup, digits=4, conf.level=0.95){
CI <- sprintf("(%s, %s)",round(inf,digits),round(sup,digits))
estimate <- round(est,digits)
df <- data.frame(estimate,CI, stringsAsFactors = FALSE)
names(df)[2]=paste("CI",100*conf.level,"%",sep="")
return(df)
}
Por último integramos todo en un cadena de “tuberías” y generamos la salida con flextable
:
tortugas %>%
filter(!is.na(cangrejos)) %>%
summarize(nE=sum(cangrejos),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
autofit()
estimate | CI95% |
0.6611 | (0.6296, 0.6913) |
Construir un intervalo de confianza al 95% para la proporción de nidos que tienen más de 100 huevos.
Primero definimos una variable que vale 1 si el nido tiene más de 100 huevos y 0 en caso contrario:
## masDe100 n percent
## 0 1166 0.91307753
## 1 111 0.08692247
## Total 1277 1.00000000
Por tanto, la proporción estimada de nidos con más de 100 huevos es el 8.69%. Calculamos el intervalo de confianza con BinomCI()
:
## est lwr.ci upr.ci
## [1,] 0.08692247 0.07268332 0.1036394
Por tanto, con un 95% de confianza podemos asegurar que entre un 7.27% y un 10.37% de los nidos de tortuga de Boa Vista tienen más de 100 huevos.
Si queremos el intervalo con mejor presentación:
tortugas %>%
mutate(masDe100=ifelse(Huevos>100,1,0)) %>%
summarize(nE=sum(masDe100),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
autofit()
estimate | CI95% |
0.0869 | (0.0727, 0.1036) |
Construir un intervalo de confianza al 95% para la proporción de tortugas que pesan más de 70 kg.
Creamos en primer lugar una variable que vale 1 si la tortuga pesa más de 70 kg y 0 en caso contrario:
## masDe70 n percent
## 0 1232 0.96476116
## 1 45 0.03523884
## Total 1277 1.00000000
Por tanto, la proporción estimada de tortugas que pesan más de 70 kg es el 3.52%. Calculamos el intervalo de confianza con BinomCI()
:
## est lwr.ci upr.ci
## [1,] 0.03523884 0.02643932 0.04682616
Por tanto, con un 95% de confianza podemos asegurar que entre un 2.64% y un 4.68% de las tortugas que anidan en Boa Vista pesan más de 70 kg.
Si queremos el intervalo con mejor presentación:
tortugas %>%
mutate(masDe70=ifelse(peso>70,1,0)) %>%
summarize(nE=sum(masDe70),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
autofit()
estimate | CI95% |
0.0352 | (0.0264, 0.0468) |
Construir una tabla resumen con las proporciones e intervalos de confianza calculados en los tres ejercicios anteriores.
Podemos construir una tabla resumen repitiendo los cálculos anteriores y “pegando” los resultados obtenidos mediante la función bind_rows()
(que significa algo así como pegar por filas):
ejerc1 <- tortugas %>%
filter(!is.na(cangrejos)) %>%
summarize(nE=sum(cangrejos),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci))
ejerc2 <- tortugas %>%
mutate(masDe100=ifelse(Huevos>100,1,0)) %>%
summarize(nE=sum(masDe100),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci))
ejerc3 <- tortugas %>%
mutate(masDe70=ifelse(peso>70,1,0)) %>%
summarize(nE=sum(masDe70),n=n()) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci))
tabla <- bind_rows(ejerc1,ejerc2,ejerc3)
tabla
## estimate CI95%
## 1 0.6611 (0.6296, 0.6913)
## 2 0.0869 (0.0727, 0.1036)
## 3 0.0352 (0.0264, 0.0468)
Ahora añadimos a la tabla una primera columna con los nombres de las variables, para poder identificar a qué corresponde cada proporción, y presentamos el resultado con flextable
:
tabla %>%
mutate(variable=c("Nidos con cangrejos","Nidos con más de 100 huevos",
"Tortugas de más de 70 kg")) %>%
select(variable, estimate, `CI95%`) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
set_caption("Proporciones estimadas e Intervalos de Confianza al 95%") %>%
autofit()
variable | estimate | CI95% |
Nidos con cangrejos | 0.6611 | (0.6296, 0.6913) |
Nidos con más de 100 huevos | 0.0869 | (0.0727, 0.1036) |
Tortugas de más de 70 kg | 0.0352 | (0.0264, 0.0468) |
También podía haberse construido la tabla de una manera más eficiente utilizando pivot_longer()
(ver sección 2 en el documento de complementos a las prácticas en http://estadistica-dma.ulpgc.es/estadFCM/complementos-Practicas.html)
tortugas %>%
mutate(masDe100=ifelse(Huevos>100,1,0),
masDe70=ifelse(peso>70,1,0)) %>%
select(cangrejos,masDe100,masDe70) %>%
pivot_longer(1:3,names_to="Variable",values_to = "Valor") %>%
group_by(Variable) %>%
summarize(nE=sum(Valor,na.rm=TRUE),n=length(which(!is.na(Valor)))) %>%
do(data.frame(BinomCI(.$nE,.$n,conf.level=0.95))) %>%
do(formatInterval(.$est,.$lwr.ci,.$upr.ci)) %>%
mutate(Variable=c("Nidos con cangrejos","Nidos con más de 100 huevos",
"Tortugas de más de 70 kg")) %>%
select(Variable, estimate, `CI95%`) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
set_caption("Proporciones estimadas e Intervalos de Confianza al 95%") %>%
autofit()
Variable | estimate | CI95% |
Nidos con cangrejos | 0.6611 | (0.6296, 0.6913) |
Nidos con más de 100 huevos | 0.0869 | (0.0727, 0.1036) |
Tortugas de más de 70 kg | 0.0352 | (0.0264, 0.0468) |
Construir un intervalo de confianza al 95% para la longitud media del caparazón de las tortugas de Boa Vista.
Para construir un intervalo de confianza para la media utilizamos la función t.test
:
##
## One Sample t-test
##
## data: tortugas$LCC
## t = 634.59, df = 1276, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 81.66133 82.16780
## sample estimates:
## mean of x
## 81.91457
Esta salida nos informa de que el valor medio de la variable LCC
observado en la muestra ha sido de 81.91457 cm, y por tanto este es el valor que estimamos para la media de la población; el intervalo de confianza al 95% es (81.66133, 82.16780), lo que nos indica que podemos estar seguros en un 95% de que el valor medio de LCC
en la población de tortugas que anida en Boa Vista está comprendido entre 81.66 y 82.17 cm.
Como vemos, esta función devuelve más cosas aparte del intervalo de confianza; ello es así porque la misma función nos sirve también para realizar contrastes de hipótesis sobre la media (lo veremos en el siguiente capítulo). Podemos analizar la estructura del objeto que devuelve t.test
mediante la función str()
:
## List of 10
## $ statistic : Named num 635
## ..- attr(*, "names")= chr "t"
## $ parameter : Named num 1276
## ..- attr(*, "names")= chr "df"
## $ p.value : num 0
## $ conf.int : num [1:2] 81.7 82.2
## ..- attr(*, "conf.level")= num 0.95
## $ estimate : Named num 81.9
## ..- attr(*, "names")= chr "mean of x"
## $ null.value : Named num 0
## ..- attr(*, "names")= chr "mean"
## $ stderr : num 0.129
## $ alternative: chr "two.sided"
## $ method : chr "One Sample t-test"
## $ data.name : chr "tortugas$LCC"
## - attr(*, "class")= chr "htest"
Así pues, el resultado de t.test
es una lista, en la que podemos ver que hay una componente (estimate
) en el que se muestra el valor medio de la variable LCC (81.9 cm) y otra componente (conf.int
) que contiene el intervalo de confianza. Podemos extraer estas componentes sencillamente mediante.
## mean of x
## 81.91457
## [1] 81.66133 82.16780
## attr(,"conf.level")
## [1] 0.95
No obstante, para poder utilizar estos valores dentro de una cadena de tuberías de tidyverse
lo más cómodo es convertir la lista producida por t.test en un data.frame
. La función tidy
de la librería broom
se encarga precisamente de esta tarea:
## # A tibble: 1 x 8
## estimate statistic p.value parameter conf.low conf.high method
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 81.9 635. 0 1276 81.7 82.2 One S…
## # … with 1 more variable: alternative <chr>
Como solo nos interesan la media y el intervalo de confianza, podemos usar nuestra función formatInterval
(téngase en cuenta que el nombre que pone tidy
a los extremos del intervalo de confianza son, como puede verse en la salida anterior, conf.low
y conf.high
):
tortugas %>%
do(tidy(t.test(.$LCC,conf.level=0.95))) %>%
do(formatInterval(.$estimate, .$conf.low, .$conf.high)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
autofit()
estimate | CI95% |
81.9146 | (81.6613, 82.1678) |
Construir intervalos de confianza al 99% para los valores medios de:
ACC: Anchura curva del caparazón.
profNido: profundidad del nido.
distancia: distancia desde el nido hasta la orilla.
Huevos: número de huevos en cada nido
Para resolver este problema podemos ir variable a variable:
## mean of x
## 77.18661
## [1] 76.89855 77.47467
## attr(,"conf.level")
## [1] 0.99
## mean of x
## 48.06335
## [1] 47.32201 48.80469
## attr(,"conf.level")
## [1] 0.99
## mean of x
## 24.00799
## [1] 23.48184 24.53414
## attr(,"conf.level")
## [1] 0.99
## mean of x
## 83.00313
## [1] 82.09878 83.90749
## attr(,"conf.level")
## [1] 0.99
Podemos usar las herramientas del tidyverse
en combinación con t.test
para sacar todos los resultados en una única tabla:
tortugas %>%
select(LCC, ACC, profNido, distancia, Huevos) %>%
pivot_longer(1:5,names_to="Variable",values_to = "Valor") %>%
group_by(Variable) %>%
do(tidy(t.test(.$Valor, conf.level=0.99))) %>%
do(formatInterval(.$estimate,.$conf.low,.$conf.high,digits=2,conf.level = 0.99)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
set_caption("Medias estimadas e Intervalos de Confianza al 99%") %>%
autofit()
Variable | estimate | CI99% |
ACC | 77.19 | (76.9, 77.47) |
distancia | 24.01 | (23.48, 24.53) |
Huevos | 83.00 | (82.1, 83.91) |
LCC | 81.91 | (81.58, 82.25) |
profNido | 48.06 | (47.32, 48.8) |
Construir los intervalos de confianza del ejercicio anterior para cada año.
Para construir un intervalo de confianza para cada año bastará simplemente con añadir Año
en el group_by()
. Para ello debemos incluir Año
en la lista de variables que seleccionamos inicialmente:
tortugas %>%
select(LCC, ACC, profNido, distancia, Huevos,Año) %>%
pivot_longer(1:5,names_to="Variable",values_to = "Valor") %>%
group_by(Año,Variable) %>%
do(tidy(t.test(.$Valor, conf.level=0.99))) %>%
do(formatInterval(.$estimate,.$conf.low,.$conf.high,digits=2,conf.level = 0.99)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
set_caption("Medias estimadas e Intervalos de Confianza al 99%") %>%
autofit()
Año | Variable | estimate | CI99% |
1999 | ACC | 77.25 | (76.55, 77.95) |
1999 | distancia | 23.92 | (22.62, 25.22) |
1999 | Huevos | 92.80 | (90.54, 95.07) |
1999 | LCC | 82.15 | (81.35, 82.95) |
1999 | profNido | 49.18 | (47.32, 51.03) |
2000 | ACC | 77.41 | (76.72, 78.09) |
2000 | distancia | 23.65 | (22.51, 24.79) |
2000 | Huevos | 82.12 | (80.26, 83.98) |
2000 | LCC | 82.45 | (81.68, 83.23) |
2000 | profNido | 48.26 | (46.54, 49.99) |
2001 | ACC | 76.99 | (76.24, 77.73) |
2001 | distancia | 24.41 | (22.96, 25.86) |
2001 | Huevos | 77.19 | (75.04, 79.34) |
2001 | LCC | 81.89 | (81.01, 82.77) |
2001 | profNido | 47.37 | (45.42, 49.33) |
2002 | ACC | 76.71 | (76.02, 77.4) |
2002 | distancia | 24.05 | (22.75, 25.35) |
2002 | Huevos | 81.73 | (79.73, 83.72) |
2002 | LCC | 81.03 | (80.19, 81.87) |
2002 | profNido | 48.33 | (46.49, 50.18) |
2003 | ACC | 77.40 | (76.63, 78.17) |
2003 | distancia | 23.71 | (22.37, 25.04) |
2003 | Huevos | 84.91 | (82.79, 87.02) |
2003 | LCC | 81.92 | (81.05, 82.79) |
2003 | profNido | 47.34 | (45.5, 49.18) |
2004 | ACC | 77.32 | (76.64, 78) |
2004 | distancia | 24.33 | (23.04, 25.63) |
2004 | Huevos | 79.87 | (77.82, 81.92) |
2004 | LCC | 81.99 | (81.22, 82.75) |
2004 | profNido | 47.86 | (46.08, 49.63) |
Podemos mejorar un poco la presentación de esta tabla si añadimos unas lineas horizontales entre los datos de los distintos años; esto se consigue con la función hline
del paquete officer
. Si además queremos evitar que se repita varias veces el valor de Año
en la tabla, usamos merge_v()
también del paquete officer
como se muestra a continuación:
library(officer)
tortugas %>%
select(LCC, ACC, profNido, distancia, Huevos,Año) %>%
pivot_longer(1:5,names_to="Variable",values_to = "Valor") %>%
group_by(Año,Variable) %>%
do(tidy(t.test(.$Valor, conf.level=0.99))) %>%
do(formatInterval(.$estimate,.$conf.low,.$conf.high,digits=2,conf.level = 0.99)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
set_caption("Medias estimadas e Intervalos de Confianza al 99%") %>%
hline(i = seq(5,25,by=5), part = "body",
border = fp_border(color = "grey70", width = 2) ) %>%
merge_v(j = ~Año) %>%
fix_border_issues() %>%
autofit()
Año | Variable | estimate | CI99% |
1999 | ACC | 77.25 | (76.55, 77.95) |
distancia | 23.92 | (22.62, 25.22) | |
Huevos | 92.80 | (90.54, 95.07) | |
LCC | 82.15 | (81.35, 82.95) | |
profNido | 49.18 | (47.32, 51.03) | |
2000 | ACC | 77.41 | (76.72, 78.09) |
distancia | 23.65 | (22.51, 24.79) | |
Huevos | 82.12 | (80.26, 83.98) | |
LCC | 82.45 | (81.68, 83.23) | |
profNido | 48.26 | (46.54, 49.99) | |
2001 | ACC | 76.99 | (76.24, 77.73) |
distancia | 24.41 | (22.96, 25.86) | |
Huevos | 77.19 | (75.04, 79.34) | |
LCC | 81.89 | (81.01, 82.77) | |
profNido | 47.37 | (45.42, 49.33) | |
2002 | ACC | 76.71 | (76.02, 77.4) |
distancia | 24.05 | (22.75, 25.35) | |
Huevos | 81.73 | (79.73, 83.72) | |
LCC | 81.03 | (80.19, 81.87) | |
profNido | 48.33 | (46.49, 50.18) | |
2003 | ACC | 77.40 | (76.63, 78.17) |
distancia | 23.71 | (22.37, 25.04) | |
Huevos | 84.91 | (82.79, 87.02) | |
LCC | 81.92 | (81.05, 82.79) | |
profNido | 47.34 | (45.5, 49.18) | |
2004 | ACC | 77.32 | (76.64, 78) |
distancia | 24.33 | (23.04, 25.63) | |
Huevos | 79.87 | (77.82, 81.92) | |
LCC | 81.99 | (81.22, 82.75) | |
profNido | 47.86 | (46.08, 49.63) |
Nótese que para colocar las lineas hemos utilizado la función seq(5,25,by=5)
que selecciona los números del 5 al 25 de 5 en 5 (es decir, selecciona los números 5, 10, 15, 20 y 25), que son las posiciones donde queremos colocar las lineas separadoras horizontales.
Construir un intervalo de confianza al 95% para la desviación típica de la longitud del caparazón de las tortugas de Boa Vista.
El procedimiento es muy similar al del cálculo del intervalo de confianza para la media, con la única diferencia de que ahora deberemos usar la función varTest()
de la librería EnvStats
:
##
## Chi-Squared Test on Variance
##
## data: tortugas$LCC
## Chi-Squared = 27151, df = 1276, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1
## 95 percent confidence interval:
## 19.71889 23.03070
## sample estimates:
## variance
## 21.27805
Téngase en cuenta que esta librería proporciona el estimador de la varianza y su intervalo de confianza; el estimador de la desviación típica y su intervalo de confianza pueden obtenerse simplemente hallando la raiz cuadrada:
## variance
## 4.612813
(aunque en realidad nos muestra la desviación típica, R sigue escribiendo variance porque ese es el nombre que tenía originalmente el valor de vtLCC$estimate
, y la función raiz cuadrada sqrt()
no altera ese nombre). Para eliminar el nombre podemos utilizar as.numeric()
:
## [1] 4.612813
## LCL UCL
## 4.440595 4.799031
## attr(,"conf.level")
## [1] 0.95
Como en el caso del t.test podemos tener una mejor presentación usando las herramientas del tidyverse
:
## # A tibble: 1 x 7
## estimate statistic p.value conf.low.LCL conf.high.UCL method alternative
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 21.3 27151. 0 19.7 23.0 Chi-Sq… two.sided
Nótese que los nombres para los extremos del intervalo son ahora conf.low.LCL
y conf.high.UCL
. Extraemos la varianza y los extremos del intervalo, hacemos su raiz cuadrada y pasamos los valores a formatInterval
y flextable
:
tortugas %>%
do(tidy(varTest(.$LCC,conf.level=0.95))) %>%
do(formatInterval(sqrt(.$estimate), sqrt(.$conf.low.LCL), sqrt(.$conf.high.UCL))) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
autofit()
estimate | CI95% |
4.6128 | (4.4406, 4.799) |
Construir intervalos de confianza al 99% para las desviaciones típicas de las variables:
ACC: Anchura curva del caparazón.
profNido: profundidad del nido.
distancia: distancia desde el nido hasta la orilla.
Huevos: número de huevos en cada nido
Para resolver este problema podemos ir variable a variable:
## [1] 3.990378
## LCL UCL
## 3.841399 4.151469
## attr(,"conf.level")
## [1] 0.95
o mejor, construir una tabla:
tortugas %>%
select(LCC, ACC, profNido, distancia, Huevos) %>%
pivot_longer(1:5,names_to="Variable",values_to = "Valor") %>%
group_by(Variable) %>%
do(tidy(varTest(.$Valor, conf.level=0.99))) %>%
do(formatInterval(sqrt(.$estimate),sqrt(.$conf.low.LCL),sqrt(.$conf.high.UCL),
digits=2,conf.level = 0.99)) %>%
flextable() %>%
fontsize(size=14, part="all") %>%
set_caption("Desviaciones típicas estimadas e Intervalos de Confianza al 99%") %>%
autofit()
Variable | estimate | CI99% |
ACC | 3.99 | (3.8, 4.2) |
distancia | 7.29 | (6.93, 7.68) |
Huevos | 12.53 | (11.92, 13.2) |
LCC | 4.61 | (4.39, 4.86) |
profNido | 10.27 | (9.77, 10.82) |