1 Ejercicios y problemas estadística no paramétrica Parte II

Lea cuidadosamente y responde las siguientes ejercicios

1.1 Prueba McNemar

1.1.1 Ejercicio 1 [Escrito por Ayrton Pablo Almada]

Como se mencionó en clase la tabla de McNemar las variables \(B\) y \(C\), hacen referencia al número de elementos que cambiaron de un estado al otro, el valor que se esperaría que cambiaran de un estado a otro es \(\frac{B+C}{2}\) y con una varianza \(\frac{B+C}{4}\). Si se supone normalidad, demostrar que:

\[\frac{(B-C)^{2}}{B+C}\sim\chi_{(1)}^{2}\] \(\underline{\text{Solución:}}\)

Suponiendo normalidad, tenemos que la variable \(B\) satisface lo siguiente:

\[B=\sum_{i=1}^{n}\mathbb{I}_{(x_i=1,y_i=0)}\sim \mathcal{N}\left(\frac{B+C}{2},\frac{B+C}{4}\right)\] Por lo que, al ocupar algunas de las propiedade que satisface la distribución normal tendremos lo siguiente:

\[ \begin{aligned} B&\sim \mathcal{N}\left(\frac{B+C}{2},\frac{B+C}{4}\right)\implies \frac{B-\frac{B+C}{2}}{\sqrt{\frac{B+C}{4}}}\sim \mathcal{N}(0,1) \\ \therefore &~\frac{\frac{B-C}{2}}{\frac{\sqrt{B+C}}{2}}\sim \mathcal{N}(0,1) \implies \frac{B-C}{\sqrt{B+C}} \sim \mathcal{N}(0,1) \\ \therefore & ~\left(\frac{B-C}{\sqrt{B+C}}\right)^2\sim \chi_{(1)}^2 \\ \therefore & ~\frac{(B-C)^2}{B+C}\sim \chi_{(1)}^2 \end{aligned} \] \(_\blacksquare\)

1.1.2 Ejercicio 2 [Escrito por Jesus Balam Rodriguez Perez]

Una casa de bolsa forma portafolios de inversión clasificados de alto y bajo riesgo. En 2005, de 100 individuos 70 elegían invertir en portafolios conformado por activos de alto riesgo. Después de la crisis de 2008, en 2010, 25 personas que elegían portafolios riesgosos cambiaron de portafolio eligiendo correr menos riesgos. Por otra parte 10 personas cambiaron de portafolios con poco riesgo a carteras riesgosas. ¿Es significativo el cambio en el numero de personas después de la crisis de 2008?

\(\underline{\text{Solución:}}\)

Se analiza con un nivel de significacia del 5%, entonces tenemos: \[ x_i= \text{1 Si la persona i invierte en acciones de alto riesgo}\\ \quad\quad\text{0 Si la persona i invierte en acciones de bajo riesgo} \]

Despues de la crisis de 2008, sea: \[ y_i= \text{1 Si la persona i invierte en acciones de alto riesgo}\\ \quad\quad\text{0 Si la persona i invierte en acciones de bajo riesgo} \]

m <- matrix(c(A=45, B=25, C=10, D=20), byrow = TRUE, nrow = 2)
colnames(m)<- c("y_i=1","y_i=0")
row.names(m)<- c("x_i=1", "x_i=0")
m
##       y_i=1 y_i=0
## x_i=1    45    25
## x_i=0    10    20

Sea \(m = 10 + 25 = 35\) el número de inversores que cambiaron ya se de inversiones de alto riego a inversiones de bajo riesgo o viseversa. Para verificar que no hay un cambio en los inversores debido a la crisis económica se esperaria que los cambios entre los inversores de alto riesgo a bajo y viceversa fuera simetrica, por lo que B que representa el numero de los inversores que cambiaron de alto riesgo a bajo riesgo tenga una esperanza \(\frac{B+C}{2}=\frac{35}{2}=17.5\) y una varianza \(\frac{B+C}{4}=\frac{35}{4}=8.75\) y si además se agrega el supuesto de que \(B\sim N(17.5, 8.75)\), asi entonces:

\[B\sim N(\frac{B+C}{2}=17.5, \frac{B+C}{4}=8.75)\\ \frac{B-\frac{B+C}{2}}{\sqrt{\frac{B+C}{4}}}\sim N(0,1)\\ (\frac{B-C}{\sqrt{B+C}})^2\sim \chi_{(1)}^2 \] Sea \(T=\frac{(B-C)^2}{B+C}\) con el arreglo de continuidad de Yates se tiene: \(T=\frac{(\mid B-C\mid-1)^2}{B+C}=\frac{(\mid 25-10\mid-1)^2}{25+10}=\frac{(14)^2}{35}=\frac{196}{35}=5.6\)

Tenemos que los cuantiles para la distribución ji-cuadrada estan dados por: \(\chi_{(1)^{2(0.25)}}=3.84\) y \(\chi_{(1)^{2(0.975)}}=5.02\)

Ya que \(T=5.6>5.02=\chi_{(1)^{2(0.975)}}\) entonces se rechaza la prueba, lo que iplica que la crisis de 2008 si tuvo una implicación en las inversiones

\(_\blacksquare\)

1.2 Pruebas de Wilcoxon / Kruskal Wallis / Medidas de correlacion

1.2.1 Ejercicio 1 [Escrito por Hugo Reyna Castañeda]

La oficina de Censo reportó que se espera que los hispanos sobrepasen a los afroamericanos como la minoría más grande en los Estados Unidos para el año 2030. Use dos pruebas diferentes para ver si hay una relacion directa entre el número de Hispanos y el porcentaje de la población del estado para los nueve estados que se presentan en la tabla siguiente:

Estado Hispanos (millones) Porcentaje de la población del estado
California 6.6 23
Texas 4.1 24
New York 2.1 12
Florida 1.5 12
Illinois 0.8 7
Arizona 0.6 18
New Jersey 0.6 8
New Mexico 0.5 35
Colorado 0.4 11

Use el nivel de significancia \(\alpha=0.05\)

\(\underline{\text{Solución:}}\)

Para fines de este problema representaremos a cada estado de la tabla anterior con un número natural de acuerdo al orden en que aparecen, es decir:

Estado \(i\)
California 1
Texas 2
New York 3
Florida 4
Illinois 5
Arizona 6
New Jersey 7
New Mexico 8
Colorado 9

Sea \(X_{i}\) la variable que representa el número de habitantes (en millones) Hispanos en el \(i\)-ésimo estado y \(Y_{i}\) la variable que representa el porcetanje de la población en el \(i\)-esimo estado.

X <-c(6.6,4.1,2.1,1.5,0.8,0.6,0.6,0.5,0.4)
Y <-c(23,24,12,12,7,18,8,35,11)
plot(X,Y, pch="x")

  • Coeficiente de correlación de Pearson

Denotaremos por \(r\) a dicho coeficiente, el cual está definido por: \[ r=\frac{\sum_{i=1}^{9}X_{i}Y_{i}-9\bar{X}\cdot\bar{Y}}{\sqrt{ \left(\sum_{i=1}^{9}X_{i}^{2} - 9 \bar{X}^{2} \right)\left( \sum_{i=1}^{9}Y_{i}^{2} - 9 \bar{Y}^{2} \right) }} \]

Cuando \(r=0\) entonces las muestras \(X\) y \(Y\) son independientes de modo que para nuestro interés deseamos realizar la siguiente prueba: \[ Ho:\,\, r=0 \,\,\,\,\,\,vs \,\,\,\,\,\,\, Ha:\,\, r \neq 0. \]

Calculamos \(r\) en \(\mathtt{R}\):

media_X = mean(X)
media_Y = mean(Y)

media_X
## [1] 1.911111
media_Y
## [1] 16.66667
m_2_X = sum(X^2)
m_2_Y = sum(Y^2)

m_2_X
## [1] 68.8
m_2_X
## [1] 68.8
S_i = sum(X*Y)
S_i
## [1] 336.5

Por lo tanto \[ r=\frac{\sum_{i=1}^{9}X_{i}Y_{i}-9\bar{X}\cdot\bar{Y}}{\sqrt{ \left(\sum_{i=1}^{9}X_{i}^{2} - 9 \bar{X}^{2} \right)\left( \sum_{i=1}^{9}Y_{i}^{2} - 9 \bar{Y}^{2} \right) }}=0.3197604 \]

numerador = S_i - 9*(media_X)*(media_Y)
numerador
## [1] 49.83333
denominador = sqrt((m_2_X - 9*(media_X)^2)*(m_2_Y - 9*(media_Y)^2))
denominador
## [1] 155.8458
r = numerador/denominador
r
## [1] 0.3197604

En \(\mathtt{R}\) realizamos la prueba con la función \(\mathtt{cor.test}\). Entonces:

cor.test(X,Y, method = "pearson",  alternative = "greater", exact = TRUE)
## 
##  Pearson's product-moment correlation
## 
## data:  X and Y
## t = 0.89288, df = 7, p-value = 0.2008
## alternative hypothesis: true correlation is greater than 0
## 95 percent confidence interval:
##  -0.3275921  1.0000000
## sample estimates:
##       cor 
## 0.3197604

Así pues, \(p-value=0.2008 > 0.05 = \alpha\) por lo que afirmamos con una confianza del 95% que no existe una relación directa entre el tamaño de la población de habitantes Hispanos y el porcentaje de población de los nueve estados.

\(_\blacksquare\)

  • \(\rho\) de Spearman

Deseamos realizar la prueba: \[ \begin{aligned} &Ho:\,\,\mbox{Las}\,X_{i}'s\,\,\mbox{y}\,\,Y_{i}'s\,\,\mbox{son mutuamente independientes}\\ &vs\\ &Ha:\,\,\mbox{Existe una tendencia entre los valores de las}\,\,X_{i}'s\,\,\mbox{y}\,\,Y_{i}'s \end{aligned} \] A lo que nos referimos en la hipótesis alternativa es que existe una tendencia para que los valores más grandes de \(X\) estén “emparejados” con los valores más grandes de \(Y\) y los valores chicos de \(X\) con los valores chicos de \(Y\) ó existe una tendencia para que los valores \(más\) grandes de \(X\) estén “emparejados” con los valores más chicos de \(Y\) y los valores chicos de \(X\) con los valores grandes de \(Y\).

Obtengamos los rangos de cada muestra:

Observación \(X_{i}\) Rango(\(X_{i}\))
1 6.6 9
2 4.1 8
3 2.1 7
4 1.5 6
5 0.8 5
6 0.6 3.5
7 0.6 3.5
8 0.5 2
9 0.4 1
Observación \(Y_{i}\) Rango(\(Y_{i}\))
1 23 7
2 24 8
3 12 4.5
4 12 4. 5
5 7 1
6 18 6
7 8 2
8 35 9
9 11 3 
R_X <-c(9,8,7,6,5,3.5,3.5,2,1)
R_Y <-c(7,8,4.5,4.5,1,6,2,9,3)

El coeficiente de Spearman está dado por: \[ \rho = \frac{\sum_{i=1}^{9}\left(R(X_i)-\frac{9+1}{2}\right)\left(R(Y_i)-\frac{9+1}{2}\right)}{\frac{9(9^2-1)}{12}}= 0.2416667 \]

rho = (sum((R_X- 5)*(R_Y - 5)))/((9*(9^2-1))/12)
rho
## [1] 0.2416667

En \(\mathtt{R}\) realizamos la prueba a través del comando:

cor.test(X,Y, method = "spearman", alternative = "greater", exact = TRUE)
## Warning in cor.test.default(X, Y, method = "spearman", alternative =
## "greater", : Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  X and Y
## S = 90.756, p-value = 0.2637
## alternative hypothesis: true rho is greater than 0
## sample estimates:
##       rho 
## 0.2436975

Dado que \(p-value = 0.2637 > 0.05 = \alpha\), asumimos con un nivel de significancia del 5% que el tamaño de la población de habitantes Hispanos y el porcentaje de población de los nueve estados guardan una relación directa. Y para tener una mejor apreciación de que las muestras \(X\) y \(Y\) son mutuamente independientes graficamos los rangos establecidos.

plot(R_X,R_Y,pch="x")

\(_\blacksquare\)

  • \(\tau\) de Kendall

Por último, recurrimos a la prueba de Kendall dada de manera similar que las anterior y la cual podemos realizar en \(\mathtt{R}\) a través del siguiente comando:

cor.test(X,Y, method="kendall",alternative="greater",exact = NULL)
## Warning in cor.test.default(X, Y, method = "kendall", alternative =
## "greater", : Cannot compute exact p-value with ties
## 
##  Kendall's rank correlation tau
## 
## data:  X and Y
## z = 0.63236, p-value = 0.2636
## alternative hypothesis: true tau is greater than 0
## sample estimates:
##       tau 
## 0.1714286

Así pues, \(p-value=0.2636 > 0.05 = \alpha\) por lo que no rechazamos \(Ho\) con un nivel de significancia del 5%, es decir, asumimos que las muestras \(X\) y \(Y\) no guardan una relación directa. \(\square\)

En todas las pruebas no rechazamos \(Ho\) con un nivel de significancia del \(5\%\) por lo que asumimos hoy que la población de Hispanos no es significativa en los estados presentados de EUA.

\(_\blacksquare\)

1.2.2 Ejercicio 2 [Escrito por Hugo Reyna Castañeda]

Un psicologo esta investigando el impacto que el divorcio de los padres tiene sobre el aprovechamiento academico de los ninios. El psicologo cuenta con las calificaciones de un grupo de ninos de escuela primaria cuyos padres tuvieron un divorcio durante el anio anterior, y las calificaciones para un grupo de ninos similares cuyos padres no se divorciaron.

\(\underline{\text{Solución:}}\)

no_divorciados divorciados
80 60
72 70
99 88
82 75
62 42
50 30
85 50

Se puede decir que: ¿hay diferencia en el aprovechamiento academico de los niños? Use \(\alpha=0.05\).

Antes de empezar con el análisis, exploremos los datos mediante un boxplot:

library(ggplot2)
library(ggthemes)
ggplot(data_cal, aes(x=Divorciado,Calificaciones))+
  geom_boxplot(fill="turquoise",alpha=.8)+theme_tufte()

A priori, es notorio que existe una diferencia en las calificaciones de estos estudiantes. Comprobemos esta aseveracion. \(\\\) Procedamos por la prueba de Wilcoxon para verificar si existen diferencias significativas en las calificaciones de estos alumnos. Esto es, queremos contrastar

\[H_0:E[X]=E[Y]~~~vs~~~H_a:E[X]\neq E[Y]\] Primeramente, asignemos rangos de manera ascendente a la muestra combinada. En caso de empate, asignaremos el punto medio entre ambos valores.

muestra_combinada<-c(no_divorciados,divorciados)
n=length(muestra_combinada)

#Ordenando la muestra de manera ascendente

rango=rank(muestra_combinada)
muestra_comb_ord<-muestra_combinada[order(muestra_combinada)]
rango_2=rank(muestra_comb_ord)

b=data.frame("Muestra combinada"=muestra_comb_ord,"Rango"=rango_2)
kable(b)
Muestra.combinada Rango
30 1.0
42 2.0
50 3.5
50 3.5
60 5.0
62 6.0
70 7.0
72 8.0
75 9.0
80 10.0
82 11.0
85 12.0
88 13.0
99 14.0

Ahora, sumemos los rangos de la primera muestra, es decir, de los alumnos con padres no divorciados y calculemos la estadistica de prueba \(T=S-\frac{n(n+1)}{2}\).

## [1] "La suma de rangos es:  64.5"
## [1] "La estadistica T es:  36.5"

Sin embargo, debido a la presencia de empates, resulta conveniente realizar el calculo de la estadistica \(T_1=\frac{T-n* \frac{N+1}{2}}{\sqrt{\frac{mn}{N(N-1)} \sum_1^N R^2_i- \frac{mn(N+1)^2}{4(N-1)}}}\). Note que \(m=n=7 \Rightarrow N=m+n=14\). De esta forma tenemos

n=length(no_divorciados)
m=length(divorciados)
N=m+n
R=sum(rango^2)
T1 =(Tw - n*(N +1)/2)/sqrt(m*n*R/( N*(N -1) ) - m*n*( N+1)^2/(4*(N -1) ) )

paste("La estadistica T1 es: ",T1)
## [1] "La estadistica T1 es:  -2.04665531689017"

Notemos que al usar un nivel de signficancia del \(\alpha=5\%\), los cuantiles \(W^{\alpha/2}=W^{0.025}\) y \(W^{1-\alpha/2}=W^{0.975}\) son los cuantiles de una normal estandar, por lo tanto, tomando en cuenta la simetria de la distribucion normal, \(W^{0.975}=1.95\) y \(W^{0.025}=-1.95\) De los resultados anteriores obtenemos que \(T1=-2.046<W^{0.025}=-1.95\), por lo que, con un nivel de confianza del \(95\%\) se rechaza la hipotesis nula de que las poblaciones son iguales, i.e. existen diferencias en las calificaciones de alumnos con padres divorciados y padres no divorciados.

Al realizar la prueba mediante la funcion wilcox_test, se obtienen los siguientes resultados.

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  data_cal$Calificaciones by data_cal$Divorciado
## V = 28, p-value = 0.02154
## alternative hypothesis: true location shift is not equal to 0

Notemos que \(p-value<0.05\), por lo que se comprueban nuestros resultados que nos permiten rechazar la hipotesis nula de igualdad de medias.

\(_\blacksquare\)

1.2.3 Ejercicio 3 [Escrito por Hugo Reyna Castañeda]

La tabla que se proporciona a continuación da el número de premios de postgraduados en ciencia médica y la razón de muerte por millón de tuberculosis para los años 1959-69 (fuente: Annual Abstract of Statistics 1970).

Año Número de premios Tasa de muerte por tuberculosis
1959 277 83
1960 318 74
1961 382 71
1962 441 65
1963 486 62
1964 597 52
1965 750 47
1966 738 48
1967 849 42
1968 932 43
1969 976 38

Demuestre que estos datos muestran una fuerte evidencia de correlación negativa entre el número de premios y la tasa de muerte por tuberculosis. Explique este “extraño” resultado. Use \(\alpha=0.05\).

\(\underline{\text{Solución:}}\)

Representaremos a cada año con un respectivo número natural de acuerdo al orden en que aparecen en la tabla anterior, es decir:

Año \(i\)
1959 1
1960 2
1961 3
1962 4
1963 5
1964 6
1965 7
1966 8
1967 9
1968 10
1969 11

Sea \(X_i\) la variable que representa el número de premios en el \(i\)-ésimo año y \(Y_{i}\) la variable que representa la tasa de mortalidad de tuberculósis en el \(i\)-ésimo año.

X_3 <- c(277,318,382,441,486,597,750,738,849,932,976)
Y_3 <- c(83,74,71,65,62,52,47,48,42,43,38)

Haremos uso del coeficiente de Spearman entonces, deseamos realizar la prueba:

\[ \begin{aligned} \begin{aligned} &Ho:\,\,\mbox{Las}\,X_{i}'s\,\,\mbox{y}\,\,Y_{i}'s\,\,\mbox{son mutuamente independientes}\\ &vs\\ &Ha:\,\,\mbox{Existe una tendencia negativa entre los valores de las}\,\,X_{i}'s\,\,\mbox{y}\,\,Y_{i}'s \,\,(\rho < 0) \end{aligned} \end{aligned} \]

A lo que nos referimos en la hipótesis alternativa es que existe una tendencia para que los valores más grandes de \(X\) estén “emparejados” con los valores más chicos de \(Y\) y que los valores más chicos de \(X\) estén “emparejados” con los valores más grandes de \(Y\).

Obtengamos los rangos de cada muestra:

Observación \(X_{i}\) Rango(\(X_{i}\))
1 277 1
2 318 2
3 382 3
4 441 4
5 486 5
6 597 6
7 750 8
8 738 7
9 849 9
10 932 10
11 976 11
Observación \(Y_{i}\) Rango(\(Y_{i}\))
1 83 11
2 74 10
3 71 9
4 65 8
5 62 7
6 52 6
7 47 4
8 48 5
9 42 2
10 43 3
11 38 1 
R_X_3 <-c(1,2,3,4,5,6,8,7,9,10,11)
R_Y_3 <-c(11,10,9,8,7,6,4,5,2,3,1)

T_3 = sum((R_X_3 - R_Y_3)^2)
T_3
## [1] 438

Dado que nuestra muestra no presenta empates entonces la estadística de prueba \(T\) está dada por: \[ T=\sum_{i=1}^{9}(R(X_i)-R(Y_i))^2=438 \]

Por tanto, el coeficiente de Spearman es: \[ \rho = 1 -\frac{6(483)}{11(11^2-1)}=-0.9909091 \]

rho_3 = 1-((6*T_3)/(11*(11^2-1)))
rho_3
## [1] -0.9909091

Obtenemos de la tabla B.5 de las notas que con un nivel de confianza del 95%, \(\omega^{0.05}= 0.536\). Así pues, \(\rho < \omega^{0.05}\) por lo que rechazamos \(Ho\) con un nivel de significancia del 5%, es decir, existe una relación negativa entre el número de premios y la tasa de muerte por tuberculosis.

En \(\mathtt{R}\) realizamos la prueba a través del comando:

cor.test(X_3 ,Y_3 , method = "spearman", alternative = "greater", exact = TRUE)
## 
##  Spearman's rank correlation rho
## 
## data:  X_3 and Y_3
## S = 438, p-value = 1
## alternative hypothesis: true rho is greater than 0
## sample estimates:
##        rho 
## -0.9909091

Podemos inferir desde las tablas de rango de las muestras que los valores grandes de \(X\) están emparejados con los valores chicos de \(Y\) y que los valores más chicos de \(X\) están emparejados con los valores más grandes de \(Y\) pues no existen empates. Sin embargo, esta relación que guardan en realidad es muy poca y es meramente por los datos numéricos que arrojan sus muestras. Podemos conocer más ejemplos en el libro de Virgen, T. “Spurious correlations” que muestra las gráficas entre muestras que numéricamente guardan una relación pero en la vida real no son casuales.

\(_\blacksquare\)

1.2.4 Ejercicio 4 [Escrito por Ayrton Pablo Almada]

El personal de un hospital mental desea saber qué clase de tratamiento es más efectivo para un tipo particular de desorden mental. Una batería de pruebas administrada a todos los pacientes delineó a un grupo de 40 pacientes quienes fueron considerados de diagnóstico similar y también personalidad, inteligencia y factores fisiológicos y proyectivos. Esta gente fue dividida en cuatro diferentes grupos de 10 cada uno para tratamiento. Durante seis meses los grupos respectivos recibieron (1) electroshock, (2) psicoterapia, (3) electroshock más psicoterapia, y (4) ningún tipo de tratamiento. Al final de este período la batería de pruebas fue repetida en cada paciente. El único tipo de medida posible para estas pruebas es un ordenamiento (ranking) de los 40 pacientes de acuerdo a su grado relativo de mejoría al final del período de tratamiento; rango 1 indica el nivel más alto de mejoría, rango 2 el segundo mejor, y así sucesivamente. De acuerdo con estos datos, ¿existe diferencia en efectividad de los tipos de tratamiento? Use \(\alpha = 0.05\).

\[\text{Grupos}\]

1 2 3 4
19 14 12 38
22 21 1 39
25 2 5 40
24 6 8 30
29 10 4 31
26 16 13 32
37 17 9 33
23 11 15 36
27 18 3 34
28 7 20 35

\(\underline{\text{Solución:}}\)

Empezamos definiendo la prueba de hipótesis.

Sean \(X_i\) la variable aleatoria del ranking del tratamiento de un individuo del grupo \(i\), la prueba de hipótesis se define como la igualdad en los valores esperados de cada una de las 4 muestras aleatorias ó grupos, en otras palabras:

\[H_0:\mathbb{E}[X_1]=\mathbb{E}[X_2]=\mathbb{E}[X_3]=\mathbb{E}[X_4] ~~~\text{vs.}~~~H_1:\exists_{j\ne i\in\{1,2,3,4\}}|~\mathbb{E}[X_i]\ne \mathbb{E}[X_j].\]

Nótese que el número total de observaciones a analizar es \(N=\sum_{i=1}^{k}n_i=40\), (4 muestras de 10 elementos cada una).

Primero obtendremos el rango de todas y cada una de las observaciones de la tabla anterior, los rangos marginales \(R_i\),por lo que ocuparemos el siguiente código:

M<-matrix(c(19,14,12,38,22,21,1,39,25,2,5,40,24,6,8,30,29,10,4,31,26,16,13,32,
37,17,9,33,23,11,15,36,27,18,3,34,28,7,20,35),nrow=10,ncol=4,byrow=T) 
#La anterior tabla la convertimos en matriz
N<-dim(M)[1]*dim(M)[2] #Obtenemos el número total de observaciones
#R<-matrix(rank(M),nrow=10,ncol=4,byrow=T) 
R<-M #Obtenemos el rango de todas las observaciones
RM<-colSums(R) #Obtenemos los rangos marginales

Seguido, calculamos la estadística de Kruskal-Wallis \(T\) (nótese que en este caso no hay empates en los rangos):

\[ \begin{aligned} T=& \frac{1}{S^2}\left(\sum_{i=1}^{k}\frac{R_i^2}{n_i}-\frac{N(N+1)^2}{4}\right)\\ :S=& \frac{N(N+1)}{12}:\text{no se presentan empates en los rangos.} \end{aligned} \]

S1<-sum((RM^2)/10) 
#Obtenemos la suma de rangos marginales al cuadrado sobre sus tamaños de muestra
S2<-sum(R^2) 
#Obtenemos la suma de todos los rangos al cuadrado
Scuad<-N*(N+1)/12
T<-(1/(Scuad))*(S1-(40*(41)^2)/(4))
paste("T = ",T)
## [1] "T =  31.8936585365854"

Dado que \(T\sim\chi^{2}_{k-1}\), le regla de desición de la prueba es rechazar \(H_0\) si \(T>\chi^{2(1-\alpha)}_{k-1}\), por ende:

cuantil<-qchisq(1-0.05,3)
paste("Cuantil = ",cuantil)
## [1] "Cuantil =  7.81472790325118"

Como \(T=31.89366>7.814728=\chi^{2(0.95)}_{3}\), rechazamos la hipótesis nula y concluímos que sí existe diferencia en efectividad de los tipos de tratamiento.

Ahora, realizando la prueba de forma automática:

x<-M[,1]
y<-M[,2]
z<-M[,3]
w<-M[,4]
l<-list(x,y,z,w)
kruskal.test(l)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  l
## Kruskal-Wallis chi-squared = 31.894, df = 3, p-value = 5.511e-07

Tendremos un \(\text{p-value}=5.511e-07<\alpha=0.05\), lo cual sólo termina por confirmar nuestra conclusión anterior.

\(_\blacksquare\)

1.2.5 Ejercicio 5 [Escrito por Jesús Balam Rodríguez Pérez]

En el archivo pregunta1.r se encuentran 1000 pares de datos \((X_i,Y_i)\). Calcule el coeficiente de correlación de Pearson, la \(\rho\) de Spearman y pruebe \(H_0\) : Las \(X_i\)’s y las \(Y_i\)’s son mutuamente independientes.

\(\underline{\text{Solución:}}\)

  • Coeficiente de correlación de Pearson
## Parsed with column specification:
## cols(
##   Xi = col_double(),
##   Yi = col_double()
## )
## Warning: 1000 parsing failures.
## row col  expected    actual          file
##   1  -- 2 columns 3 columns 'pregunta1.r'
##   2  -- 2 columns 3 columns 'pregunta1.r'
##   3  -- 2 columns 3 columns 'pregunta1.r'
##   4  -- 2 columns 3 columns 'pregunta1.r'
##   5  -- 2 columns 3 columns 'pregunta1.r'
## ... ... ......... ......... .............
## See problems(...) for more details.
## 
##  Pearson's product-moment correlation
## 
## data:  X1 and Y1
## t = -0.82627, df = 998, p-value = 0.7956
## alternative hypothesis: true correlation is greater than 0
## 95 percent confidence interval:
##  -0.07808601  1.00000000
## sample estimates:
##        cor 
## -0.0261463

El coeficiente de correlación de Pearson es -0.0261463

  • \(\rho\) de Spearman
cor.test(X1, Y1, method="spearman",alternative="greater",exact = TRUE)
## 
##  Spearman's rank correlation rho
## 
## data:  X1 and Y1
## S = 170732778, p-value = 0.7796
## alternative hypothesis: true rho is greater than 0
## sample estimates:
##         rho 
## -0.02439769

\(H_0\) : Las \(X_i\)’s y las \(Y_i\)’s son mutuamente independientes con un grado de significancia del 95%, \(w_{0.025}=0.255\) y \(w_{0.975}=0.745\) ya que \(\rho\) esta entre estos dos valores se acepta la hipotesis nula es decir las \(X_i's\) y \(Y_i's\) son mutuamente independientes.

\(_\blacksquare\)

1.2.6 Ejercicio 6 [Escrito por Carlos Alberto Gómez Correa]

En el archivo pregunta2.r se encuentran 1000 pares de datos \((X_i, Y_i)\). Calcule el numero de parejas concordantes, el numero de parejas discordantes, el numero de empates, la estadıstica \(T = N_c − N_d\) y la \(\tau\) de Kendall y pruebe la hipotesis nula de independencia de las \(X_i's\) y las \(Y_i's\)

\(\underline{\text{Solución:}}\)

Antes de realizar la prueba, efectuemos una grafica:

## Parsed with column specification:
## cols(
##   Xi = col_double(),
##   Yi = col_double()
## )
## Warning: 500 parsing failures.
## row col  expected    actual          file
##   1  -- 2 columns 3 columns 'pregunta2.r'
##   2  -- 2 columns 3 columns 'pregunta2.r'
##   3  -- 2 columns 3 columns 'pregunta2.r'
##   4  -- 2 columns 3 columns 'pregunta2.r'
##   5  -- 2 columns 3 columns 'pregunta2.r'
## ... ... ......... ......... .............
## See problems(...) for more details.

Podemos notar que a prior hay una correlacion positiva entre los datos. \(\\\) Calculemos el numero de parejas concordantes, discordantes y empates mediante la funcion Pares() definida como sigue:

#Funcion que calcula las parejas concordantes, discordantes y empates
Pares <- function(x,y){
  n <- length(x)
  ix <- order(x)
  x <- x[ix]
  y <- y[ix]
  
  #Concordantes
  Nc <- sum(sapply(1:(n-1),function(i) sum(x[i]<x[(i+1):n] & y[i]<y[(i+1):n])))
  
  #Discordantes
  Nd <- sum(sapply(1:(n-1),function(i) sum(x[i]<x[(i+1):n] & y[i]>y[(i+1):n])))
  
  #Empates por la derecha
  Nr<- sum(sapply(1:(n-1),function(i) sum(y[i]==y[(i+1):n])))
 #Empares por la izquiera
   Nl<- sum(sapply(1:(n-1),function(i) sum(x[i]==x[(i+1):n])))
 #Empates en ambas entradas
    Nb<- sum(sapply(1:(n-1),function(i) sum(x[i]==x[(i+1):n] & y[i]==y[(i+1):n])))
    
    return(list(Nc=Nc,Nd=Nd,Nl=Nl,Nr=Nr,Nb=Nb))
}

Pares(x_i,y_i)
## $Nc
## [1] 124748
## 
## $Nd
## [1] 0
## 
## $Nl
## [1] 0
## 
## $Nr
## [1] 2
## 
## $Nb
## [1] 0

Notemos que el numero de parejas concordantes es \(Nc=101168\); el numero de parejas discordantes \(N_d=23578\); el numero de empates en la primera entrada es \(Nl=2\); el numero de empates en la segunda entrada es \(Nr=2\), y el numero de empates en ambas entradas es \(N_b=0\).

Procedamos ahora al calculo de la estadistica \(T= N_c-N_d\). Por lo anterior, tenemos \[T=N_c-N_d=101168-23578=77590\] \[\therefore~ T=77590\]

Ahora, calculemos la \(\tau\) de Kendall dada por \(\tau=\frac{N_c-N_d}{n(n-1)/2}\)

n<-length(x_i)
Nc<-as.numeric(Pares(x_i,y_i)[1])
Nd<-as.numeric(Pares(x_i,y_i)[2])
Te<-Nc-Nd
tau_1=(Nc-Nd)

tau_1=Te/(n*(n-1)/2)
tau_1
## [1] 0.999984

Por lo tanto \(\tau=0.62196\).

\(\\\)

Calculemos ahora la version dos de la tau de Kendall para presencia de empates, a saber \(\tau=\frac{N_c-N_d}{N_c+N_d}\)

## [1] 1

Por lo que esta version da como resultado \(\tau=0.62198\). Notemos que ambas versiones de esta medida son casi identicas.

\(\\\)

Probemos la hipotesis de independencia de \(X_i\) y \(Y_i\), de una cola, es decir, buscamos probar \[H_0:Las~X_i's~y~las~Y_i's~son~mutuamente~independientes\] contra la hipotesis alternativa \(H_a:\) Existe una tendencia para que los valores mas grandes de X esten emparejados con los valores mas chicos de Y y que los valores mas chicos de X esten emparejados con los valores mas grandes de Y (\(\tau < 0\)). Por lo tanto, rechazaremos \(H_0\) a un nivel de significancia usual de \(\alpha=5\%\) si \(\tau>\omega_{1-\alpha/2}\), donde \(\omega_p\) es el p-esimo cuantil de la distribucion de \(\tau\), que es aproximadamente:

\[\omega_p=z_p\frac{\sqrt {2(2n+5)}}{3 \sqrt{n(n-1)}}\] donde \(z_p\) es el p-esimo cuantil de una distribucion normal estandar.

\(\\\)

Notemos que para nuestro nivel de significancia \(\alpha=0.05\), por la simeria de la distribucion normal \(z_{\alpha/2}=-1.959964\) y \(z_{1-\alpha/2}=-z_{0.025}=1.959964\).

## [1] -0.05863941
## [1] 0.05863941

Por lo tanto \[\omega_{\alpha/2}=z_{\alpha/2}\frac{\sqrt {2(2(500)+5)}}{3 \sqrt{500(500-1)}}=-1.959\frac{\sqrt {2(2(500)+5)}}{3 \sqrt{500(500-1)}}=-0.05863941\]

Y \[\omega_{1-\alpha/2}=-\omega_{\alpha/2}=0.05863941\]

\[\therefore~ \tau=0.612>\omega_{1-\alpha/2}=0.05863941\] Por lo que tenemos evidencia suficiente para rechazar la hipotesis nula de independencia entre las variables, a un nivel de significancia del \(5\%\) \(\\\)

Nuestros resultados se comprueban realizando la prueba mediante la funcion cor.test() como sigue:

## 
##  Kendall's rank correlation tau
## 
## data:  x_i and y_i
## z = 33.423, p-value < 2.2e-16
## alternative hypothesis: true tau is greater than 0
## sample estimates:
##      tau 
## 0.999992

Por lo anterior, \(p-value<0.05=\alpha\). En definitiva, se rechaza la hipotesis nula de independencia entre las variables \(X_i's\) y \(Y_i's\)

\(_\blacksquare\)

1.2.7 Ejercicio 7 [Escrito por Jesús Balam Rodríguez Pérez]

Para las dos muestras que se encuentran en el archivo pregunta3.r pruebe utilizando la prueba de Wilcoxon \(H_0\) : Las medias son iguales.

\(\underline{\text{Solución:}}\)

 X<- c(96.589,98.916,98.89,101.784,98.196,96.747,99.49,100.696,100.219,101.651,
       95.934,96.609,96.585,97.836,99.787,99.607,100.128,103.619,97.402,104.652,
       94.773,98.149,93.087,101.407,94.227,95.529,98.306,99.359,96.456,93.067,
       100.063,97.363,99.172,94.568,94.247,97.981,101.759,98.767,100.46,95.077,
       98.442,97.322,101.76,96.212,92.216,99.805,94.68)

Y<- c(101.293,97.878,99.522,101.097,97.558,101.048,92.328,96.704,98.756,96.915,
      102.158,99.452,98.813,101.063,102.584,102.817,96.901,92.854,104.189,
      100.312,99.581,99.216,103.894,100.83,104.165,97.844,104.793,102.402,
      101.861,101.045,98.69,98.28,102.539,98.939,95.834,103.356,102.969,
      103.858,99.38,98.856,96.54,99.719,96.926,99.665,101.866,99.334,100.548,
      94.712,102.494,97.843,103.684,97.97,101.549,101.026,98.174,95.375,
      101.682,97.606)

n=length(X)
m=length(Y)
N=n+m

a<- c(X,Y)
rango<- rank(a)

S<- sum(rango[1:47])

Te=S-n*(n+1)/2

R= sum(rango^2)

T1=(Te-n*(N+1)/2)/sqrt(n*m*R/(N*(N-1)) - m*n*(N+1)^2/(4*(N-1))) 
T1
## [1] -10.29797

El resultado final es que \(T_1=-10.29797<-1.96=W_{0.025}\) con lo que se rechaza la hipotesis nula es decir concluimos que las medias no son iguales.

\(_\blacksquare\)

1.2.8 Ejercicio 8 [Escrito por Erick Arturo Cortés Galán]

Realice la prueba de Friedman para los datos que se dan en el archivo pregunta4.r

\(\underline{\text{Solución:}}\)

Debemos tener presente que la prueba de Friedman nos es útil para comparar si el promedio de datos entre \(k\) elementos tiene diferencias o no. Por lo que establecemos:

\(H_0\): Los promedios entre los datos de los k elementos es parecida.
\(H_a\): Hay algunas diferencias entre los promedios.

Primeramente creamos una matriz con los datos

A=matrix(c(83,77,88,87,92,96,110,109,124,119,129,125,
           127,123,132,129,124,117,122,123),nrow = 5,ncol=4,byrow=TRUE)
A
##      [,1] [,2] [,3] [,4]
## [1,]   83   77   88   87
## [2,]   92   96  110  109
## [3,]  124  119  129  125
## [4,]  127  123  132  129
## [5,]  124  117  122  123

Podemos observar que tenemos \(k=4\) elementos a comparar y \(n=5\) renglones de datos.
Ahora llevaremos acabo la asignación de rangos por renglón:

#Asignamos rangos
c<-c()
for (k in 1:dim(A)[1]){
  reng<-A[k,]
  b<-sort(reng)
  for (i in 1:dim(A)[2]){
    for(j in 1:dim(A)[2]){
      if(reng[i]==b[j]){
        c=c(c,j)
      }
    }
  }
}
matrang=matrix(c,nrow=dim(A)[1],ncol=dim(A)[2],byrow=TRUE)
matrang
##      [,1] [,2] [,3] [,4]
## [1,]    2    1    4    3
## [2,]    1    2    4    3
## [3,]    2    1    4    3
## [4,]    2    1    4    3
## [5,]    4    1    2    3

También es necesario crear el vector que contiene la suma de los rangos y media los rangos por elemento, es decir, por columna

suma_rangos<-c()
for (i in 1:dim(A)[2]){
  suma_rangos<-c(suma_rangos,sum(matrang[,i]))
  
}
suma_rangos
## [1] 11  6 18 15
media_rangos=suma_rangos/dim(A)[1]
media_rangos
## [1] 2.2 1.2 3.6 3.0

En la prueba de Friedman se comparan los promedios de los rangos con el promedio general. Así que es necesario obtener el promedio general

#La media general R de los rangos es:
media_gral_rangos=(dim(A)[2]+1)/2

Procedemos al cálculo de la suma de los cuadrados de las diferencias de la media de los rangos y la general \[ \begin{aligned} S=& \sum_{j=1}^{k}( \bar{R_j} -\bar{R} )^2 \end{aligned} \]

S=0
for(i in 1:dim(A)[2]){
  S=S+(media_rangos[i]-media_gral_rangos)^2
}
S
## [1] 3.24

Si \(H_0\) es cierta, nuestra estadística debe tomar valores pequeños, pues significa que las medias de rangos tienen valores cercanos a \(\bar{R}\). Caso contrario \(H_a\) es verdadera.

Ahora obtenemos la estadística de Friedman, la cual es un múltiplo de \(S\) de la forma: \[ \begin{aligned} M=& \frac{12n}{k(k+1)}S \end{aligned} \]

##Obtenemos la estadística de Friedman
M=(12*dim(A)[1]/(dim(A)[2]*(dim(A)[2]+1)))*S
M
## [1] 9.72

Se compara este estadístico con los valores en la tabla, pues los valores de k y n, observaciones y elementos, si aparecen en ella. En este caso, \(n=5\) y \(k=4\), por lo que con un nivel de significancia \(\alpha=.01\), nuestra región crítica es: \[ \begin{aligned} M\geq9.96 \end{aligned} \] Por lo tanto, no tenemos evidencia suficiente para rechazar \(H_0\).

Comparamos con la prueba cargada en R:

A
##      [,1] [,2] [,3] [,4]
## [1,]   83   77   88   87
## [2,]   92   96  110  109
## [3,]  124  119  129  125
## [4,]  127  123  132  129
## [5,]  124  117  122  123
friedman.test(A)
## 
##  Friedman rank sum test
## 
## data:  A
## Friedman chi-squared = 9.72, df = 3, p-value = 0.0211

Observamos que nuestros resultados coinciden con la prueba que desarrollamos, no rechazamos nuestra hipótesis nula con un nivel de significancia \(\alpha=.01\).

\(_\blacksquare\)

1.2.9 Ejercicio 9 [Escrito por Ayrton Pablo Almada]

Realice la prueba de Bartlett para los datos que se dan en el archivo pregunta5.r. Realice la prueba dividiendo la población en 3 grupos del mismo tamaño, después realice la prueba dividiendo la población en 4 grupos del mismo tamaño y finalmente realice la prueba con \(n_i = \{49, 82, 103, 66\}\). Use \(\alpha = 0.05\).

\(\underline{\text{Solución:}}\)

Vamos a probar la homogeneidad en la varianza de \(r\in\{3,4\}\) muestras aleatorias del mismo tamaño \(c\), tendremos la siguiente prueba de hipótesis:

\[H_0:\sigma_1^2=\sigma_2^2=\cdots=\sigma_r^2~\text{ vs. }~H_a:\neg H_0.\] Para realizar la prueba vamos a ocupar la estadística \(T\) que sigue una distribución \(\chi^2_{(r-1)}\) la cual es de la forma:

\[ \begin{aligned} T&=\frac{(N-r)\text{ln}(S_p^2)-(c-1)\sum_{i=1}^r\text{ln}(S_i^2)}{1+\frac{1}{3(r-1)}\left(\frac{r}{c-1}-\frac{1}{N-r}\right)} :\\ S_p^2&=\frac{\sum_{i=1}^{r}SC_i}{N-r},~S_i^2=\frac{SC_i}{c-1},~SC_i=\sum_{j=1}^cx_{ij}^2-\frac{\left(\sum_{j=1}^{c}x_{ij}\right)^2}{c}. \end{aligned} \] Obtendremos esta estadística valuada en las muestras por medio del siguiente código:

#Cada columna de la Matriz con la que trabajaremos correspone a una 
#de las r muestras de c observaciones cada una. 
MatrixBase1<-matrix(as.numeric(M),ncol=3)
#i\in\{1,2,3\}
#c=100
MatrixBase2<-matrix(as.numeric(M),ncol=4)
#i\in\{1,2,3,4\}
#c=75
N1<-dim(MatrixBase1)[1]*dim(MatrixBase1)[2]
N2<-dim(MatrixBase2)[1]*dim(MatrixBase2)[2]
SC_i<-function(m,i){
  sum(m[,i]^2)-(1/dim(m)[1])*(sum(m[,i]))^2
}
S_i<-function(m,i){
  SC_i(m,i)/(dim(m)[1]-1)
}
S_p<-function(m){
  s<-0
  for (j in 1:dim(m)[2]) {
    s<-s+SC_i(m,j)
  }
  s/(dim(m)[1]*dim(m)[2]-dim(m)[2])
}
T<-function(m){
  N<-dim(m)[1]*dim(m)[2]
  e<-(N-dim(m)[2])*log(S_p(m))
  r<-0
  for (j in 1:dim(m)[2]) {
    n<-S_i(m,j)
    r<-r+-(dim(m)[1]-1)*log(n)
  }
  t<-1+(1/(3*(dim(m)[2]-1)))*(dim(m)[2]/(dim(m)[1]-1)-1/(N-dim(m)[2]))
  (e+r)/(t)
}

De forma tal que al evaluar \(T\) en cada muestra tendremos y al calcular los correspondientes cuantiles tendremos que:

a<-0.05
T1<-T(MatrixBase1)
T2<-T(MatrixBase2)
chi1<-qchisq(1-a,dim(MatrixBase1)[2]-1)
chi2<-qchisq(1-a,dim(MatrixBase2)[2]-1)
paste("T1=",T1,'<',chi1,'=chi1')
## [1] "T1= 1.7568778609265 < 5.99146454710798 =chi1"
paste("T2=",T2,'>',chi2,'=chi2')
## [1] "T2= 8.41361551906536 > 7.81472790325118 =chi2"

Por lo que para el caso donde la muestra sea de 3 grupos de mismo tamaño no tendremos razón para rechazar la hipótesis nula de homogeneidad de varianzas, por otro lado sí tenemos evidencia necesaria para rechazar la hipótesis nula de homogeneidad de varianzas en el caso de sea de 4 grupos de mismo tamaño.

Ahora, para el caso de tener una muestra aleatoria de cuatro grupos de tamaños \(n_i=\{49, 82, 103, 66\}\):

N<-300
R<-4
M1<-M[1:49]
M2<-M[50:131]
M3<-M[132:234]
M4<-M[235:300]
C<-c(49, 82, 103, 66)
SC_i<-function(m){
  sum(m^2)-(1/length(m))*(sum(m))^2
}
S_i<-function(m){
  SC_i(m)/(length(m)-1)
}
S_p<-(1/(N-R))*((C[1]-1)*S_i(M1)+(C[2]-1)*S_i(M2)+(C[3]-1)*S_i(M3)+(C[4]-1)*S_i(M4))
U<-(C[1]-1)*log(S_i(M1))+(C[2]-1)*log(S_i(M2))+(C[3]-1)*log(S_i(M3))+(C[4]-1)*log(S_i(M4))
V<-1/(C[1]-1)+1/(C[2]-1)+1/(C[3]-1)+1/(C[4]-1)
Tr<-((N-R)*log(S_p)-U)/(1+(1/(3*(R-1)))*(V-(1/(N-R))))
a<-0.05
chi1<-qchisq(1-a,dim(MatrixBase1)[2]-1)
paste("T=",Tr,'<',chi1,'=chi1')  
## [1] "T= 5.89802962875278 < 5.99146454710798 =chi1"

Por lo que para el caso donde la muestra sea de 4 grupos de tamaños \(n_i=\{49, 82, 103, 66\}\) no tendremos razón para rechazar la hipótesis nula de homogeneidad de varianzas.

\(_\blacksquare\)

2 Ejercicio Extra [Ejercicios.pdf] [Escrito por Ayrton Pablo Almada]

2.1 Pruebas de Bondad de Ajuste

2.1.1 Ejercicio 1

A continuación se muestra un registro de 100 mediciones de la velocidad de la luz en el aire. Este experimento fue realizado por un marino de Estados Unidos, el maestro Albert Michelson en 1879. ¿Los datos parecen tener distribución Normal? Utilice las pruebas de Kolmogorov, Anderson-Darling y Shapiro-Wilk. Compare sus respuestas.

Contamos con los siguientes datos:

D<-c(299.85,299.74,299.90,300.07,299.93,299.85,299.95,299.98,300.00,299.98,
     299.93,299.65,299.76,299.81,300.00,300.00,299.96,299.94,299.96,299.94,
     299.88,299.80,299.85,299.88,299.83,299.79,299.81,299.88,299.88,299.83,
     299.80,299.79,299.88,299.88,299.88,299.86,299.72,299.72,299.62,299.86,
     299.88,299.91,299.85,299.87,299.84,299.84,299.85,299.84,299.89,299.81,
     299.81,299.82,299.80,299.77,299.76,299.74,299.91,299.92,299.89,299.86,
     299.88,299.72,299.84,299.85,299.89,299.84,299.78,299.81,299.76,299.81,
     299.79,299.81,299.87,299.87,299.81,299.74,299.81,299.94,299.95,299.80,
     299.98,299.88,299.96,299.96,299.90,299.84,299.76,299.80,299.97,299.95,
     299.84,299.84,299.75,299.76,299.85,299.78,299.82,299.85,299.81,299.87)
s<-var(D)
m<-mean(D)
M<-matrix(D,nrow=10,ncol=10)
B<-data.frame(M)
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
299.85 299.93 299.88 299.80 299.88 299.81 299.88 299.79 299.98 299.84
299.74 299.65 299.80 299.79 299.91 299.82 299.72 299.81 299.88 299.84
299.90 299.76 299.85 299.88 299.85 299.80 299.84 299.87 299.96 299.75
300.07 299.81 299.88 299.88 299.87 299.77 299.85 299.87 299.96 299.76
299.93 300.00 299.83 299.88 299.84 299.76 299.89 299.81 299.90 299.85
299.85 300.00 299.79 299.86 299.84 299.74 299.84 299.74 299.84 299.78
299.95 299.96 299.81 299.72 299.85 299.91 299.78 299.81 299.76 299.82
299.98 299.94 299.88 299.72 299.84 299.92 299.81 299.94 299.80 299.85
300.00 299.96 299.88 299.62 299.89 299.89 299.76 299.95 299.97 299.81
299.98 299.94 299.83 299.86 299.81 299.86 299.81 299.80 299.95 299.87

Kolmogorov-Smirnoff

\(\underline{\text{Solución:}}\)

Deseamos probar

\[H_0: \text{La muestra std} \sim N(\mu=0, \sigma^2=1) ~\text{ vs. }~ H_a:\text{La muestra std} \nsim N(\mu=0, \sigma^2=1)\]

Dado que ya tenemos la muestra ordenada, calculemos la función de distribución empírica \(S_n(x_i)=\frac{i}{n}\), la función de distribucion empírica retrasada \(S_n^-(x_i)=\frac{i-1}{n}\), y la distribución propuesta, que en este caso de una normal estándar ocupando \(\mu=\bar{X}\) y \(\sigma^2=S^2\), donde \(\bar{X}\) es la media muestral y \(S^2\) es la varianza muestral.

s<-var(D) #Varianza muestral
m<-mean(D) #Media Muestral
D<-(sort(D)-m)/sqrt(s) #Arreglo de datos sorteado y estadarizado
D_<-sort(D,decreasing = TRUE) #Arreglo de datos sorteado decreciente y estadarizado
n<-length(D) #Cantidad de datos


#Funcion de distribucion empirica

Sn<-function(x){
  sum(ifelse(D<=x,1,0))/n
}


#Funcion de distribucion teorica normal

Fn<-function(x){
  pnorm(x)
}

SN<-c()
for (i in 1:n) {
  SN<-c(SN,Sn(D[i]))
}

FN<-Fn(D)

Mediante la comparación,se pueden realizar las diferencias en valor absoluto de la distribucion teórica y la distribucion empírica, es decir \(|S_n(X_i)-F(X_i)|\) y \(|S_n(X_{i-1})-F(X_i)|\)

diff_F<-c()
diff_F_retraso<-c()
for (i in 1:n){
  diff_F[i]=abs(Fn(D[i])-Sn(D[i]))
  
  }

diff_F_retraso[1]=Fn(D[1])
for (i in 2:n){
 
    diff_F_retraso[i]=abs(Fn(D[i])-Sn(D[i-1]))
  }

Una vez efectuados estos cálculos, se obtiene la siguiente tabla:

library(knitr)
resumen<-cbind(D,SN,FN,diff_F,diff_F_retraso)
colnames(resumen)<-c("X_i","Sn(z_i)","F(z_i)",
                     "Di+","Di-")
X_i Sn(z_i) F(z_i) Di+ Di-
-2.941379 0.01 0.0016338 0.0083662 0.0016338
-2.561683 0.02 0.0052083 0.0147917 0.0047917
-1.675726 0.05 0.0468960 0.0031040 0.0268960
-1.675726 0.05 0.0468960 0.0031040 0.0031040
-1.675726 0.05 0.0468960 0.0031040 0.0031040
-1.422595 0.08 0.0774268 0.0025732 0.0274268

Procedamos a calcular las estadísticas \(D^-=\max\{D_i^-\}\) y \(D^+=\max\{D_i^+\}\) para el contraste con los cuantiles de la distribución:

## [1] "D+ = 0.0834243742739519"
## [1] "D- = 0.0834243742739519"

Por lo tanto, \(D=\max \{D^+,D^-\}=\max \{0.0834243,0.0834243\}=0.0834243\). Este resultado debe ser comparado con los valores críticos de Lilliefors. Para un nivel de significancia \(\alpha=5\%\) y el tamaño de muestra \(n=100\), tenemos que \(W_{\alpha}=W_{0.05}=\frac{0.875897}{\sqrt{100}}=0.0875897\)

\[\therefore~ D=0.0834243<0.0875897=W_{0.05} \] De esta manera,tenemos evidencia suficiente para no rechazar la hipótesis nula de normalidad de la muestra.

Ahora, efectuemos la prueba mediante la función lillie.test():

lillie.test(D)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  D
## D = 0.083424, p-value = 0.08289

Nótese que hemos obtenido un \(\text{p-value} = 0.08289>\alpha=0.05\), por lo que tenemos mayor razón para no rechazar la hipótesis nula de normalidad de la muestra.

\(_\blacksquare\)

Anderson-Darling

\(\underline{\text{Solución:}}\)

Nuevamente, ocuparemos los datos estandarizados para facilitar los cálculos al ocupar una distribución normal std. Dado que ya tenemos los datos ajustados, procedemos a calcular la estadística de Anderson-Darling, (donde \(F^*\) es la distribución teórica propuesta) con \(\alpha=0.05\):

\[A_n^2=-n-\sum_{i=1}^n\left(\frac{2i-1}{n}\right)\left[\text{ln}(F^*(X_i))+\text{ln}(1-F^*(X_{n-i+1}))\right]\]

A=-n-sum((2*(seq(1,n,1)-1)/n)*(log(Fn(D))+log(1-Fn(D_))))
paste("A_n^2 =",A)
## [1] "A_n^2 = -1.53183227882967"

Debido a que \(A_n^2=-1.531832<2.492=W^{0.95}\), entonces tenemos evidencia suficiente para no rechazar la hipótesis nula de normalidad de la muestra.

Ahora, efectuemos la prueba mediante la función ad.test():

ad.test(D)
## 
##  Anderson-Darling normality test
## 
## data:  D
## A = 0.46076, p-value = 0.255

Nótese que hemos obtenido un \(\text{p-value} = 0.255>\alpha=0.05\), por lo que tenemos mayor razón para no rechazar la hipótesis nula de normalidad de la muestra.

\(_\blacksquare\)

Shapiro-Wilk

\(\underline{\text{Solución:}}\)

Nuevamente, ocuparemos la misma prueba de hipótesis de los anteriores incisos de este ejercicio.

Debido a la complejidad de la simulación de las medias de los estadísticos de orden, aplicaremos la prueba de hipótesis por medio de su función en R.

shapiro.test(D)
## 
##  Shapiro-Wilk normality test
## 
## data:  D
## W = 0.98807, p-value = 0.5137

Como tenemos un \(\text{p-value}=0.5137>0.05=\alpha\), podemos concluír que para un nivel de significancia \(\alpha=0.05\) no rechazamos la hipótesis nula de normalidad de la muestra.

\(_\blacksquare\)

2.2 Pruebas de Proporciones

2.2.1 Ejercicio 1

Se sabe que el 20% de ciertas especies de insectos exhiben una característica particular A. 18 insectos de esa especie se obtienen de un ambiente inusual, y ninguno de ellos tiene la característica A. ¿Es razonable suponer que los insectos de ese ambiente tienen la misma probabilidad de 0.20 que tiene la especie en general? Utilice una prueba de dos colas.

\(\underline{\text{Solución:}}\)

Vamos a considerar una muestra aleatoria \(\{X_i\}_{i=1}^{n=18}\) con distibución Bernoulli de parámetro \(p\), (suponemos esta distribución ya que se nos habla de dos clases de insectos: ‘aquellos con carácterística A’ ó ‘aquellos sin carácterística A’, por lo que \(p\) representa la probabilidad de un insecto presente la característica A), por lo tanto nuestra prueba de hipótesis será de la forma:

\[H_o:~p= p^*=0.2~\text{vs.}~H_a:~p\ne p^*=0.2\] Por lo que trabajaremos con el siguinete estadístico \(T\):

\[\forall_{j,i=1,...,n}~X_i\sim\text{Blli}(p)~,~X_i\perp X_j~~\text{si}~i\ne j\implies T=\sum_{i=1}^{n}\mathbb{I}_{\{X_i=1\}}\sim\text{Binom}(n,p)\]

Por ende, suponiendo válida la hipóteis nula tendremos que \(T=\sum_{i=1}^{n}\mathbb{I}_{\{X_i=1\}}\sim\text{Binom}(18,0.2)\). Por lo que, sea \(Y\sim\text{Binom}(18,0.2)\), dado que \(T=0\) y como \(\alpha=0.10\), vamos a buscar cuantiles \(t_1,t_2\in\mathbb{Z}\) tale que:

\[\mathbb{P}[Y\le t_1]=0.05=\alpha_1~,~\mathbb{P}[Y\le t_2]=1-\alpha_1=0.95\] Obteniendo así:

t_1<-qbinom(0.05,18,0.2)
t_2<-qbinom(0.95,18,0.2)
paste("t_1= ",t_1)
## [1] "t_1=  1"
paste("t_2= ",t_2)
## [1] "t_2=  7"

Rechazaremos la hipótesis nula con un nivel de confianza del \(0.90\) si \(T\le t_1\) ó si \(T > t_2\); dado que \(T=0<i=t_1\), rechazamos la idea de suponer que los insectos de ese ambiente tienen la misma probabilidad de 0.20 que tiene la especie en general. Automatizando lo anteriormente realizado:

prop.test(0,18,p=0.2,alternative=c("two.sided"),conf.level=0.90)
## Warning in prop.test(0, 18, p = 0.2, alternative = c("two.sided"),
## conf.level = 0.9): Chi-squared approximation may be incorrect
## 
##  1-sample proportions test with continuity correction
## 
## data:  0 out of 18, null probability 0.2
## X-squared = 3.3368, df = 1, p-value = 0.06775
## alternative hypothesis: true p is not equal to 0.2
## 90 percent confidence interval:
##  0.0000000 0.1751337
## sample estimates:
## p 
## 0

Dado que el p-value calculado es de \(\text{p-value} = 0.06775<0.1=\alpha\), tenemos razón suficiente para descartar la hipótesis nula con un nivel de condianza de 0.9.

\(_\blacksquare\)

2.3 Pruebas de Rachas

2.3.1 Ejercicio 2

Los siguientes son 30 lapsos en minutos entre las erupciones del volcán Old Faithful Geyaer en el Parque Nacional de Yellowstone, fueron grabados entre las 8:00 a.m. y las 10:00 p.m. de un día determinado, medidos desde el comienzo de una erupción hasta el inicio de la siguiente.

68 63 66 63 61 44 60 62 71 62 62 55 62 67 73 72 55 67 68 65 60 61 71 60 68 67 72 69 65 66

Un investigador desea utilizar estos datos con fines de inferencia, pero le preocupa si es razonable tratar dichos datos como una muestra aleatoria. ¿Debería considerarla aleatoria? Justifica tu respuesta.

\(\underline{\text{Solución:}}\)

Dado que se tiene múltiples datos y la prueba de corridas es dicotómica, se procede a categorizar los resultados para formar dos subclases ambas excluyentes la una de la otra. Para se escoge como medida la mediana, en este caso la mediana de la muestra es 65.

Los valores que son menores a la mediana se les asignará la letra ‘p’, y los que son mayores a la mediana se les asignará la letra ‘q’, de esta manera con la nueva asignación se genera una muestra aleatoria dicotómica.

Lapsos<-c(68,63,66,63,61,44,60,62,71,62,62,55,62,67,73,72,
          55,67,68,65,60,61,71,60,68,67,72,69,65,66)
Mediana<-median(Lapsos)#65
Clases<-dplyr::case_when(Lapsos>Mediana~'q',Lapsos<=Mediana~'p',TRUE ~ as.character(Lapsos))
Clases
##  [1] "q" "p" "q" "p" "p" "p" "p" "p" "q" "p" "p" "p" "p" "q" "q" "q" "p"
## [18] "q" "q" "p" "p" "p" "q" "p" "q" "q" "q" "q" "p" "q"

Se tienen contabilizan las frecuencias de las variables dicotómicas, de forma tal que \(n_1=14,n_2=16\).

Con un nivel de significancia \(\alpha=0.05\) y con las subrachas \(n_1=14\) y \(n_2=16\), se procede a buscar las estadísticas \(W^{0.025}\) y \(W^{0.975}\), de manera tal que (si \(n_1,n_2<20\)):

W_1=randtests::qruns(0.025,14,16)
W_2=randtests::qruns(0.975,14,16)

paste("W_1 = ",W_1)
## [1] "W_1 =  11"
paste("W_2 = ",W_2)
## [1] "W_2 =  21"

Dado que \(n_1=14,n_2=16\), tenemos que:

\[R=15>11=W^{0.025}~\text{y}~R=15<21=W^{0.975}\]

tseries::runs.test(as.factor(Clases))
## 
##  Runs Test
## 
## data:  as.factor(Clases)
## Standard Normal = -0.34844, p-value = 0.7275
## alternative hypothesis: two.sided

De esta manera la estadística R no cae en la región de rechazo, por lo que no se rechaza la hipótesis nula, asumiéndose así con un nivel de significancia del 5 %, que el organizador efectivamente realizo la rifa de manera aleatoria. Eso y además de que el p_value es del \(0.7275>\alpha=0.05\) \(_\blacksquare\)

2.4 Medidas de Correlación

2.4.1 Ejercicio 1

Dos enólogos, X y Y, se les pidió que clasificarán 8 vinos catados de mejor a peor (rango 1=mejor, rango 8= peor). Encuentre el coeficiente de Correlación de Spearman entre los expertos. Si el tamaño de la muestra incrementa a N=80 y se encuentra que \(\hat{\rho}\) es 10 veces más pequeño que el que encontró antes, ¿Cuál sería el p-value para la prueba de hipótesis de dos colas?

Marca de Vino A B C D E F G H
Experto X 1 2 3 4 5 6 7 8
Experto Y 2 3 1 4 7 8 5 6

\(\underline{\text{Solución:}}\)

Vamos a obtener el l coeficiente de Correlación de Spearman entre los expertos:

\[\hat{\rho}=\frac{\sum_{i=0}^{n}\left(R(X_i)-\frac{n+1}{2}\right)\left(R(Y_i)-\frac{n+1}{2}\right)}{\left[\frac{n(n^2-1)}{12}\right]}\] Es decir:

A<-c(1,2,3,4,5,6,7,8)
B<-c(2,3,1,4,7,8,5,6)
I<-rbind(A,B)
n<-dim(I)[1]*dim(I)[2]
i<-sum((I[1,]-(n+1)/2)*(I[2,]-(n+1)/2))
k<-(n)*(n^2-1)/12
rho<-i/k
paste('rho^= ',rho)
## [1] "rho^=  0.467647058823529"

Ahora, dado que se supuso que si i el tamaño de la muestra incrementa a \(n=80\) entonces \(\hat{\rho}\) es 10 veces más pequeño que el que encontró antes, eso es \(\hat{\rho}=\frac{0.294117647058824}{10}=0.02941176\). Procederemos a analizar la siguiente prueba de hipótesis:

\[H_0:\text{Las } X_i\text{'s y las }Y_i\text{'s son mutuamente independientes vs. } H_1:\neg H_0\] Dado que la misma \(\hat{\rho}\) de Spearman funge como estadístico de esta prueba, obtendremo los cuantiles con un nivel de significancia \(\alpha=0.05\)

q1<-qSpearman(0.025,80)
q2<-qSpearman(0.975,80)
paste('w_{0.025}= ', q1)
## [1] "w_{0.025}=  -0.219878105954055"
paste('w_{0.975}= ', q2)
## [1] "w_{0.975}=  0.2199484294421"

Ahora, debido a que \(-0.2198781=w_{0.025}\le \hat{\rho}=0.02941176\le w_{0.975}=0.2199484\), entonces no tenemos mayor razón para rechazar la hipótesis de independencia mutua entre las muestras. Procederemos a calcular el p-value para este caso:

\[\text{p-value}=\mathbb{P}\left[|\hat{\rho}|>\hat{\rho}_{eva}=0.02941176\right]=2\min{\{\mathbb{P}\left[\hat{\rho}\ge 0.02941176\right],\mathbb{P}\left[\hat{\rho}\le0.02941176\right]\}}\]

p1<-pSpearman(0.02941176,80)
p2<-pSpearman(0.02941176,80,lower.tail = F)
p_value=2*min(p1,p2)
paste('p-value= ',p_value)
## [1] "p-value=  0.795922848322558"

Por lo tanto, el p-value para la prueba de hipótesis de dos colas y el tamaño de la muestra incrementa a N=80 y se encuentra que \(\hat{\rho}\) es 10 veces más pequeño que el que encontró antes es de 0.7959228.\(_\blacksquare\)