Ejercicio 1.a

Dados los datos de la tabla calcular las medidas de centralización y de dispersión.

Lo primero cargamos y mostramos los datos:

library(readr)
tabla1 <- read_csv("tabla1.txt")
tabla1$Alumnos <- c(1:nrow(tabla1))
tabla1 <- tabla1[,c(2,1)]
Primer set de datos
Alumnos Estatura
1 1,2
2 1,3
3 1,3
4 1,2
5 1,2
6 1,3
7 1,3
8 1,2
9 1,3
10 1,3
11 1,2
12 1,3
13 1,3
14 1,2
15 1,3
16 1,3
17 1,2
18 1,2
19 1,2
20 1,3
21 1,2
22 1,3
23 1,3
24 1,2
25 1,3
26 1,3
27 1,3
28 1,2
29 1,2
30 1,2

Y ahora calculamos las medidas pedidas. Las medidas de posición centrales y no centrales:

media_arit <- mean(tabla1$Estatura)
cat("La media aritmética del primer set de datos es:", media_arit)
## La media aritmética del primer set de datos es: 1.253333
media_geom <- prod(tabla1$Estatura)^(1/nrow(tabla1))
cat("La media geométrica del primer set de datos es:", media_geom)
## La media geométrica del primer set de datos es: 1.252927
middle <- median(tabla1$Estatura)
cat("La mediana del primer set de datos es:", middle)
## La mediana del primer set de datos es: 1.26
max.freq <- max(table(tabla1$Estatura))
max.freq.names <- names(which(table(tabla1$Estatura) == max.freq))
cat("Las modas de los datos son", max.freq.names, "repitiéndose", max.freq, "veces")
## Las modas de los datos son 1.21 1.22 1.28 repitiéndose 4 veces
cuartiles <- quantile(tabla1$Estatura)
show(cuartiles)
##   0%  25%  50%  75% 100% 
## 1.20 1.22 1.26 1.28 1.30
deciles <- quantile(tabla1$Estatura, prob=seq(0, 1, length = 11))
show(deciles)
##    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
## 1.200 1.210 1.220 1.227 1.246 1.260 1.270 1.280 1.282 1.291 1.300
percentiles <- quantile(tabla1$Estatura, prob=seq(0, 1, length = 101))
show(percentiles)
##     0%     1%     2%     3%     4%     5%     6%     7%     8%     9% 
## 1.2000 1.2029 1.2058 1.2087 1.2100 1.2100 1.2100 1.2100 1.2100 1.2100 
##    10%    11%    12%    13%    14%    15%    16%    17%    18%    19% 
## 1.2100 1.2100 1.2100 1.2100 1.2106 1.2135 1.2164 1.2193 1.2200 1.2200 
##    20%    21%    22%    23%    24%    25%    26%    27%    28%    29% 
## 1.2200 1.2200 1.2200 1.2200 1.2200 1.2200 1.2200 1.2200 1.2212 1.2241 
##    30%    31%    32%    33%    34%    35%    36%    37%    38%    39% 
## 1.2270 1.2299 1.2300 1.2300 1.2300 1.2315 1.2344 1.2373 1.2402 1.2431 
##    40%    41%    42%    43%    44%    45%    46%    47%    48%    49% 
## 1.2460 1.2489 1.2500 1.2500 1.2500 1.2505 1.2534 1.2563 1.2592 1.2600 
##    50%    51%    52%    53%    54%    55%    56%    57%    58%    59% 
## 1.2600 1.2600 1.2600 1.2600 1.2600 1.2600 1.2624 1.2653 1.2682 1.2700 
##    60%    61%    62%    63%    64%    65%    66%    67%    68%    69% 
## 1.2700 1.2700 1.2700 1.2700 1.2700 1.2700 1.2714 1.2743 1.2772 1.2800 
##    70%    71%    72%    73%    74%    75%    76%    77%    78%    79% 
## 1.2800 1.2800 1.2800 1.2800 1.2800 1.2800 1.2800 1.2800 1.2800 1.2800 
##    80%    81%    82%    83%    84%    85%    86%    87%    88%    89% 
## 1.2820 1.2849 1.2878 1.2900 1.2900 1.2900 1.2900 1.2900 1.2900 1.2900 
##    90%    91%    92%    93%    94%    95%    96%    97%    98%    99% 
## 1.2910 1.2939 1.2968 1.2997 1.3000 1.3000 1.3000 1.3000 1.3000 1.3000 
##   100% 
## 1.3000

Podemos ver que la media aritmética y geométrica aparecen con valores muy similares, conociendo la AM-GM inequality corroboramos que la media geométrica es menor que la aritmética, pero por poco. Dado que solo serían iguales si todos los valores fueran idénticos podemos esperar que los datos sigan una distribución leptocúrtica, lo que veremos en el análisis de las medidas de dispersión. No obstante ya mirando los deciles podemos ver que hay un 40% de los datos entre 1,246 y 1,282, luego esta suposición cobra fuerza. Para asegurar nuestra palabras miramos ahora las medidas de dispersión.

rango <- max(tabla1$Estatura) - min(tabla1$Estatura)
cat("El rango del primer set de datos es:", rango)
## El rango del primer set de datos es: 0.1
varianza <- var(tabla1$Estatura)
cat("La varianza del primer set de datos es:", varianza)
## La varianza del primer set de datos es: 0.001050575
std <- var(tabla1$Estatura)^(1/2)
cat("La desviación típica del primer set de datos es:", std)
## La desviación típica del primer set de datos es: 0.03241257
varcoef <- std/media_arit
cat("El coeficiente de variación de Pearson del primer set de datos es:", varcoef)
## El coeficiente de variación de Pearson del primer set de datos es: 0.02586109
library(moments)
asimetria <- skewness(tabla1$Estatura)
cat("El coeficiente de asimetría de los datos es", asimetria)
## El coeficiente de asimetría de los datos es -0.1130675
curtosis <- kurtosis(tabla1$Estatura)
cat("La curtosis de los datos es", curtosis)
## La curtosis de los datos es 1.605299

Todo el rango de nuestras medidas es 0,1 y vemos como la desviación estándar de toda la muestra es apenas 0,03, un 30% del rango, y su coeficiente de variación de Pearson 0,0258. Podemos ver además que los datos son ligeramente asímetricos hacia valores menores que la media y que como habíamos aventurado la distribución es leptocúrtica. Con toda esta información ya sólo nos queda mirar el diagrama de caja y bigotes:

Ejercicio 1.b

Dados los siguientes datos dibujar la tabla de correspondencias y las medidas de centralización

Lo primero cargamos y mostramos los datos:

library(readr)
tabla2 <- read_csv("tabla2.txt")
tabla2$Alumnos <- c(1:nrow(tabla1))
tabla2 <- tabla2[,c(3,1,2)]
Segundo set de datos
Alumnos Estatura Peso
1 1,23 32
2 1,28 33
3 1,27 31
4 1,21 34
5 1,22 32
6 1,29 31
7 1,30 34
8 1,24 32
9 1,27 32
10 1,29 35
11 1,25 31
12 1,28 35
13 1,27 34
14 1,21 33
15 1,22 33
16 1,29 31
17 1,30 35
18 1,24 32
19 1,27 31
20 1,29 33
21 1,25 33
22 1,28 32
23 1,27 34
24 1,21 34
25 1,22 35
26 1,29 31
27 1,30 34
28 1,24 33
29 1,27 35
30 1,29 34

Ahora lo que hago es elegir los intervalos en los que categorizaré cada variable. Para ello voy a usar la media (\(\mu\)) y la desviación típica (\(\sigma\)) y definiré las categorías: muy bajo/delgado: \([min, \mu - \sigma]\), bajo/delgado: \((\mu - \sigma, \mu]\), alto/gordo: \((\mu, \mu + \sigma]\) y muy alto/gordo: \((\mu + \sigma, max]\). Convertimos los datos a estos factores y realizamos la tabla de correspondencias:

tabla2.cat <- tabla2
colnames(tabla2.cat) <- c("Alumno", "EstaturaCat", "PesoCat")
factor_4 <- function (values, desired.labels){
   lim0 <- min(values)
   lim1 <- mean(values) - sqrt(var(values))
   lim2 <- mean(values)
   lim3 <- mean(values) + sqrt(var(values))
   lim4 <- max(values)   
   category <- cut(values, c(lim0,lim1,lim2,lim3,lim4), include.lowest=c(T,F,F,F), labels=desired.labels)
}
tabla2.cat$PesoCat <- factor_4(tabla2$Peso, c("muy delgado", "delgado", "gordo", "muy gordo"))
tabla2.cat$EstaturaCat <- factor_4(tabla2$Estatura, c("muy bajo", "bajo", "alto", "muy alto"))
tabla2.contin <- table(tabla2.cat[,c(2,3)])
Correspondencias segundo set de datos
muy delgado delgado gordo muy gordo
muy bajo 0 2 4 1
bajo 1 2 2 0
alto 5 2 5 3
muy alto 0 0 2 1

Ejercicio 2.a

Con los datos del primer apartado del primer ejercicio, dividir la muestra en 2 submuestras (cada una de 15 alumnos) y obtener un intervalo de confianza para la diferencia de medias de las dos submuestras

Suponemos que ambas muestras son independientes y que las estaturas se encuentran distribuidas normalmente, con varianzas desconocidas pero iguales. Suponemos además que queremos determinar la diferencia de medias con un intervalo de confianza del 95%. Para ello aplicaremos la siguiente fórmula, en la que ya se incluye la confianza y el tamaño de las muestras: \[ IC_{0.05}(\mu_1 - \mu_2) = \left( \bar{x} - \bar{y} - t_{28, 0.025}\sqrt{\frac{\sigma_{x}^{2}\ +\ \sigma_{y}^{2}}{14}}\ ,\ \bar{x} - \bar{y} + t_{28, 0.025}\sqrt{\frac{\sigma_{x}^{2}\ +\ \sigma_{y}^{2}}{14}} \right) \] Así vemos que solo necesitamos las medias y varianzas de las muestras, y el valor de la t de Student para nuestro caso. Lo calculamos con R.

set1<- tabla1[seq(1,15),2]$Estatura
set2<- tabla1[seq(16,30),2]$Estatura
student.test <- t.test(set1, set2, aternative="two.sided", var.equal=T, conf.level=0.95)
t.value<- 2.0484
termino1<- mean(set1) - mean(set2)
termino2<- sqrt((var(set1) + var(set2))/15)
lim.inf<- termino1 - t.value*termino2
lim.sup<- termino1 + t.value*termino2
cat("El intervalo de confianza al 95% es: [",lim.inf, ",", lim.sup, "]")
## El intervalo de confianza al 95% es: [ -0.01079349 , 0.03746016 ]
student.test$conf.int
## [1] -0.01079358  0.03746024
## attr(,"conf.level")
## [1] 0.95
# OFFTOPIC
#Diferencias entre la varianza poblacional y muestral
var.pobla<- sum((set1 - mean(set1))^2)/length(set1)
var.muest<- sum((set1 - mean(set1))^2)/(length(set1) -1)
var.R<- var(set1)
cat("Varianzas: poblacional:", var.pobla, "; muestral:", var.muest, "; con R:", var.R)
## Varianzas: poblacional: 0.0009333333 ; muestral: 0.001 ; con R: 0.001

Hay que ver que los resultados obtenidos con la función t.test solo difieren a partir de la séptima cifra decimal de los resultados obtenidos analíticamente, esto seguramente se deba a que el valor usado en la parte analitica procede de una tabla, y será menos preciso. Por otro lado conviene aclarar que en el segundo término aparece un factor \(14/15\) para lo que hay dentro de la raíz. Esto se debe a que R calcula la varianza muestral \(s^2\ =\ \tfrac{n}{n-1}\ \sigma^2\), si usamos esta varianza en el denominador debe aparecer el tamaño total de ambas muestras. Si usáramos la varianza poblacional en el denominador tendríamos la suma de los tamaños de ambas muestras restando 1 (pérdida de un grado de libertad) a cada una de ellas. Creo que en los apuntes aparece la varianza muestral y el denominador asociado a la varianza poblacional (también creo que esto se comentó en clase), hecho de esa forma los resultados diferirían de los del t.test. Mas allá de estos rigores asociados al pequeño tamaño de las muestras, hemos obtenido un intervalo de confianza al 95% para la diferencia entre las medias de ambas muestras valiéndonos de funciones predefinidas en R y usando las expresiones estudiadas en clase.

Ejercicio 2.b

La hípotesis nula es: \(H_0\ :\ \mu_x \ =\ \mu_y\). Lo que hacemos es rechazar la hipótesis nula de que ambas medias coinciden si: \[|t_{28, 0.025}|\ \lt\ \frac{\bar{x}\ -\ \bar{y}}{\sqrt{\frac{\sigma_{x}^{2}\ +\ \sigma_{y}^{2}}{n - 1}}} \rightarrow P\left(|t_{28, 0.025}|\ \gt\ \frac{\bar{x}\ -\ \bar{y}}{\sqrt{\frac{\sigma_{x}^{2}\ +\ \sigma_{y}^{2}}{n - 1}}}\right)\ \lt\ \alpha\] Lo hago en R:

# Ya extraemos t(28,0.025) de la distribución incluida en R
t.value<- qt(0.975,28)
termino1<- abs(mean(set1) - mean(set2))
termino2<- sqrt((var(set1) + var(set2))/15)
comparator<- termino1/termino2
cat("Rechazamos la hipótesis nula si:", t.value, "<", comparator, "lo que es:", t.value < comparator)
## Rechazamos la hipótesis nula si: 2.048407 < 1.132018 lo que es: FALSE
# Comprobamos que el p-valor obtenido con el test coincide con el obtenido analíticamente
p.anal.value<- 2*(1 - pt(comparator, 28))
cat("p- valor procedente del test:", student.test$p.value, "y procedente del análisis: ", p.anal.value)
## p- valor procedente del test: 0.2672288 y procedente del análisis:  0.2672288
cat("Rechazamos la hipótesis nula si:", p.anal.value, "<", 0.05, "lo que es:", p.anal.value < 0.05)
## Rechazamos la hipótesis nula si: 0.2672288 < 0.05 lo que es: FALSE

Cabe señalar que de nuevo llegamos a las mismas conclusiones analíticamente y con los test predefinidos de R. También debemos tener ahora en cuenta la relacion entre varianza muestral y poblacional para conservar la coherencia de los resultados. Como hemos visto: No rechazamos la hipótesis nula, toda la comparativa es consistente, ya que ambas muestras procedían de la misma muestra en origen.

Y hasta aquí la entrega de estadística, ha sido un placer y muchas gracias por todo.

\[~\] \[~\] \[~\]