DISTRIBUCIÓN JI CUADRADO, Apuntes de Psicología

¡Descarga DISTRIBUCIÓN JI CUADRADO y más Apuntes en PDF de Psicología solo en Docsity! ANALISIS DE DATOS EN PSICOLOGIA 1. La distribución Prácticamente todas las pruebas estadísticas se basan en la distribución de probabilidad de ciertas variables aleatorias. La razón para que esto sea así tiene que ver de que muchas variables y experimentospresenta parecerse a alguna distribución conocida número suficientemente grande de veces población de 100.000 individuos donde una determinada variable aleatoria una distribución normal (N( sería la forma que adoptaría el histograma realizado con todos los valores de la población? Crearemos con rnorm( ) > poblacion=rnorm(100000,100,15) > head(poblacion) # veamos los 6 primeros valores [1] 89.37865 98.44241 111.18166 110.39095 121.08057 91.67609 > summary(poblacion) Min. 1st Qu. Median Mean 3rd Qu. 29.10 89.84 99.98 99.96 110.10 163.30 El comando summary( )como ya sabes características de una variable.Si mínimo, máximo así como la media y los tres cuartiles de la variable.Podemos ver que los valores de la población de 100000 individuos generadas mediante encuentran entre 29.1 y 163.3. La media es prácticamente 100 y con los argumentos entregados (100 y 15 para la media y la desviación típica respectivamente).Recuerda que si replicas estas instrucciones tus valores serán algo diferentes. > library(MASS) #Multivariate Analysis in Social Sciences > truehist(poblacion) . El Análisis de la Varianza. Háztelo EN R. La distribución 27 χ2 n un histograma que tiende a cuando se realizan de forma repetida un . Veamos un ejemplo:si disponemos de una µ=100,σ=15)con media 100 y desviación típica 15 dicha población ficticia. # mu=100 sigma=15 Max. nos permite conocer de forma rápida las ésta es cuantitativa la función nos entrega los ín , por tanto χ2 con el hecho xsiguiese )¿cuál dices rnorm se , coherente , ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 28 La librería MASS nos permite mediante la función truehist( ) generar un histograma con mejor apariencia visual que el generado con hist( ).En ese histograma el mayor volumen de valores se encuentra entre 80 y 120 centrado en la media 100. Así mismo,es evidente el patrón característico de la distribución normal que produce una caída exponencial de la probabilidad de verificar un valor en los extremos de la curva, así como la igualdad de probabilidad de encontrar un valor tan extremo como 160 o 40. Ese histograma está compuesto por un número muy grande de barras cuyo perfil nos recuerda una campana. En élhemos proyectado unos200 intervalos entre los valores 29 y 163 (mínimo y máximo de nuestra población). Como sabemos, muchas variables en la naturaleza tienen esa distribución, pero lo realmente interesante es que hay una función f(x) que para valores de x entre - ∞ y +∞ presenta precisamente un perfil muy parecido. ( ) 2 1 2 2 1 2 x f x e µ σ πσ − −    = Es decir, para todo valor de x cuando le aplicamos esta función f(x), obtenemos un nuevo valor que es posible proyectar en una gráfica de infinitos puntos entre el mínimo y el máximo de los valores de x. A pesar de que R dispone de una función que nos da precisamente f(x) para cualquier valor (dnorm( )), podemos programar una a la que llamaremos dist.normal.fncque tenga x, mu y sigmacomo argumentos y f(x) como resultado. > dist.normal.fnc=function(x,mu,sigma) { fx=(1/sqrt(2*pi*sigma^2)) * exp(-0.5 * ( (x-mu)/sigma)^2) return(fx) } Observa que en el cuerpo de la función hemos replicado la ecuación de la distribución normal anterior. En R, si escribes pi en consola obtendrás el número π con seis decimales. Por otra parte, la función sqrt( ) nos devuelve la raíz cuadrada del número entregado como argumento. Si escribes en consola sqrt(100), obtendrás su raíz cuadrada 10. Por último, tenemos la función exp( ) que devuelve el resultado de elevar el número e (2,718282; base de los logaritmos naturales)al valor entregado como argumento.Veamos algunos ejemplos: > exp(1) [1] 2.718282 El exponencial de 1 es precisamente el número e, base de los logaritmos naturales porque 2,718 elevado a 1 es justamente el mismo valor. Ahora veremos que el exponencial de 2 es el resultado de elevar el número e a 2. > exp(2) [1] 7.389056 > 2.718282^2 # ^ es el símbolo que debemos usar para indicar “elevado a” [1] 7.389057 Veamos qué nos devuelve nuestra función (dist.normal.fnc ) para valores de xiguales a la media (100)y para puntuaciones tan extremas como 40 y 160 (Observa la gráfica anterior). ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 31 Dado que pnorm nos devuelve acumulados si queremos saber la probabilidad de verificar un valor mayor o igual a un valor dado (140) en una normal N(100,15) utilizaremos 1-pnorm. > 1-pnorm(140,100,15) [1] 0.003830381 Vamos a crear un data.frameque contenga75 valores en 2 variables aleatorias normales N(0,1) e independientes entre sí. > dos_variables=data.frame(z1=rnorm(75), z2=rnorm(75)) > dos_variables z1 z2 1 -0.640241005 0.949772281 2 -0.605080181 0.122949020 3 0.205633778 0.838481907 4 1.611343894 0.376440745 5 -0.256572474 -0.800790192 6 1.125630811 -0.707707316 7 -0.151457957 0.514694586 8 0.420414936 0.529758516 9 -0.115752298 0.135029506 10 -1.045094089 1.272266925 11 1.526055628 0.529717040 12 -2.100581597 -1.362282553 13 0.504094571 -1.035498466 14 -0.702547809 -0.078609846 15 -1.551576610 1.482090645 16 0.031089891 1.686167992 17 0.200068352 -0.005960142 18 0.498078090 -1.485345237 19 0.827876523 2.155310481 20 1.915948821 -0.474269820 21 -0.433034042 0.716459112 22 0.374536983 0.854207799 23 0.536532284 0.397706575 24 0.170282135 -0.374137626 25 -0.170479086 1.126974792 26 -1.347319765 -1.138682948 27 -0.501448203 0.480081399 28 0.126890408 0.481606155 29 1.011449703 1.964979006 30 0.245374901 0.994981030 31 -0.383765275 0.431260453 32 -0.808922706 -1.792226374 33 -0.353292627 -0.713855028 34 0.261855834 -0.458816393 35 0.270660978 -0.436194671 36 1.052240159 0.023039492 37 0.297722514 0.847457629 38 -0.666795036 -0.645646635 39 1.250119451 0.319064105 40 1.176559174 0.596246787 41 0.507264789 -0.012697805 42 0.158522273 -0.753936839 43 -0.492712710 -0.543223544 44 0.501809659 0.132020474 45 0.757784549 0.779599153 46 1.719091459 -0.688964546 47 -0.229805197 0.633040776 48 -0.222022092 0.065636855 49 0.582306383 -0.420986152 50 0.602491563 0.391791902 51 0.039485933 -1.978763687 52 -0.456459456 -0.002024598 53 0.036114871 -0.441314409 54 -0.182839654 -1.017308094 55 -0.486635659 0.793879675 56 -2.040744500 -0.512958438 57 0.512500581 -0.791240980 58 -0.364448093 -0.392036683 59 1.168076300 -1.173471217 60 -0.353265901 0.230136794 61 -1.155100382 0.153244734 62 -1.243772449 -1.292660796 63 0.062715191 0.441689122 64 -0.759394115 -1.599020380 65 -1.164941870 -0.246008959 66 0.798828772 0.164427141 67 -1.049685253 -0.318435232 68 -0.003586949 0.130505304 69 0.560274498 -1.381932470 70 0.676730793 -0.896303023 71 -1.677147054 1.020670931 72 0.117764007 -0.985139769 73 -0.608309122 0.602398317 74 -0.362590415 0.044328542 75 0.056791274 -0.483672851 Si observamos línea a línea cada una de los 75 valores de z1 y z2,veremos que es muy poco probable que en la misma línea se den simultáneamente dos valores grandes de z1y z2. De hecho, podemos preguntarle a R cuántos de los 75 valen en la misma línea más de 1.96 en valor absoluto para ambas variables. Esto lo haremos mediante el comando subset( ) que tiene la siguiente sintaxis: subset(objeto,condición a satisfacer). Veamos cómo se aplicaría a nuestro ejemplo. > subset(abs(dos_variables),(z1 >= 1.96 & z2 >= 1.96)) [1] z1 z2 <0 rows> (or 0-length row.names) R responde que ninguna línea de dos_variables supera nuestro criterio. Veamos qué ocurre si en vez de 75 valores lo hacemos sobre 1000. Volvemos a crear el objeto dos_variables pero esta vez con 1000 observaciones. ANALISIS DE DATOS EN PSICOLOGIA > dos_variables=data.frame(z1=rnorm(1000), z2=rnorm(1000)) > subset(abs(dos_variables),(z1 >= 1.96 z1 z2 295 2.209519 2.235239 983 2.299818 2.850362 Sólo un 2/1000 de los registros presentan valores simultáneos tan grandes para z2.Es decir, que cuando agrupamos poco probable verificar valores elevados simultáneamente. Qué ocurriría si para cada sujeto (fila de dos variables) generamos una nueva variable de nombre sumaque sea de hecho >dos_variables$suma=apply # como queremos que sea la suma del cuadrado de # hemos definido como función # solamente sum, R obviamente nos daría el sumatorio de cada # fila pero no la suma de La función apply(datos, función (ya existente o creada por nosotros) sus elementos de fila (marginal =1) o de columna ( queremos que se realice para cada fila que separa al objeto (dos_variables deseamoscalcular. Hagamos un histograma de la nueva variable creada > truehist(dos_variables$suma) Hay un patrón claro en este gráfico: una fuerte asimetría p la "probabilidad" de observar valores mayores de veremos que es prácticamente 2 (número de variables cuadráticas sumadas). Vamos a intentar profundizar en la información mejor repetiremos el experimento con 5000 valores presenten perfiles más suaves. > dos_variables=data.frame(z1=rnorm( . El Análisis de la Varianza. Háztelo EN R. La distribución 32 & z2 >= 1.96)) dos variables normales e independientes es muy (en valor absoluto) para ambas la suma del cuadrado de z1 y z2. (dos_variables,1,function(x)sum(x^2) cada valor, lo “al vuelo”. De haber puesto cuadrado de cada elemento de fila. que.marginal,funcion), permite apli a una base de datos o data.frame marginal=2). En nuestro caso , lo hemos indicado con un 1 después de la coma ) dela función de suma de cuadrados : ositiva y una 6.Si pedimos la media de que encierra este gráfico y para comprenderlo , de tal forma que los h 5000), z2=rnorm(5000)) χ2 z1 y ) car cualquier sobre , como que intensa caída de esta variable istogramas ANALISIS DE DATOS EN PSICOLOGIA > dos_variables$suma=apply > par(mfrow=c(1,2)) > truehist(dos_variables$z1) > truehist(dos_variables$z2) > par(mfrow=c(1,1)) # Re > truehist(dos_variables$suma) > abline(v=mean(dos_variables$suma)) El comando par( ), permite modificar muchos aspectos de la forma en que se realizan las gráficas. Entre ellas, con el subcomando que queremos dividir la ventana gráfica en irá llenando la pantalla fila a fila regrese a su aspecto original (sin división) Esta gráfica conjunta nos va a permitir entender la distribución del histograma de nuestra variable suma. En primer lugar, observamos z1 y z2 están en la vecindad del 0 y entre 1 y frecuentes a la media (cero). ellos mismos en valor absoluto correspondan a la barra del 0 en el histograma inferior segu solicitado una línea vertical en la media de (abline(v=mean(suma)) variables z al cuadrado sumadas cuadráticas grandesa la media 0, . El Análisis de la Varianza. Háztelo EN R. La distribución 33 (dos_variables,1,function(x)sum(x^2)) #Divide la pantalla en una fila y 2 columnas -establece la pantalla a 1 gráfico por ventana mfrow=c(1,2) le estamos indicando a R dos partes: 1 fila y 2 columnas. Cada gráfico , a medida que los solicitamos. Para que la ventana la restauramos a c(1,1). que los valores más frecuentes en -1. Es decir, 0, 1 y -1 son las distancias más Los cuadrados de esos números son precisamente iguales a . Esto explica que cuando se suman la ida de la barra del 1. sumapara el último histograma y vemos cómo coincide esta media con el número de (2).Fíjate en ese histograma que las cuando se suman dos variables tipificadas e χ2 s alturas mayores Se ha distancias ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 36 > points(20,0,col='brown',pch=19) > points(20,0.031,col='magenta',pch=19) En la página 30 veíamos las funciones para la distribución normal que nos permitían conocer la función de probabilidad dnorm, de densidad pnorm, etc. Para χ2, R sigue el mismo paralelismo sólo que al prefijo d,p,q y r le sigue chisq. dchisq(x,gl) función de probabilidad χ2 con gl, grados de libertad pchisq(x,gl) función de densidad χ2con gl, grados de libertad qchisq(p,gl) valor de χ2 que verifica una probabilidad acumulada rchisq(n,gl) Genera una distribución de n observaciones con una distribución χ2 con gl grados de libertad. Para demostrar nuestra afirmación de que el sumatorio de k variables tipificadas e independientes al cuadrado siguen una distribución χ2 con k grados de libertad, crearemos una función que precisamente haga esta tarea. crea.xi.fnc=function(n,gl) { z_dat=mvrnorm(n, mu=rep(0,gl),Sigma=diag(1,gl,gl)); z_dat=scale(z_dat); z2=z_dat^2; return(xi=apply(z2,1,sum)); } # mvrnorm( ), es una función de la librería MASS que nos permite crear # distribuciones normales multivariadas (muestras de n observaciones # en p variables con una matriz de correlaciones o covarianzas # poblacional Sigma). Como argumentos la función necesita el vector de # medias poblacional y la matriz Sigma de varianzas y # covarianzas o correlaciones poblacional de la que queremos # extraer muestras multivariadas Esta función es el equivalente a la que ya existe en R rchisq( ). Sin embargo, merece la pena crearla y describirla porque contiene algunas funciones interesantes. Vemos que lo primero que hacemos es asignar az_datel resultado de aplicar mvrnorm( ) con argumentos mu=rep(0,gl)(repite 0 gl veces) y Sigma=diag(1,gl,gl), una matriz diagonal, porque consideramos que las gl variables son independientes entre sí (ausencia de relación lineal entre ellas). > Sigma=diag(1,gl,gl) ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 37 > Sigma [,1] [,2] [,3] [,4] [,5] [1,] 1 0 0 0 0 [2,] 0 1 0 0 0 [3,] 0 0 1 0 0 [4,] 0 0 0 1 0 [5,] 0 0 0 0 1 Otra función muy importante es scale( ). Con ella provocamos la estandarización de las variables z generadas. Cuando solicitamos una muestra con distribución normal, es el valor poblacional de µel que vale realmente cero. En una muestra cualquiera, la media de zestará próxima a 0 y su desviación típica a 1.scale( ) centra la variable exactamente en 0 y desviación típicaigual a1. Llamemos a nuestra función con los siguientes argumentos: > xi_5=crea.xi.fnc(n=1000,gl=5) > summary(xi_5) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.1566 2.7310 4.2110 4.9950 6.6000 35.1800 # la media es casi exactamente 5 (las variables sumadas) > truehist(xi_5,xlab=’’) # xlab=’’ no pone etiqueta en x > x1=seq(0,35,0.001) > lines(x1,dchisq(x1,gl),lwd='2',type='l',col='blue') Con el comando lines anterior dibujamos la χ2 teórica con 5 grados de libertad sobre el histograma previamente graficado con truehist. Ahora generaremos 20 muestras aleatorias mediante un bucle for( )y proyectaremos en esa gráfica un histograma suavizado de cada una (gráfica derecha). for (i in 1: 20) { lines(density(rchisq(1000,gl)),lwd=1.5,lty=2,col='red') abline(v=mean(rchisq(1000,gl))) } Efectivamente, la suma del cuadrado de las puntuaciones tipificadas de k variables independientes genera una distribución empírica que sigue la χ2teórica con k grados de libertad. En ambas gráficas vemos que el valor más frecuente es 5, el mayor valor obtenido es 35, pero de carácter excepcional (1 entre 1000). Si utilizamos el criterio percentílicopara definir un valor anormalmente grande en la suma de 5 variables típicas cuadráticas e independientes diremos que: ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 38 > c(qchisq(0.95,5), qchisq(0.99,5), qchisq(0.999,5)) [1] 11.07050 15.08627 20.51501 Para 5 grados de libertad si utilizamos como criterio de distancia “anormalmente grande”a 0 el p95, p99 o p99.9 tenemos que dichos valores son 11, 15 y 20.5 respectivamente. Dibujemos nuevamente la misma distribución, pero ahora proyectaremos tres verticales relativas a los tres criterios percentílicos estimados. Asimismo, dibujaremos un punto en la posición 12.7 para valorar si una distancia cuadrática de 12.7 puntos es grande (mayor que cero) o, por el contrario, pequeña o común en una distribución χ2 con 5 grados de libertad. > x1=seq(0,35,0.001) > gl=5 > plot(x1,dchisq(x1,gl), lwd='2',type='l',col='blue') > text(15, 0.15, col='black',font=2,'La distribución ') > text(20.5, 0.153, col='black',font=2, cex=1.5,expression(chi^2)) > abline(v=gl); abline(h=0) > q95=qchisq(0.95,gl); p95=dchisq(q95,gl); > q99=qchisq(0.99,gl); p99=dchisq(q99,gl); > q999=qchisq(0.999,gl); p999=dchisq(q999,gl); > lines(c(q95,q95),c(0,p95),col='red',lwd=2) > lines(c(q99,q99),c(0,p99),col='red',lwd=2) > lines(c(q999,q999),c(0,p999),col='red',lwd=2) > points(12.7,0.005,pch=19,col='red',cex=2) El código anterior puede parecer complicado, pero es el que nos permite situar las 3 barras verticales rojas hasta la altura definida por la función en cada valor de x (11, 15 y 20.5). Con el comando points(x,y), situamos un punto en la coordenada deseada de la gráfica. Es importante resaltar que no podemos ver la línea vertical relativa al percentil 99.9 porque aunque presente en la gráfica su altura es muy pequeña como revela el código que presentamos a continuación. > dchisq(qchisq(0.999,5),5) [1] 0.0004336301 Fíjate que solicitamos el valor de la función χ2 (valor en el eje y) (dchisq( )) para aquel valor que deja por debajo de si al 99.9% de todos los valores en una distribución con 5 grados de libertad qchisq( ). Como puedes observar este valor de altura es demasiado pequeño como para que su representación sea visible en la gráfica. Compáralo con los que se corresponden con la primera y segunda línea roja vertical: > dchisq(qchisq(0.95,5),5) [1] 0.01932464 > dchisq(qchisq(0.99,5),5) [1] 0.004127806 ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 41 A la vista de estas diferenciasnos podemos preguntar. ¿Han cambiado los psicólogos sus preferencias vocacionales entre el año 2000 y 2004? Para responder a esta pregunta vamos a calcular la distancia cuadrática que hay entre lo observado y esperado aplicando la siguiente ecuación. ( )22 ˆ Observada Esperada Esperada χ  −  =     ∑ Para comprender el significado de estealgoritmo de cómputo de distancia observada- esperada (en adelante obs-esp) debemos preguntarnos qué resultado obtendríamos para cada especialidad si la frecuencia observada fuese en realidad exactamente la esperada. En ese caso el numerador del sumatorio anterior sería cero para cada especialidad, como también lo sería el cociente resultante. Por ello, el sumatorio a través de todas las especialidades sería también cero. En nuestro ejemplo, es obvio que dicha distancia no será cero. Sin embargo, sí que podemos plantearnos si es lo suficientemente distante de cero como para considerar que los psicólogos han cambiado realmente sus preferencias vocacionales entre los años 2000 y 2004. > xi=sum((observada-esperada)^2/esperada) > xi [1] 3.908945 Hemos obtenido una distancia cuadrática de 3.9 puntos. Para valorar si es una distancia a cero grande o pequeña, debemos compararla con la mayor distancia que cabe esperar en una distribución χ2 con grados de libertad igual al número de categorías-1. En nuestro ejemplo, son cuatro especialidades y,por tanto, 3 grados de libertad (4 - 1).¿Por qué 3 grados de libertad? Debemos plantearnos las cuatro categorías desde el siguiente punto de vista. ¿Cuántas de las cuatro podemos modificar su valor garantizando que el valor total sean 500 observaciones? Obviamente tres, porque si asignamos cualquier valor arbitrario a cada una de esas categorías, la cuarta deberá ajustarse a un valor fijo que garantice el sumatorio final de 500 observaciones. Luego nuestro vector de categorías profesionales tiene tres grados de libertad. Por ese motivo acudimos a una distribución χ2 con 3 grados de libertad para valorar si 3.9 es una distancia frecuente o infrecuente cuando 4 categorías fluctúan libremente (en realidad 4 menos 1) Dibujemos una χ2 teórica con 3 grados de libertad. > x=seq(0,20,0.01); gl=3 > plot(x,dchisq(x,gl), lwd='2',type='l',col='blue',axes=F) axis(1,at=seq(0,20,1)) axis(2,at=seq(0,0.26,0.02)) text(6, 0.20, col='black',font=2,'La distribución ') text(8.7, 0.202, col='black',font=2, cex=1.5,expression(chi^2)) text(10.5,0.20,col='black',font=2, cex=1,'con 3 gl.') > abline(v=3) # media de la χ2 con gl=3 > abline(h=0) > q95=qchisq(0.95,gl) > lines(c(q95,q95),c(0,0.08),col='red',lwd=2) # Línea roja vertical del p95 Situamos en la gráfica nuestra distancia cuadrática de 3.9 puntos. > points(3.9,0.01,pch=19,col='red',cex=2); box() ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 42 Vemos que esta distancia es un valor muy común para esta distribución χ2 (está próximo al valor más probable que es la media 3). De hecho, podemos plantearnos algunaspreguntas importantes:¿cuán probable es una distancia menor o igual a 3.9 puntos en esa distribución? ¿ y mayor o igual a 3.9 puntos?¿cuál sería la distancia que comenzaría a ser indicativo de distancia anormalmente grande o poco común? > pchisq(xi,3) # probabilidad de χ2 ≤ 3.9 en χ2(3) [1] 0.7285342 > 1-pchisq(xi,3) # probabilidad de χ2 ≥ 3.9 en χ2(3) [1] 0.2714658 > qchisq(0.95,3) # Valor de χ2 con [p(x) ≥ 7.81] ≤ 0.05 [1] 7.814728 Nuestro valor de 3.9 puntos deja por debajo de si el 72.9% (0,7285342) de todos los valores, lo cual nos indica que la probabilidad de verificar un valor mayor o igual a ese es del 0.271 (27,2%). Si definimos el criterio de distancia “anormalmente grande” como aquel valor que supere al 95% de todos los valores, tenemos que 7.8 sería aquelcon una probabilidad de ser verificado menor o igual al 5%. Luego nuestra decisión de si los psicólogos han cambiado sus preferencias vocacionales sería: dado que una distancia obs-esp mayor o igual a 3.9 puntos es muy frecuente (más frecuente que el criterio de excepcionalidad prefijado del 5%) cuando en la población la frecuencia observada es igual a la frecuencia esperada,diremos que los psicólogos no han cambiado sus preferencias vocacionales (p > 0.05). A lo largo de este capítulo vamos a necesitar realizar la bondad de ajuste a χ2 en múltiples ocasiones. Por este motivo, parece rentable crear una función propia que devuelva el valor del estadístico de distancia cuadrática obs-esp así como la probabilidad asociada (P(x ≥χ21-α)) cuando entregamos como argumentoslos vectores de frecuencias absolutas observadas y de probabilidades poblacionales. > bondad_ajuste_xi=function(obser,p) { esper=sum(obser)*p xi=sum((obser-esper)^2/esper) p=1-pchisq(xi,length(obser)-1) return(c(xi=xi, p=p)) } ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 43 Una de las fortalezas de R es que permite llevar a cabo simulaciones de datos mediante extracciones muestrales de poblaciones simuladas de forma muy sencilla y rápida. En estadística es muy común este tipo de trabajos de investigación donde se pretende estudiar mediante métodos intensivos de replicación los efectos que tiene la violación de los supuestos sobre las distintas pruebas estadísticas. En general, este tipo de estrategia recibe el nombre de Simulación de Monte Carlo. A lo largo de este manual haremos muchas porque es una herramienta excelente para aprender la estadística de modelos. En los siguientes párrafos vamos a explicar el proceso para extraer muestras aleatorias de una población de la que conocemos los parámetros a estudiar en las distintas muestras. De nuestra población de 20.000 psicólogos guardada en el objeto R población,vamos a extraer una muestra de 500 psicólogos y calcularemos, mediante la función creada anteriormente, el estadístico de distancia y suprobabilidad asociada. El comando sample(x,n) extrae muestras de tamaño n del total de observaciones x. Veamos un ejemplo que nos permita comprender el uso de dicha instrucción. Supongamos que tenemos un vector con 20 valores y queremos extraer aleatoriamente 5 de ellos. > sample(20,5) #De 20 elementos posibles extrae 5 [1] 12 7 18 3 9 > sort(sample(20,5)) # Extrae 5 de 20 posibles y ordénalos [1] 6 8 9 19 20 El comando sample( ) lo que hace es darme en este ejemplo 5 valores enteros y aleatorios entre 1 y 20. Si por ejemplo asigno esos 5 valores aleatorios a un objeto de nombre arbitrario indice, puedo extraer una muestra aleatoria de ese tamaño de una determinada población de nombre por ejemplo pobla simplemente indicando: pobla[indice,]. Con esta instrucción estaremos declarando que queremos los casos completos (todas las columnas) de pobla cuyos números de filas coincidan con los guardados en indice. Cada vez que ejecuto sample( ) se extrae una tabla de índices diferente.En el ejemplo anterior de los 20 elementos de un objeto cualquiera,habrían sido seleccionados los valores de las filas6, 8, 9, 19 y 20. Ya podemos aplicar esta función para seleccionar una muestra de 500 psicólogos de nuestra población. > indice=sort(sample(20000,500)) # Extrae 500 de 20000 y los pone en indice indice contiene el número de las 500 filas seleccionadas al azar y ordenadas de población que formaran parte de la muestra. > muestra=data.frame(poblacion[indice,]) #Extrae las filas de población que coincidan con índice y todas #sus columnas ([indice,]) > names(muestra)='ocupa' > observada=table(muestra$ocupa); > observada clinica educativa laboral otras 140 165 124 71 > esperada=500*p; ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 46 Ya podemos preguntarnos cuántas de las mil replicaciones han generado un valor de distancia superior al p95 de una χ2 con 3 grados de libertad. > sum(almacen_xi >= qchisq(0.95,3))/1000 [1] 0.046 Efectivamente, utilizando como criterio de distancia distinta de cero el percentil 95 de una teórica con 3 grados de libertad, decidimos en aproximadamente 5 de cada 100 que es una distancia grande (cuando en realidad no debe serlo, dado que las 1000 muestras han sido extraídas de la misma población). En esta simulación hemos demostrado algo muy importante. Si utilizamos el criterio percentílico para decidir mediante el contraste de hipótesis la pertenencia o no de una muestra a una misma población nos equivocaremos como máximo en un porcentaje previamente conocido y que llamaremos de ahora en adelante α (alfa). Hemos creado una función de nombre bondad_ajuste_xi para calcular la distancia cuadrática observa a esperada para un vector de frecuencias absolutas observadas. En realidad R ya dispone de una función para ello de nombre chisq.test. Esta función para el caso de bondad de ajuste requiere como argumentos el vector de frecuencias absolutas observadas y el de probabilidades poblacionales (imprescindible para el cálculo de las frecuencias esperadas). Veamos su uso con los datos de nuestro ejemplo: > observada = c(172,136,130,62) > p.pobla = c(0.306,0.289,0.264,0.141) >chisq.test(observada, p = p.pobla) Chi-squared test for given probabilities data: observada X-squared = 3.9146, df = 3, p-value = 0.2708 Vemos que efectivamente obtenemos exactamente el mismo resultado que con nuestra función bondad_ajuste_xi. Recordemos cuales eran los datos con los que iniciamos la búsqueda de la respuesta de si nuestra muestra de 500 psicólogos del año 2004 habían cambiado o no sus preferencias vocacionales desde el año 2000. > observada=c(172,136,130,62) > esperada=p.pobla * 500 > esperada clinica educativa laboral otras 153.000 144.500 132.000 70.475 Ahora vamos a hacer un nuevo experimento. Vamos a solicitar la distancia obs-esp para esta muestra, tal y como hicimos anteriormente. > # Muestra de 500 psicólogos >dist=bondad_ajuste_xi(observada,p.pobla); dist xi p 3.9089448 0.2714658 Si multiplicamos ambos vectores por 1.5, las distancias entre las especialidades seguirán siendo proporcionalmente iguales pero ahora tendremos una muestra de 750 psicólogos (500 * 1.5 =750) > # 750 psicólogos ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 47 > dist=bondad_ajuste_xi(observada*1.5,p.pobla); dist xi p 5.8634171 0.1184477 Por último repetiremos el experimento multiplicando ambos vectores por un factor de 2 (1000 psicólogos). > dist=bondad_ajuste_xi(observada*2,p.pobla); dist xi p 7.8178895 0.0499292 Ha ocurrido algo sorprendente y es que con la muestra de 500 la distancia obs-esp era común (pequeña, frecuente) ya que la probabilidad de la misma era mayor o igual que el criterio de excepcionalidad prefijado del 5%. En el segundo experimento con 750,nuestra conclusión es idéntica a pesar de que ahora la probabilidad es menor. Sin embargo, con 1000 psicólogos, y proporcionalmente la misma distancia obs-esp, concluimos que sí han cambiado sus preferencias vocacionales (p < 0.05). ¿Cómo es posible que la misma distancia nos conduzca a cambiar la decisión dependiendo del tamaño de la muestra de la que surge el estadístico? En el próximocapítulo veremos que a partir de un cierto tamaño de muestra, si éste es suficientemente grande, nuestra decisión mediante el contraste de hipótesis tenderá a ser que es una distancia grande a 0, lo cual hace inservible el instrumento de bondad de ajuste para tomar dicha decisión. Comprenderemos además que para que esta herramienta sea realmente útil, debemos contar con un número adecuado de sujetos muestrales. Si tenemos demasiados, siempre diremos “distancia grande y rara”, pero si tenemos menos de los necesario, diremos justo lo contrario “distancia pequeña y, por tanto, frecuente”. 1.2. La prueba de independencia de variables mediante χ2 En la vida cotidiana la mayoría de las situaciones que anticipamos y vivimos, pueden concluir con múltiples resultados diferentes. El concepto de lo incierto es inherente a la vida, y la incertidumbre, como dificultad para predecir con cierto margen de seguridad lo que sucederá, está incorporada a nuestra cotidianeidad. ¿Subirán los tipos de interés este año? ¿Cuántos toxicómanos abandonarán el programa de deshabituación antes de finalizarlo completamente? ¿Si tienen un fuerte vínculo familiar y afectivo disminuirá el riesgo de abandono prematuro? ¿Cuántas asignaturas aprobaré dado que dedico 3 horas diarias al estudio? Es evidente que estudiar es una condición necesaria para aprobar, pero no garantiza la ausencia de suspensos. En el proceso de superar un curso académico intervienen otras variables que pueden alterar el devenir deseable de los acontecimientos. Si como terapeutas recibimos la visita de un paciente con fobia específica a volar,esperamos que su recuperación sea más rápida y eficaz que si padeciese un trastorno de tipo agorafóbico. El concepto de lo incierto tiene además un ámbito de aplicación de vital importancia en ciencia. En ciencias humanas y de la salud, los tratamientos deben demostrar su eficacia a través de la evidencia empírica de que el o los grupos tratados presentan un nivel de mejora superior a lo que cabría esperar por puro azar. En todos los ejemplos anteriores, si bien no tenemos certeza de qué ocurrirá, manejamos un criterio de certidumbre que nos ayuda a predecir el resultado más “probable”.En este sentido, sabemos que cuanto más estudiemos más difícil haremos el resultado de suspender y viceversa. No fumar no garantiza la inmunidad contra el cáncer pero todas las ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 48 investigaciones demuestran que una manera exitosa de incrementar la probabilidad de contraerlo es adquirir el hábito de fumar unos 20 cigarrillos diarios. En ciencia es muy frecuente el interés por conocer si existe una asociación o dependencia estadística entre dos características de una determinada población. Así, podemos estar interesados en conocer si el nivel de ingresos y el nivel educativo son estadísticamente independientes o interdependientes. Para entender el proceder de una prueba de independencia vamos a partir de un sencillo ejemplo de una población formada por 20.000 personas, 50% mujeres y 50% de hombres, donde la mitad de cada género fuma y la otra mitad no. > smoke=data.frame(sexo=rep(c('mujer','hombre'),each=10000), fuma=rep(rep(c('si','no'),each=5000),2)) > head(smoke) sexo fuma 1 mujer si 2 mujer si 3 mujer si 4 mujer si 5 mujer si 6 mujer si > tabla=table(smoke); tabla; fuma sexo no si hombre 5000 5000 mujer 5000 5000 Vamos a completar esta tabla calculando sus marginales de fila y de columna. Es decir, cuántos hombres y mujeres y cuántas personas fuman y no fuman.Cada celdilla de intersección género x conducta de fumar ha sido dividida en 2 partes. La izquierda está formada por el valor real del número de hombres y mujeres que fuman o no en esa población. En la parte derecha vamos a situar el valor que le correspondería a esa celdilla si las 20.000 personas se repartieran en la intersección (4 celdillas) de forma azarosa (género y conducta de fumar son absolutamente independientes). La probabilidad de caer en cada celdilla deberá entonces estar exclusivamentecondicionada por el peso de sus marginales (número de hombres o mujeres, y cuantas personas fuman y no fuman) en relación al total de observaciones de la tabla. ¿Cuántos hombres de las 20.000 personas esperamos que fumen? Esta pregunta incluye a los 10.000 hombres y a las 10.000 personas que fuman (marginales de fila y de columna de la celdilla ser hombre y fumar). Si multiplicamos ambos marginales 10.000 x 10.000 y al resultado lo dividimos por el total de observaciones de la tabla (20.000) obtenemos un valor de 5.000 hombres esperados en la conducta de fumar si ambas variables son absolutamente independientes entre sí. Dado que los marginales de fila y de columna son idénticos para el resto de las celdillas, concluiremos que el número de personas esperadas por azar en cada una de ellas es el mismo (5.000). Fuma No fuma Marginal Columna Hombre 5000 ? 5000 ? 10000 Hombres Mujer 5000 ? 5000 ? 10000 Mujeres Marginal Fila 10000 Fuman 10000 No Fuman 20000 ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 51 para el histograma de esas distancias? ¿cuál será la distancia observada-esperada más frecuente? > alma_xi=c() #Con c() generamos un almacén inicialmente vacio > for (i in 1:10000){ # Extraeremos 10000 muestras de la población indice=sort(sample(20000,100)) # cada muestra de tamaño 100 muestra=smoke[indice,] tabla=table(muestra$sexo,muestra$fuma); tabla=with(muestra,table(sexo,fuma)) # Otra forma de pedir la tabla xi=chisq.test(tabla) alma_xi[i]=xi$statistic } > x=seq(0,10,0.01) > truehist(alma_xi) > rug(alma_xi,col='blue') # añade marcas verticales de color azul > lines(x,dchisq(x,1),col='red',lwd=2) El promedio de los 1000 valores está próximo a 1. Los valores mayoritarios de distancia están en la vecindad del 0 y por puro azar ha habido algunos valores de distancia muy grandes en concreto 1 de casi 15 y algunos un poco superiores a 10. Este histograma es absolutamente coherente con lo esperado dado que hemos extraído muestras de la misma población. > qchisq(0.95,1) # P95 para una χ2 de 1 grado de libertad [1] 3.841459 Cuántas distancias de las 1000 han sido superiores a este criterio del p95? > sum(alma_xi >= 3.84)/10000 [1] 0.039 Sólo un 3.9% de las distancias nos invitarían a concluir que ambas variables son interdependientes a pesar de que esas muestras también han sido extraídas de la misma población donde el género y la conducta de fumar son realmente independientes. Veamos ahora un ejemplo diferente y con valores observados más plausibles y ecológicos. En la tabla adjunta tenemos las frecuencias observadas o empíricas ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 52 registradas de 818 toxicómanos que asisten a un programa de deshabituación de la heroína. Transcurrido un año de la finalización del mismo, 637 recaen mientras que181 no lo hacen. Tal y como podemos observar la mitad había conseguido un trabajo estable y la otra mitad no.En este sentido, nos planteamos si la variable cualitativa recaer en la heroína (con dos posible valores si y no) es independiente de la variable tener trabajo estable con dos alternativas: si y no. > heroina=data.frame(no_recaida=c(76,105),recaida=c(335,302)) > row.names(heroina)=c('No_trabajo','Si_trabajo') > heroina no_recaida recaida No_trabajo 76 335 Si_trabajo 105 302 > sum(heroina) # número total de toxicómanos [1] 818 >m_fila=apply(heroina,1,sum) # marginales de fila Si_trabajo No_trabajo 411 407 >m_columna=apply(heroina,2,sum) # marginales de columna no_recaida recaida 181 637 ¿Es independiente recaer en la heroína de tener trabajo estable? > contraste=chisq.test(heroina) > contraste Pearson's Chi-squared test with Yates' continuity correction data: heroina X-squared = 5.9196, df = 1, p-value = 0.01497 > contraste$expected no_recaida recaida No_trabajo 90.94 320.1 Si_trabajo 90.06 316.9 A partir de este resultado dado que se cumplen los supuestos tenemos que concluir que la distancia observada-esperada es distinta de cero. Es decir, los valores observados por celdilla son diferentes a los que cabría esperar si recaer en la heroína y tener trabajo estable fuesen dos variables independientes entre sí. Luego hemos de concluir que ambas variables son interdependientes (p < 0.05).Para comprender el razonamiento que nos conduce a esta decisión tenemos nuevamente que valorarlo desde el punto de vista de la decisión en el contraste por probabilidad asociada al estadístico. Hemos obtenido una distancia cuadrática de 5.9196 puntos. La probabilidad de verificar un valor mayor o igual a este en una distribución χ2 con 1 grado de libertad cuando las variables de fila y de columna son independientes entre si es de 0.0149. Como esa probabilidad es inferior al criterio α menos exigente (5%) diremos que recaer en la heroína y tener o no trabajo estable son dos sucesos interdependientes (p < 0.05). Una vez que tomamos la decisión de interdependencia, es imprescindible describir la naturaleza de la misma. Para ello debemos valorar los porcentajes de fila y columna de ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 53 la tabla de contigencia contrastada. Con ese objetivo, hemos creado la función porcentaje.tabla.fnc que requiere como único argumento una tabla de doble entrada (filas y columnas) y devuelve como resultado múltiples tablas con los porcentajes marginales de fila, columna, así como los de la intersección de ambas. porcentaje.tabla.fnc=function(tabla){ tabla=as.matrix(tabla) por.fila=prop.table(tabla,1) # 1 marginal de fila por.colu=prop.table(tabla,2) # 2 marginal de columna por.tabla=prop.table(tabla) # sin marginal es el % de la tabla tabla.c= cbind(tabla,margin.table(tabla,1)) tabla.c= rbind(tabla.c, cbind(t(margin.table(tabla,2)),margin.table(tabla))) return(list(tabla.original=tabla.c,por.fila=por.fila, por.columna=por.colu, por.tabla=por.tabla)) } porcentaje.tabla.fnc(tabla) #EJECUTA LA FUNCION HASTA AQUI > porcentaje.tabla.fnc(tabla) $tabla.original no_recaida recaida No_trabajo 76 335 411 Si_trabajo 105 302 407 181 637 818 $por.fila no_recaida recaida No_trabajo 0.1849148 0.8150852 Si_trabajo 0.2579853 0.7420147 $por.columna no_recaida recaida No_trabajo 0.4198895 0.5259027 Si_trabajo 0.5801105 0.4740973 $por.tabla no_recaida recaida No_trabajo 0.09290954 0.4095355 Si_trabajo 0.12836186 0.3691932 Los porcentajes marginales de fila nos indican que hay igualdad de proporción entre los que tienen trabajo estable y los que no. Sin embargo, vemos que en lo que a la recaída se refiere, recae casi un 78% de los toxicómanos frente al 22% que permanecen “limpios”.En los porcentajes de fila para la intersección (interior de la tabla) vemos que de los que no trabajan recaen el 81% frente al 74% de los que si tienen trabajo estable.La siguiente tabla de porcentajes de columna para la intersección nos informa de que del total de los que recaen el 52% no tiene trabajo estable frente al 47% de los que si tienen. Podemos concluir, por tanto, que el tener trabajo establetiene un efecto protector sobre el riesgo de recaer en la heroína después de haber asistido a un programa de desintoxicación. ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 56 Prueba con los siguientes contrastes en el simulador: simula.p.asociada(gl=8,valor=15.51,0.95) simula.p.asociada(gl=8,valor=15.51,0.99) simula.p.asociada(gl=8,valor=27,0.999) Intenta después de algunos ensayos creados por ti adivinar la decisión que tomará el simulador a partir del valor, la confianza y los grados de libertad simulados. ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 57 Ejercicios 1.- En la página 40 de este capítulo nos planteábamos si los psicólogos habían cambiado sus preferencias vocacionales entre los años 2000 y 2004 a partir de una muestra obtenida en este último año de 500 psicólogos. Otra pregunta que puede ser respondida con la misma técnica podría ser: ¿Es igual la probabilidad para las cuatro especialidades? Obviamente solo debemos modificar el vector de probabilidades poblacionales (probabilidad esperada) dividiendo la probabilidad total por el número de especialidades (0.25 para cada especialidad) . > observada.año.2004=c(172,136,130,62) > p.pobla=rep(0.25,4) > p.pobla [1] 0.25 0.25 0.25 0.25 > chisq.test(observada.año.2004, p= p.pobla) Chi-squared test for given probabilities data: observada.año.2004 X-squared = 50.592, df = 3, p-value = 5.976e-11 La distancia cuadrática obtenida de 50.6 puntos nos indica que para 3 grados de libertad la probabilidad de obtener una distancia como esa cuando Ho es cierta es prácticamente cero (6-10). Por ello rechazamos Ho (p < 0.001) y declaramos que las cuatro especialidades no son equiprobables en la población de psicólogos de la que ha sido extraída la muestra de 500. 2.- Loftus y Palmer (1974). Journal of Verbal and Learning & Verbal Behaviour, 13, 585-589.Se presentó a 150 personas una secuencia de 20 segundos de duración de un accidente leve de tráfico. Una semana más tarde a 50 individuos se les preguntó si habían visto cristales rotos en el suelo después de que los vehículos se “estrellaran”. A los siguientes 50 se les realizó la misma pregunta cambiando el verbo por “golpearan”. A los últimos 50 simplemente se les preguntaba si recordaban haber visto cristales rotos en el suelo. Evidentemente no había cristales en el suelo en la película que se les presentaba. ¿Es independiente la respuesta afirmativa de la forma en que se pregunta? > accidente=data.frame( Estrellarse=c(16,34),Golpearse=c(7,43), Control=c(6,44)) > row.names(accidente)=c('Si','No') > accidente Estrellarse Golpearse Control Si 16 7 6 No 34 43 44 Respuesta Condición verbal Estrellarse Golpearse Control Si 16 7 6 No 34 43 44 ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 58 > chi_accidente=chisq.test(accidente); chi_accidente Pearson's Chi-squared test data: accidente X-squared = 7.78, df = 2, p-value = 0.02045 > chi_accidente$expected Estrellarse Golpearse Control Si 9.667 9.667 9.667 No 40.333 40.333 40.333 #Se cumplen los supuestos Hemos obtenido para esta tabla una distancia cuadrática observada-esperada de 7.78 puntos que para dos grados de libertad tiene una probabilidad asociada de 0.02. Es decir, que esa distancia ocurre como máximo un 2% de las veces que Ho es correcta y por lo tanto independientes entre sí. Dado lo escaso de su ocurrencia, concluiremos que ambas variables son interdependientes y que la respuesta no es independiente de la forma en que se pregunta (p < 0.05).Si solicitamos las tablas de porcentajes de fila y de columna, podremos dimensionar adecuadamente el alcance de dicha interdependencia. Es imprescindible ejecutar primero la función porcentaje.tabla.fnc si no se ha hecho antes (página 53) > porcentaje.tabla.fnc(accidente) $tabla.original Estrellarse Golpearse Control Si 16 7 6 29 No 34 43 44 121 50 50 50 150 $por.fila Estrellarse Golpearse Control Si 0.5517241 0.2413793 0.2068966 No 0.2809917 0.3553719 0.3636364 $por.columna Estrellarse Golpearse Control Si 0.32 0.14 0.12 No 0.68 0.86 0.88 $por.tabla Estrellarse Golpearse Control Si 0.1066667 0.04666667 0.0400000 No 0.2266667 0.28666667 0.2933333 En la tabla de porcentajes de fila podemos ver que del total que han respondido SI el 55.1 % pertenece a la condición Golpearse y Control que es menos de la mitad (24 y 20% respectivamente. En la tabla de porcentajes de columna, vemos que los sujetos de la condición Estrellarse han declarado haber visto cristales rotos en un 32% porcentaje que duplica el observado para las condiciones Golpearse y Control (4.6 y 4% respectivamente). Resulta claro desde estos datos que la interdependencia encontrada entre la respuesta a la pregunta y la condición experimental de pertenencia tiene que ver con que estar en la condición Estrellarse duplica la probabilidad de falso recuerdo frente a las otras dos condiciones experimentales. ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 61 Laboratorio de Estadística 1.A LA DISTRIBUCION XI CUADRADO OBJETIVO GENERAL Entender por qué si una distribución empírica no sigue una teórica se producirá un incremento de los valores máximos esperados para esa distribución. *************************************************************************** Objetivos específicos. 1. Saber generar una distribución normal multivariada con media y matriz de correlaciones dada (pag. 36). 2. Saber utilizar la función matrix( ) para crear una matriz de datos a partir de un vector de entrada con un número de filas y columnas dado (pag. 18 y 19). 3. Cambiar los nombres de las variables en un data.frame mediante c( )(pag. 23 y en esta práctica). 4. Aplicar la función apply( ) con función al vuelo, para conseguir el sumatorio cuadrado de los elementos de la fila de un data.frame (pag. 32). 5. Solicitar un histograma de una variable continua y dibujar encima su distribución teórica f(x) para demostrar el ajuste de la primera a la segunda (pag. 29). 6. Solicitar percentiles de una variable continua mediante el uso combinado de quantile( ) y c( ) (pag. 7 y 34). 7. Solicitar los valores de superficie (acumulados) para una distribución chi cuadrado con gl grados de libertad (pag. 42). 8. Aprender a utilizar la función de conectores de fila rbind( ) para añadir filas a una tabla (en esta práctica). 9. Saber interrogar una base de datos (data.frame) para conocer el numero de valores que superan una condición (pag. 10 y 31). *************************************************************************** PLANTEAMIENTO DEL PROBLEMA Demuestra mediante una pequeña simulación de Monte Carlo que si dos variables tipificadas Z no son independientes entre sí (cor(z1,z2)≠ 0). La variable generada a partir de ellas con dos grados de libertad no seguirá una distribución χ2centrada con gl= 2. *************************************************************************** CONSEJO: Cada vez que veas un comando o función nueva. Escríbela en consola con una interrogación delante. Por ejemplo ?head. Eso producirá que R te presente la ayuda para que veas los argumentos y características de la función. Por otra parte al final de la ayuda vienen ejemplos de uso de esa función que al copiarlos y pegarlos en consola te permite comprender y mejorar tu dominio de las funciones de R. Al principio puede parecer complicado pero poco a poco te será más fácil entenderlo. Sobre todo con las funciones que mas utilizas. ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 62 # Cargaremos en primer lugar la librería MASS que nos permite mediante la función mvrnorm( ) generar muestras de tamaño n en p variables con media y matriz de correlaciones deseada así como usar histogramas de mayor calidad mediante truehist( ). > library(MASS) # Creo la matriz de correlaciones poblacional que quiero que tengan los datos que voy a generar para dos variables z tipificadas. En este caso quiero que las dos variables z no sean independientes (r = 0.85). La llamaré sigma. > sigma=matrix(c(1,0.85,0.85,1),nrow=2,ncol=2); sigma # Compruebo que la matriz de correlaciones tiene la estructura deseada. Normalmente significa comprobar que en la diagonal sólo hay 1, que la submatriz inferior a la diagonal es igual a la superior y que el valor de la correlación es aproximadamente el simulado (0.85). En nuestro caso al ser una 2 x 2 es muy fácil verificar que es correcta. # genero el data.frame de nombre "data" a partir de la función mvrnorm( ) que significa crea una distribución empírica normal multivariada. Los argumentos de mvrnorm son n, vector de medias y matriz de correlaciones que quiero que tengan mis datos. Como son dos variables con la misma media 0 introduzco como mu=c(0,0) y sigma que es la matriz de correlaciones anteriormente generada. > data=data.frame(mvrnorm(n=5000,mu=c(0,0),sigma)) > head(data) # Veo los 6 primeros valores de data # Cambiamos el nombre de las variables de data a z1 y z2 para que sea más coherente con nuestro objetivo de manipular variables z y no directas (X) > names(data)=c('z1','z2') > head(data) # Chequeamos que la matriz de la muestra "data" generada tiene en realidad la correlación esperada entre las variables z mediante el uso de la función cor( ). > cor(data) # Ahora vamos a crear nuestra variable chi cuadrado empírica con dos grados de libertad, sumando para cada fila de data el cuadrado de z1 y z2. Esto lo haremos mediante apply( ), que como ya sabes tiene como argumentos: # apply(objeto de datos, marginal, función a aplicar). Como queremos el cuadrado de las puntuaciones de z1 y z2 debemos definirla como función al vuelo: function(x) sum(x^2). Fíjate que no es lo mismo que decir sum(x)^2. En este último caso estarías planteando el cuadrado de la suma total por fila, mientras que lo que deseamos es el cuadrado de cada valor en la fila. # Fijate que el 1 a la derecha de data, representa el marginal de fila > data$chi2=apply(data,1,function(x) sum(x^2)) # Pedimos el histograma de nuestra variable cuadrática generada. ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 63 # Como es truehist( ) y no hist( ) necesitamos tener cargada la librería MASS mediante el comando library(librería) > truehist(data$chi2) # recuerda objeto$variable # Fíjate en la gráfica y veras que los valores posibles de nuestra variable chi2 estan entre0 y 25. Por ese motivo si queremos superponer una chi cuadrado teórica con dos grados de libertad. Es importante que generemos una variable de nombre x (es un nombre arbitrario) que reproduzca la secuencia entre 0 y 30 con saltos muy pequeños (0.001) para que los puntos de la función proyectada parezcan en realidad una línea. La función seq(inicio,fin,incremento) nos permite crear esa secuencia. > x=seq(0,25,0.001) # pide en consola ?seq # Con lines(que en eje x?, que en eje y?) dibujamos el valor de la función teórica chi cuadrado con dos grados de libertad en el eje y para todo valordel vector x creado anteriormente. Esto nos permitirá valorar gráficamente si el histograma de nuestra variable chi2 empírica sigue en realidad la distribución teórica chi cuadrado con dos grados de libertad (gl=2) esperada. > lines(x,dchisq(x,2)) #recuerday repasa la función dchisq(x,gl) de la página 36 del manual. # Aparentemente el histograma se ajusta bien a la teórica. Pero sabemos que nuestras variables z no son independientes entre sí. Luego si z1 tiene por azar un valor grande. Al estar muy correlacionada con z2 ser muy probable que z2 también lo tenga. Eso debería conducir a valores cuadráticos grandes más frecuentes que cuando z1 y z2 son en realidad independientes (esta idea se explica ampliamente en la página 33 del manual, repásalo si no lo recuerdas). # Vamos a demostrarlo mediante la comparación de los percentiles de nuestra chi2 generada con los acumulados para una chi teórica con 2 grados de libertad. # Si nuestro análisis es correcto deberíamos esperar que la proporción de valores grandes sea siempre mayor en la variable chi2 generada por nosotros que en la teórica esperada bajo el supuesto de independencia de las variables implicadas. > PCentiles=quantile(data$chi2, c(0.95,0.99,0.999)); PCentiles # Aqui solicitamos los percentiles 95,99 y 99.9 de la variable chi2 generada. Observa que debo poner dentro de un compose c( ), la secuencia de percentiles que quiero. # Ahora solicitamos los valores de una chi cuadrado gl=2 (teórica) que dejan por debajo de si al 95,99 y 99.9% de todos los valores. Los llamaremos teoricos. > teoricos= qchisq(c(0.95,0.99,0.999),df=2); teoricos # df es degrees of freedom (grados de libertad) # Hay una conclusión evidente que podremos ver mejor si unimos ambos resultados en una misma tabla. Lo haremos mediante la función rbind( ). Que significa conecta filas (r de row, que es fila en ingles). Lo que hace ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 66 # cambiaremos los nombres de las dos variables por sus relativos en castellano. > names(cars)= c('velocidad','distancia') # Con los comandos combinados cut( ) y quantile( ) (Pag. 13) pedimos ver la variable velocidad recodificada en cuatro niveles percentílicos. Si deseas ver los puntos de corte que utilizará cut para recodificar la variable ejecuta: >quantile(cars$velocidad) >with(cars, cut(velocidad,quantile(velocidad))) # Observa que las dos primeras observaciones son NA (no disponibles). Esto ocurre porque por defecto la instrucción cut no incluye los valores mínimos. Dado que queremos que los incluya debemos modificar el subcomando que actúa por defecto (include.lowest = F). > with(cars, cut(velocidad,quantile(velocidad), include.lowest = T)) # Ahora si tenemos todos los casos incluidos en algunas de las cuatro categorías de recodificación. Una vez comprobado que la instrucción es correcta, crearemos la nueva variable en el data.frame. >cars$velocidad.r= with(cars, cut(velocidad,quantile(velocidad),include.lowest = T)) # Haremos lo mismo para la variable distancia. >cars$distancia.r= with(cars, cut(distancia,quantile(distancia),include.lowest = T)) # Es el momento de comprobar que las dos variables creadas presentan la estructura deseada. > summary(cars) # Ahora crearemos la tabla de doble entrada para el cruce de ambas variables recodificadas con el comando table( ). > tabla=with(cars, table(velocidad.r,distancia.r)); tabla # Es el momento de plantearnos el contraste de independencia con la función chisq.test( ). >contraste=chisq.test(tabla); contraste # Aparentemente debemos concluir que ambas variables son interdependientes (p < 0.001). Sin embargo no debemos valorar el contraste hasta que no comprobemos el cumplimiento de los dos supuestos de la técnica. Por ello solicitaremos los valores esperados que están guardados en el objeto contraste donde hemos situado el resultado de la ejecución de chis.test. Si tenemos dudas acerca del nombre que debemos utilizar simplemente solicitamos el nombre de las variables contenidas en contraste. ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 67 > names(contraste) # Efectivamente es expected el nombre que debemos usar y como ya sabes precedido del nombre del objeto contenedor y el carácter $ >contraste$expected # Vemos que se incumple claramente el supuesto que declara que como máximo el 20% de las celdillas serán inferior a 5 (lo son el 100%). Esto hace que no podamos utilizar esta prueba de independencia con las variables así recodificadas. Vamos a ejecutar de nuevo la recodificación y el contraste pero esta vez lo haremos creando sólo dos niveles (inferior o superior al percentil 50) para ambas variables. # Como la instrucción quedará más grande y por tanto más confusa crearemos en primer lugar los objetos cortes.velocidad y cortes.distancia. >cortes.velocidad= quantile(cars$velocidad,c(0,0.50,0.999)) >cortes.distancia= quantile(cars$distancia,c(0,0.50,0.999)) # Se crearan dos niveles para cada variable: primer nivel entre 0 y el percentil 50 y el segundo nivel entre el percentil 50 y el 99.9 > cars$velocidad.r= with(cars, cut(velocidad,cortes.velocidad,include.lowest = T)) > cars$distancia.r= with(cars, cut(distancia,cortes.distancia,include.lowest = T)) > tabla=with(cars, table(velocidad.r,distancia.r)); tabla # Por último llevamos a cabo el contraste y la comprobación de los supuestos. > contraste=chisq.test(tabla); contraste > contraste$expected # Dado que se cumplen los dos supuestos concluimos que velocidad y distancia (en su versión recodificada) son dos variables interdependientes (p < 0.001). # Si quisiéramos responder a la pregunta por punto crítico simplemente debemos conocer cuál es el valor que en una chi cuadrado para 1 grado de libertad deja por debajo de si al 95, 99 y 99.9 de los valores. >qchisq(c(0.95,0.99,0.999), 2) # Como el valor de distancia obtenido en el contraste de 18.6 es superior a 13.81 (pc 99.9) diremos por punto crítico que ambas variables son interdependientes (p < 0.001). # Debemos ahora cargar la función porcentaje.tabla.fnc( ) para poder valorar en qué consiste la interdependencia encontrada en la tabla. Nos trasladamos a la página 53 del manual y ejecutamos la función. > porcentaje.tabla.fnc(tabla) # Si observas la tabla que incluye los porcentajes de fila y de columna veras que esa interdependencia se refiere al hecho obvio de que a mayor ANALISIS DE DATOS EN PSICOLOGIA. El Análisis de la Varianza. Háztelo EN R. La distribución χ2 68 velocidad mayor distancia recorrida. # Prueba a solicitar la correlación de las dos variables contínuas velocidad y distancia. > cor(cars$velocidad, cars$distancia) # Alternativamente puedes utilizar esta otra instrucción > with(cars, cor(velocidad,distancia)) # O incluso esta otra que te dará la matriz de correlaciones completa. > cor(cars[,c(1,2)]) # Estamos pidiendo la correlación de las dos primeras columnas del objeto cars a través de todas las filas, que son precisamente las variables continuas velocidad y distancia. En el tema 3 veras que dicha interdependencia entre velocidad y distancia puede ser planteada de otra forma, para concluir algo muy similar aunque con un abordaje estadístico diferente. *************************************************************************** # EJERCICIO PROPUESTO 1.B *************************************************************************** En las pasadas elecciones de un determinado ayuntamiento de Canarias los distintos grupos políticos (CC,PP,PSOE,AFV,IUC,OTROS) obtuvieron los siguientes porcentajes de los votos 51.09,27.28,6.75,5.93,2.72,6.23 respectivamente. Ahora nos han encargado una encuesta de tendencia de voto en ese municipio y hemos preguntado a 509 encuestados no sólo por qué partido o sigla votaría en las próximas elecciones sino que también lo hemos hecho por el recuerdo de voto de las pasadas. Las frecuencias absolutas de recuerdo de voto obtenidas para cada partido fueron las siguientes: 248,145,36,20,20,40. ¿Podemos considerar que la proporción de recuerdo de voto de los encuestados es igual al voto realmente emitido en las pasadas elecciones?. En realidad podemos formular la pregunta de otra manera: ¿se ajusta la distribución de voto por partido registrada en la muestra de 509 sujetos a el voto realmente observado en las pasadas elecciones en esa misma población de la que se extrajo dicha muestra? Si esto es así la variable tendencia de voto investigada tendría credibilidad y por tanto fiabilidad demoscópica. Si por el contrario demostramos que el recuerdo y el voto real ocurrido difieren, habría que modificar la tendencia de voto observada mediante métodos de ponderaciones de la prospección de tendencia por recuerdo, como única alternativa razonable a la evidencia de que voto ocurrido y su recuerdo difieren. Luego la muestra no parecería ser fiable. # Pista: Observa que los valores de referencia de voto real imprescindibles para el cálculo de p están en porcentajes y no en probabilidades que es como deben estar.