Lección 7 Contrastes de independencia y homogeneidad
El test \(\chi^2\) explicado en la lección anterior permite contrastar, en situaciones adecuadas, si una muestra proviene de una determinada distribución y por lo tanto también si una muestra sigue la misma distribución que otra muestra. Como, en última instancia, la independencia de dos variables cualitativas se puede describir en términos de igualdades de probabilidades, el test \(\chi^2\) también nos permitirá contrastar si dos variables cualitativas son independientes, tanto en el sentido de que las probabilidades conjuntas sean el producto de las probabilidades marginales (con un contraste de independencia) como en el sentido de que las distribuciones condicionadas de una respecto de los valores de la otra sean todas iguales (con un contraste de homogeneidad).
Aunque, como veremos, los contrastes de independencia y homogeneidad son idénticos desde el punto de vista matemático e incluso utilizan el mismo estadístico \(\chi^2\) y la misma definición de p-valor, provienen de diseños experimentales diferentes:
En un contraste de independencia se toma una muestra transversal de la población, es decir, se selecciona al azar una cierta cantidad de individuos de la población, se observan las dos variables sobre cada uno de ellos, y se contrasta si las probabilidades conjuntas son iguales al producto de las probabilidades marginales de cada variable. Formalmente, si \(X\) e \(Y\) son las dos variables, se contrasta si para cada par de posibles valores \(x\) de \(X\) e \(y\) de \(Y\) se tiene que \[ P(X=x,Y=y)=P(X=x)\cdot P(Y=y) \] o si por el contrario hay algún par de valores \(x\), \(y\) para los que esta igualdad sea falsa.
En un contraste de homogeneidad se escoge una de las variables y para cada uno de sus posibles valores se toma una muestra aleatoria, de tamaño prefijado, de individuos con ese valor para esa variable; su unión forma una muestra estratificada en el sentido de la Sección 2.1. A continuación, se observa sobre cada uno de estos individuos la otra variable. En esta situación contrastamos si la distribución de probabilidades de la segunda variable es la misma en los diferentes estratos definidos por los niveles de la primera variable. Formalmente, si \(Y\) es la variable que usamos en primer lugar para clasificar los individuos de la población y tomar una muestra de cada clase, con posibles valores \(y_1,\ldots,y_k\), y \(X\) es la variable que medimos en segundo lugar sobre los individuos escogidos, se contrasta si, para cada posible valor \(x\) de \(X\), \[ P(X=x|Y=y_1)=P(X=x|Y=y_2)=\cdots=P(X=x|Y=y_k) \] o si por el contrario existen \(x\), \(y_i\), \(y_j\) tales que \(P(X=x|Y=y_i)\neq P(X=x|Y=y_j)\).
En ambos contrastes, la hipótesis nula es que las variables son independientes, bajo una u otra formulación matemática, y la hipótesis alternativa es que son dependientes (o que hay asociación entre ellas). La hipótesis nula se rechaza si se obtiene evidencia que hace inverosímiles las igualdades de probabilidades que se contrastan.
Para ilustrar esta lección, hemos generado una muestra aleatoria de cadenas formadas por las bases “a”, “c”, “g” y “t”. En concreto, hemos generado cadenas de longitud 100 de tres tipos: A, B y C. Estos tipos se distinguen por los vectores de probabilidades que han determinado las frecuencias de las cuatro bases en las secuencias. Queremos investigar si hay relación entre el tipo (A, B o C) de una cadena, y la base de frecuencia máxima en ella.
Los datos y el método de generación se encuentran en el repositorio siguiente: https://github.com/biocom-uib/Experimento-Cadenas. Este directorio contiene:
El fichero LeemeGeneracionDatos, que explica cómo se han generado las muestras.
El fichero MuestraTotalBases.txt, que contiene una tabla de datos de 10000 observaciones de las dos variables siguientes sobre cadenas: el
tipo, que es un factor con los nivelesA,ByC, ymax.frec, que es otro factor que indica qué base tiene mayor frecuencia en la cadena. Este fichero es de formato texto, con una primera fila con el nombre de las variables y sus columnas separadas por comas
El código siguiente carga la tabla de datos, comprueba que no ha habido problemas, y extrae tres subtablas, una para cada tipo de cadena.
poblacion=read.table("https://raw.githubusercontent.com/biocom-uib/Experimento-Cadenas/master/MuestraTotalBases.txt",
header=TRUE,sep=",")
str(poblacion)## 'data.frame': 10000 obs. of 2 variables:
## $ tipo : Factor w/ 3 levels "A","B","C": 2 1 1 3 1 1 1 3 2 2 ...
## $ max.frec: Factor w/ 4 levels "a","c","g","t": 3 1 2 4 4 1 1 3 2 2 ...head(poblacion)## tipo max.frec
## 1 B g
## 2 A a
## 3 A c
## 4 C t
## 5 A t
## 6 A apoblacionA=subset(poblacion,tipo=="A")
poblacionB=subset(poblacion,tipo=="B")
poblacionC=subset(poblacion,tipo=="C")7.1 Tablas de contingencia
Ya estudiamos en la Lección ?? de la primera parte las tablas de contingencia. En esta sección vamos a repasar y ampliar algunas de las funciones de R para el manejo de esta clase de tablas.
La tabla de contingencia de frecuencias absolutas conjuntas de las dos variables del data frame poblacion se calcula de la manera siguiente:
tabla=table(poblacion$tipo,poblacion$max.frec)
tabla##
## a c g t
## A 1027 1003 998 998
## B 1259 1297 209 232
## C 245 602 1453 677Su tabla de frecuencias relativas conjuntas en el total de la muestra, en términos de proporciones (tantos por uno), es:
prop.table(tabla)##
## a c g t
## A 0.1027 0.1003 0.0998 0.0998
## B 0.1259 0.1297 0.0209 0.0232
## C 0.0245 0.0602 0.1453 0.0677Para añadir las distribuciones marginales de la tabla de contingencia (o márgenes de
la tabla), se añade una nueva fila con las sumas de cada columna y una nueva columna con las sumas de cada fila. Con R, esto se puede llevar a cabo fácilmente con la función
addmargins. Su sintaxis básica es
addmargins(tabla, margin=..., FUN=...)donde:
tablaes unatable.margines un parámetro que puede tomar los valores siguientes:1si queremos una nueva fila con las marginales de cada columna.2si queremos una nueva columna con las marginales de cada fila.c(1,2), que es el valor por defecto para tablas de contingencia bidimensionales (y por lo tanto no hace falta especificarlo), si queremos las marginales por filas y por columnas.7
FUNes la función que se aplica a las filas o columnas para obtener el valor marginal. Por defecto es la suma, que es la función que nos interesa en esta lección, y por tanto tampoco hace falta especificarlo.
El resultado es otro objeto de la clase table al que se le han añadido una
o varias filas o columnas. Éstas contienen los márgenes resultantes
de aplicar la función indicada por FUN. La etiqueta de las nuevas filas o columnas es la
función que se aplica.
Por ejemplo, para obtener las tablas marginales completas de frecuencias absolutas y relativas en nuestro ejemplo, haríamos:
addmargins(tabla) ##
## a c g t Sum
## A 1027 1003 998 998 4026
## B 1259 1297 209 232 2997
## C 245 602 1453 677 2977
## Sum 2531 2902 2660 1907 10000addmargins(prop.table(tabla)) ##
## a c g t Sum
## A 0.1027 0.1003 0.0998 0.0998 0.4026
## B 0.1259 0.1297 0.0209 0.0232 0.2997
## C 0.0245 0.0602 0.1453 0.0677 0.2977
## Sum 0.2531 0.2902 0.2660 0.1907 1.0000También podemos calcular la tabla de proporciones por filas y con su marginal por filas comprobar que efectivamente la suma de cada fila es 1:
addmargins(prop.table(tabla,margin=1),margin=2)##
## a c g t Sum
## A 0.25509190 0.24913065 0.24788872 0.24788872 1.00000000
## B 0.42008675 0.43276610 0.06973640 0.07741074 1.00000000
## C 0.08229762 0.20221700 0.48807524 0.22741014 1.00000000Y viceversa, podemos calcular la tabla de proporciones por columnas y con su marginal por columnas comprobar que efectivamente la suma de cada columna es 1:
addmargins(prop.table(tabla,margin=2),margin=1)##
## a c g t
## A 0.40576847 0.34562371 0.37518797 0.52333508
## B 0.49743185 0.44693315 0.07857143 0.12165705
## C 0.09679968 0.20744314 0.54624060 0.35500787
## Sum 1.00000000 1.00000000 1.00000000 1.00000000Observad que el significado de margin en addmargins es diferente de, por ejemplo, en prop.table o en apply: en estas dos últimas instrucciones indica la dimensión en la que calculamos las proporciones o aplicamos la función, mientras que en addmargins indica la dimensión en la que añadimos el margen, que, por lo tanto, se calcula aplicando la función en la otra dimensión.
Si sólo nos interesa la fila o la columna de marginales, podemos usar las instrucciones
colSums y rowSums, que suman una tabla por columnas y por filas, respectivamente. Por ejemplo, para obtener los vectores de marginales por columnas y por filas, respectivamente, podríamos entrar:
colSums(tabla)## a c g t
## 2531 2902 2660 1907rowSums(tabla)## A B C
## 4026 2997 2977También podemos obtener estos márgenes extrayendo los márgenes de la tabla con márgenes, obtenida aplicando addmargins a la tabla original, por medio de las instrucciones usuales para extraer filas y columnas.
addmargins(tabla)["Sum",-dim(addmargins(tabla))[2]]## a c g t
## 2531 2902 2660 1907addmargins(tabla)[-dim(addmargins(tabla))[1],"Sum"]## A B C
## 4026 2997 2977Observad, por ejemplo, la construcción
addmargins(tabla)["Sum",-dim(addmargins(tabla))[2]]Con addmargins(tabla)["Sum",] obtendríamos la fila Sum de la tabla con las marginales, incluyendo la última entrada, correspondiente a la columna Sum. Por lo tanto, hay que eliminarla: como dim(addmargins(tabla))[2] es el número de columnas de la tabla addmargins(tabla), es decir, la longitud del vector addmargins(tabla)["Sum",], la última entrada (correspondiente a la última columna) se puede eliminar especificando -dim(addmargins(tabla))[2] en las columnas al extraer la fila Sum.
7.2 Contraste de independencia
El contraste de independencia para tablas de contingencia bidimensionales consiste en decidir si las dos variables de la tabla tienen distribuciones independientes, es decir, si la distribución de probabilidades conjunta es igual al producto de las probabilidades marginales. En nuestro ejemplo, queremos decidir si podemos aceptar que las variables tipo y max.frec son independientes o si por el contrario hay evidencia de que la distribución de las bases de máxima frecuencia depende del tipo de cadena.
Vamos a extraer una muestra aleatoria simple de la población y observar los valores de las dos variables. En concreto, seleccionaremos una muestra transversal de 150 filas, al azar y con reposición, de entre las 10000 filas del data frame poblacion. El código es el siguiente (fijamos la semilla de aleatoriedad para que sea reproducible):
set.seed(42)
n=150
indices.muestra=sample(1:10000, size=n, replace=TRUE)
muestra.test.indep= poblacion[indices.muestra, ] #Las filas que forman la muestraAhora calculamos la tabla de contingencia con sus marginales.
tabla.ind=table(muestra.test.indep$tipo, muestra.test.indep$max.frec)
tabla.ind##
## a c g t
## A 8 17 16 12
## B 23 22 3 7
## C 4 9 17 12tabla.ind.marg=addmargins(tabla.ind)
tabla.ind.marg##
## a c g t Sum
## A 8 17 16 12 53
## B 23 22 3 7 55
## C 4 9 17 12 42
## Sum 35 48 36 31 150Extraemos sus dos márgenes. Las frecuencias marginales de las filas:
frec.abs.tipo=tabla.ind.marg[-dim(tabla.ind.marg)[1],"Sum"]
frec.abs.tipo## A B C
## 53 55 42Las frecuencias marginales de las columnas:
frec.abs.max.frec=tabla.ind.marg["Sum",-dim(tabla.ind.marg)[2]]
frec.abs.max.frec## a c g t
## 35 48 36 31El test de independencia usa las frecuencias absolutas esperadas bajo la hipótesis nula de independencia, que se obtienen, para cada celda (i,j), multiplicando la frecuencia marginal de la fila i por la de la columna j y dividiendo por el tamaño de la muestra. En nuestro ejemplo estas frecuencias esperadas son \[ \begin{array}{l|cccc} & \mbox{a} & \mbox{c} & \mbox{g} & \mbox{t} \\ \hline \mbox{A} & 53\cdot 35/150 & 53\cdot 48/150 & 53\cdot 36/150 & 53\cdot 31/150 \\ \mbox{B} & 55\cdot 35/150 & 55\cdot 48/150 & 55\cdot 36/150 & 55\cdot 31/150 \\ \mbox{C} & 42\cdot 35/150 & 42\cdot 48/150 & 42\cdot 36/150 & 42\cdot 31/150 \end{array} \] y podemos obtenerlas fácilmente mediante un producto de matrices: \[ \frac{1}{150}\cdot\left(\begin{array}{c} 53 \\ 55 \\ 42\end{array}\right)\cdot \big( 35 ,48 , 36 , 31\big). \] Por lo tanto, con R obtenemos esta tabla de frecuencias esperadas de la manera siguiente:
frec.esperadas=frec.abs.tipo%*%t(frec.abs.max.frec)/n
frec.esperadas## a c g t
## [1,] 12.36667 16.96 12.72 10.95333
## [2,] 12.83333 17.60 13.20 11.36667
## [3,] 9.80000 13.44 10.08 8.68000Aunque vayamos a realizar el test de independencia con una función de R, es necesario comprobar que todas estas frecuencias esperadas son mayores o iguales que 5, por lo que no podemos evitar este cálculo. En este caso vemos que se cumple esta condición.
Si queremos realizar el test \(\chi^2\) de independencia a mano, podemos calcular el estadístico de forma directa con
chi2.estadistico=sum((tabla.ind-frec.esperadas)^2/frec.esperadas)
chi2.estadistico## [1] 32.12115y el p-valor del contraste, con
p.valor=1-pchisq(chi2.estadistico,df=(dim(tabla.ind)[1]-1)*(dim(tabla.ind)[2]-1))
p.valor## [1] 1.546778e-05Para realizar el test \(\chi^2\) de independencia con R, es suficiente aplicar la función chisq.test a la tabla de contingencia de frecuencias absolutas:
chisq.test(tabla.ind)##
## Pearson's Chi-squared test
##
## data: tabla.ind
## X-squared = 32.121, df = 6, p-value = 1.547e-05Como el p-valor es muy pequeño, podemos rechazar la hipótesis de que las variables objeto de estudio sean independientes: hemos obtenido evidencia estadísticamente significativa de que la distribución de las bases de frecuencia máxima sí que depende del tipo de cadena.
Si algunas frecuencias absolutas esperadas fueran inferiores a 5, la
aproximación del p-valor por una distribución \(\chi^2\) podría no ser adecuada.
En este caso, al ser las variables cualitativas, no podemos recurrir al
agrupamiento de valores consecutivos, puesto que no tienen orden. Si se da esta
situación, lo mejor es recurrir a simular el p-valor usando el parámetro
simulate.p.value=TRUE.
Por ejemplo, consideremos la situación siguiente, en la que tomamos una muestra de solo 100 cadenas y con otra semilla de aleatoriedad:
set.seed(300)
n2=100
indices.muestra2=sample(1:10000,size=n2,replace=TRUE)
muestra.test.indep2= poblacion[indices.muestra2,]
tabla.ind2=table(muestra.test.indep2$tipo,muestra.test.indep2$max.frec)
tabla.ind2##
## a c g t
## A 13 10 11 12
## B 9 18 3 2
## C 4 3 10 5Si aplicamos a esta tabla la función chisq.test, obtenemos:
chisq.test(tabla.ind2)## Warning in chisq.test(tabla.ind2): Chi-squared approximation may be
## incorrect##
## Pearson's Chi-squared test
##
## data: tabla.ind2
## X-squared = 21.843, df = 6, p-value = 0.001293¡Vaya! Veamos la tabla de frecuencias esperadas:
frec.abs.tipo2=rowSums(tabla.ind2)
frec.abs.max.frec2=colSums(tabla.ind2)
frec.esperadas2=frec.abs.tipo2%*%t(frec.abs.max.frec2)/n2
frec.esperadas2## a c g t
## [1,] 11.96 14.26 11.04 8.74
## [2,] 8.32 9.92 7.68 6.08
## [3,] 5.72 6.82 5.28 4.18Hay una frecuencia esperada inferior a 5. Por lo tanto, lo recomendable es estimar el p-valor del test \(\chi^2\) de independencia mediante simulaciones. Pero ahora tenemos que ir con cuidado en una cosa: hemos fijado la semilla de aleatoriedad para obtener una muestra de cadenas con frecuencias esperadas inferiores a 5. Lo recomendable es reiniciar esta semilla a un valor aleatorio con set.seed(NULL).
set.seed(NULL)
chisq.test(tabla.ind2,simulate.p.value=TRUE,B=5000)$p.value## [1] 0.00119976chisq.test(tabla.ind2,simulate.p.value=TRUE,B=5000)$p.value## [1] 0.00239952chisq.test(tabla.ind2,simulate.p.value=TRUE,B=5000)$p.value## [1] 0.00079984chisq.test(tabla.ind2,simulate.p.value=TRUE,B=5000)$p.value## [1] 0.0009998El p-valor es sistemáticamente pequeño, lo que nos permite rechazar la hipótesis de que las variables son independientes.
7.3 Contraste de homogeneidad
Como ya hemos dicho, la diferencia entre el contraste de homogeneidad y el de independencia está en el diseño del experimento: en cada contraste se selecciona la muestra de una manera diferente.
En nuestro caso, para contrastar si la distribución de probabilidades de la base de mayor frecuencia es la misma para cada tipo de cadena o no, lo que vamos a hacer es tomar una muestra aleatoria de 50 cadenas de cada tipo, juntarlas en una sola muestra estratificada, y aplicar el test \(\chi^2\) a esta muestra. El código siguiente realiza el muestreo en cada subpoblación de tipo y guarda la muestra total en el data frame muestra.test.homo. Fijamos de nuevo la semilla de aleatoriedad (otra), para que el test sea reproducible.
set.seed(100)Generamos los vectores de índices de las muestras:
n3=50
indices.muestraA=sample(1:dim(poblacionA)[1],size=n3,replace=TRUE)
indices.muestraB=sample(1:dim(poblacionB)[1],size=n3,replace=TRUE)
indices.muestraC=sample(1:dim(poblacionC)[1],size=n3,replace=TRUE)Finalmente, tomamos las filas de cada muestra y las combinamos en un data frame:
muestraA.50=poblacionA[indices.muestraA,]
muestraB.50=poblacionB[indices.muestraB,]
muestraC.50=poblacionC[indices.muestraC,]
muestra.test.homo=rbind(muestraA.50,muestraB.50,muestraC.50)
str(muestra.test.homo)## 'data.frame': 150 obs. of 2 variables:
## $ tipo : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
## $ max.frec: Factor w/ 4 levels "a","c","g","t": 1 3 1 4 4 3 3 1 4 2 ...Calculamos la tabla de contingencia de la muestra:
tabla.homo=table(muestra.test.homo$tipo,muestra.test.homo$max.frec)
tabla.homo##
## a c g t
## A 13 12 16 9
## B 26 18 3 3
## C 3 10 28 9Añadimos los márgenes:
addmargins(tabla.homo)##
## a c g t Sum
## A 13 12 16 9 50
## B 26 18 3 3 50
## C 3 10 28 9 50
## Sum 42 40 47 21 150Confirmamos que hemos tomado 50 cadenas de cada grupo. Ahora calculamos las frecuencias esperadas bajo la hipótesis nula, para comprobar si son todas mayores o iguales que 5:
frec.abs.tipo=rowSums(tabla.homo)
frec.abs.tipo## A B C
## 50 50 50frec.abs.max.frec=colSums(tabla.homo)
frec.abs.max.frec## a c g t
## 42 40 47 21frec.esperadas=frec.abs.tipo%*%t(frec.abs.max.frec)/sum(frec.abs.tipo)
frec.esperadas ## a c g t
## [1,] 14 13.33333 15.66667 7
## [2,] 14 13.33333 15.66667 7
## [3,] 14 13.33333 15.66667 7Todas las frecuencias son mayores o iguales que 5, así que aplicamos la función chisq.test sin simular el p-valor:
chisq.test(tabla.homo)##
## Pearson's Chi-squared test
##
## data: tabla.homo
## X-squared = 44.986, df = 6, p-value = 4.71e-08El p-valor es muy pequeño, por lo que podemos rechazar que las distribuciones de los valores de las bases de máxima frecuencia sean la misma para cada valor de la variable tipo.
En definitiva, el tipo de cadena afecta a la distribución de la base de mayor frecuencia. Es la misma conclusión a la que habíamos llegado con el test de independencia, solo que ahora hemos realizado un tipo de experimento diferente.
7.4 Potencia de un contraste \(\chi^2\)
La potencia de un contraste \(\chi^2\), tanto de bondad de ajuste como de independencia o de homogeneidad, se puede calcular de manera similar a cómo lo hacíamos en otros tipos de contrastes de uno o dos parámetros. La instrucción para llevarlo a cabo es pwr.chisq.test del paquete pwr. Su sintaxis básica es
pwr.chisq.test(N=..., df=..., sig.level=..., w=..., power=...)donde:
Nes el tamaño de la muestra.dfes el número de grados de libertad del estadístico (recordad que en un test de bondad de ajuste es el número de clases menos 1 y menos el número de parámetros estimados, y en un test de independencia o de homogeneidad es el número de niveles de una variable menos 1 por el número de niveles de la otra variable menos 1).sig.leveles el nivel de significación \(\alpha\).wes la magnitud del efecto, que en este tipo de tests se define como \(\sqrt{X^2/N}\), siendo \(X^2\) el valor del estadístico de contraste y \(N\) el tamaño de la muestra completa.poweres la potencia \(1-\beta\).
Si se especifican todos estos parámetros menos uno, la función da el valor del parámetro que falta. Normalmente, querremos saber la potencia de un contraste a posteriori o el tamaño de la muestra necesario para tener la potencia deseada para una magnitud del efecto esperada concreta.
Veamos algunos ejemplos de uso de esta función.
Ejemplo 7.1 Vamos a calcular la potencia del contraste del Ejemplo 6.3. En ese ejemplo, el tamaño de la muestra fue \(N=40\), el número de grados de libertad fue 5 y obtuvimos que \(X^2=7.7\), por lo que la magnitud del efecto fue \(w=\sqrt{7.7/40}\). Tomaremos el nivel de significación usual, \(\alpha=0.05\).
library(pwr)
pwr.chisq.test(N=40, df=5, sig.level=0.05, w=sqrt(7.7/40))##
## Chi squared power calculation
##
## w = 0.4387482
## N = 40
## df = 5
## sig.level = 0.05
## power = 0.545751
##
## NOTE: N is the number of observationsLa potencia del contraste ha sido de, aproximadamente, un 55%.
Para obtener solo la potencia, podemos usar el sufijo $power:
pwr.chisq.test(N=40, df=5, sig.level=0.05, w=sqrt(7.7/40))$power## [1] 0.545751Ejemplo 7.2 Vamos a calcular la potencia del contraste de normalidad de las longitudes de los sépalos de flores iris del Ejemplo 6.8. En ese ejemplo el tamaño de muestra fue \(N=150\); como usamos 7 clases, pero estimamos 2 parámetros, el número de grados de libertad fue 4; obtuvimos que \(X^2=11.0637\), por lo que \(w=\sqrt{11.0637/150}\); y ahora, por variar, tomaremos \(\alpha=0.1\).
pwr.chisq.test(N=150, df=4, sig.level=0.1, w=sqrt(11.0637/15))$power## [1] 1La potencia da 1: la probabilidad de que aceptáramos que la muestra seguía una distribución normal si no fuera verdad es prácticamente 0.
Ejemplo 7.3 En el contraste de homogeneidad de la Sección 7.3 hemos tomado tres muestras de 50 individuos cada una, en total 150 individuos. El estadístico de contraste ha valido \(X^2=28.59\), por lo que la magnitud del efecto en ese test ha sido de \(w=\sqrt{28.59/150}=0.4366\), entre mediana y grande según la función cohen.ES, que nos da las magnitudes del efecto que por convención se entienden como pequeñas, medianas o grandes para los diferentes tests considerados en el paquete pwr:
cohen.ES(test="chisq", size="medium")$effect.size## [1] 0.3cohen.ES(test="chisq", size="large")$effect.size## [1] 0.5¿De qué tamaño deberíamos haber tomado las muestras para garantizar una potencia del 90%, suponiendo que esperásemos una magnitud del efecto mediana y tomásemos un nivel de significación \(\alpha=0.05\)?
pwr.chisq.test(df=6, sig.level=0.05, w=0.3, power=0.9)##
## Chi squared power calculation
##
## w = 0.3
## N = 193.5425
## df = 6
## sig.level = 0.05
## power = 0.9
##
## NOTE: N is the number of observationsHubiéramos necesitado como mínimo un total de unos 194 individuos: si queríamos tomar las tres muestras del mismo tamaño, esto significa tres muestras de como mínimo 65 individuos cada una.
Para obtener solo el tamaño de la muestra, se puede añadir el sufijo $N:
pwr.chisq.test(df=6, sig.level=0.05, w=0.3, power=0.9)$N## [1] 193.54257.5 Guía rápida
tablecalcula tablas de contingencia de frecuencias absolutas.prop.tablecalcula tablas de contingencia de frecuencias relativas.addmarginssirve para añadir a unatableuna fila o una columna obtenidas aplicando una función a todas las columnas o a todas las filas de la tabla, respectivamente. Sus parámetros principales son:margin: igualado a1, se aplica la función por columnas, añadiendo una nueva fila; igualado a2, se aplica la función por filas, añadiendo una nueva columna; igualado ac(1,2), que es su valor por defecto, hace ambas cosas.FUN: la función que se aplica a las filas o columnas; su valor por defecto essum.
colSumscalcula un vector con las sumas de las columnas de una matriz o una tabla.rowSumscalcula un vector con las sumas de las filas de una matriz o una tabla.chisq.testsirve para realizar tests \(\chi^2\) de independencia y homogeneidad. El resultado es unalistformada, entre otros, por los objetos siguientes:statistic(el valor del estadístico \(X^2\)),parameter(los grados de libertad) yp.value(el p-valor). Sus parámetros principales en el contexto de esta lección son:simulate.p.value: igualado aTRUE, calcula el p-valor mediante simulaciones.B: en este último caso, permite especificar el número de simulaciones.
pwr.chisq.test, del paquete pwr, sirve para calcular uno de los parámetros siguientes a partir de los otros cuatro:N: el tamaño de la muestra.df: el número de grados de libertad del contraste.sig.level: el nivel de significación \(\alpha\).power: la potencia \(1-\beta\).w: la magnitud del efecto.
7.6 Ejercicios
Modelo de test
(1) Hemos observado dos variables cualitativas en una muestra de una población. Cada variable tiene 3 niveles. La tabla de contingencia resultante ha sido la siguiente: \[ \begin{array}{c|ccc} &X&Y&Z\cr\hline A & 2 & 17 & 11\cr B & 8 & 10 & 25\cr C & 3 & 14 & 5 \end{array} \] ¿Es verdad que, si estas variables aleatorias fueran independientes, las frecuencias esperadas de cada combinación de niveles, uno de cada variable, serían todas \(\geq 5\)? Tenéis que contestar SI, en mayúsculas y sin acento, o NO.
(2) Hemos observado dos variables cualitativas en una muestra de una población. Una variable tiene 4 niveles y la otra 3. La tabla de contingencia resultante ha sido la siguiente:
\[
\begin{array}{c|cccc} &A&B&C & D\cr\hline X & 50 & 19&17 & 21\cr Y &69 & 47 & 56 & 37 \cr Z &33 & 23 & 18 & 21 \end{array}
\]
Emplead la función chisq.test para contrastar si estas dos variables son independientes o no. Tenéis que dar el p-valor del test (redondeado a 3 cifras decimales, sin ceros innecesarios a la derecha) y decir (contestando SI o NO) si, con un nivel de significación \(\alpha=0.1\), podríamos rechazar la hipótesis nula de que estas dos variables son independientes. Dad las dos respuestas en este orden y separadas por un único espacio en blanco.
(3) Hemos realizado un test \(\chi^2\) de independencia sobre una muestra de 200 individuos, con un nivel de significación de 0.1. Las variables objeto de estudio tenían 5 y 6 niveles, respectivamente. El estadístico de contraste ha valido \(16.56\). ¿Cuál es el p-valor del contraste? ¿Cuál es la potencia del contraste realizado? Tenéis que dar ambos valores en este orden, redondeados a 3 cifras decimales sin ceros innecesarios a la derecha, y separados por un único espacio en blanco.
Ejercicios
(1) En un estudio para determinar si la tensión arterial depende del grupo sanguíneo, se tomó una muestra aleatoria de 1500 personas, se determinó su grupo sanguíneo y se les midió la tensión, clasificándola en baja, normal o alta. Los resultados fueron los de la tabla siguiente:
\[ \begin{array}{l} \hphantom{TensionAAa}\text{Grupo sanguíneo}\\ \begin{array}{l|cccc} \text{Tensión} & \text{A} & \text{B} & \text{AB} & 0\\ \hline \text{Baja} & 28 & 9 &7 & 31\\ \text{Normal} & 543 & 211 & 90 & 476\\ \text{Alta} & 44 &22 &8 & 31\\ \hline \end{array} \end{array} \]
Como hemos dicho, a partir de estos datos queremos decidir si hay asociación entre la tensión arterial y el grupo sanguíneo.
(a) Según el diseño del experimento, ¿qué tipo de contraste estamos realizando: de independencia o de homogeneidad?
(b) Comprobad si se dan las condiciones para poder efectuar un test \(\chi^2\). En caso afirmativo, llevadlo a cabo. En caso negativo, usad el método de Montecarlo.
Respuestas al test
(1) NO
Nosotros lo hemos mirado con
A=matrix(c(2,17,11,8,10,25,3,14,5),nrow=3,byrow=TRUE)
RS=rowSums(A)
CS=colSums(A)
frec.esperadas=RS%*%t(CS)/sum(A)
frec.esperadas## [,1] [,2] [,3]
## [1,] 4.105263 12.947368 12.947368
## [2,] 5.884211 18.557895 18.557895
## [3,] 3.010526 9.494737 9.494737(2) 0.128 NO
Nosotros lo hemos resuelto con
X=matrix(c(50,19,17,21,69,47,56,37,33,23,18,21),nrow=3,byrow=T)
round(chisq.test(X)$p.value,3)## [1] 0.128(3) 0.681 0.777
Nosotros lo hemos resuelto con
p.valor=1-pchisq(16.56,(5-1)*(6-1))
potencia=pwr.chisq.test(N=200, df=(5-1)*(6-1), sig.level=0.1, w=sqrt(16.56/200))$power
round(c(p.valor,potencia),3)## [1] 0.681 0.777Soluciones sucintas de los problemas
(a) De independencia, porque la muestra es transversal
(b)
Datos=rbind(c(28,9,7,31),c(543,211,90,476),c(44,22,8,31))
Datos| 28 | 9 | 7 | 31 |
| 543 | 211 | 90 | 476 |
| 44 | 22 | 8 | 31 |
n=sum(Datos)
n## [1] 1500freq.esp.tipo=colSums(Datos)
freq.esp.tensión=rowSums(Datos)
frec.esperadas=freq.esp.tensión%*%t(freq.esp.tipo)/n
frec.esperadas| 30.75 | 12.10 | 5.25 | 26.90 |
| 541.20 | 212.96 | 92.40 | 473.44 |
| 43.05 | 16.94 | 7.35 | 37.66 |
chisq.test(Datos)##
## Pearson's Chi-squared test
##
## data: Datos
## X-squared = 5.1163, df = 6, p-value = 0.529El valor por defecto de
margines el vector de todas las dimensiones de la tabla. Hay que recordar que, aunque ahora sólo tratamos con tablas bidimensionales, contablese pueden especificar tablas de contingencia de un número arbitrario de dimensiones.↩︎