#####################################################################
###
###                         ATENCIÓN:
###
###  Los rectángulos grises como este contienen código que puedes 
###  copiar y pegar directamente en RStudio. Pero te recomendamos 
###  que, sobre todo al principio, te esfuerces en teclear los
###  comandos para que te resulte más fácil recordarlos.
###
###  Debajo de los rectángulos grises aparecerán las respuestas de R, 
###  en rectángulos blancos o en forma de figuras. No es necesario,
###  por supuesto, que copies esas respuestas a R, basta con ejecutar  
###  el código de los rectángulos grises.
###
#####################################################################

Práctica 1. ESTADÍSTICA DESCRIPTIVA

Objetivos

  • Resumir, ordenar y analizar conjuntos de datos.
  • Calcular diversas características de una variable estadística cuantitativa o categórica.
  • Representar gráficamente la distribución de frecuencias.
  • Realizar análisis exploratorio de datos.

Cuestión 1.

Se está estudiando la efectividad de dos tipos de antitérmicos, A y B, en la reducción de la temperatura corporal. El fichero adjunto:

contiene una tabla con datos sobre el tipo de antitérmico que se le suministró a cada uno de treinta niños (segunda columna) y la reducción de temperatura (primera columna), en grados centígrados, que se consiguió:

Apartado a.

Asegúrate de que has seleccionado el directorio de trabajo adecuado para esta sesión con R. Usa el menú Session-> Set Working Directory de RStudio o el comando setwd de R para hacer esto.

Indicación: Si tienes dudas de cómo hacer esto mira la Sección 3, pág. 12, del Tutorial02.

Apartado b.

Descarga el fichero y guárdalo en la carpeta datos de tu directorio de trabajo. Si no existe esa carpeta, créala antes.

Indicación: Para que todo el código que aparece más abajo funcione la carpeta de trabajo no es la carpeta datos, sino la carpeta que contiene a la carpeta datos.

Apartado c.

Lee ese fichero con R, usando el comando read.table y guárdalo en un data.frame con el nombre practica1c01. Asegúrate de usar las opciones header=TRUE, as.is=TRUE,sep="," al leer el fichero.

Solución:

# Usamos este comando para leer el fichero. 
practica1c01 = read.table("./datos/practica01c01.csv", header = TRUE, as.is=TRUE, sep=",")

Apartado d.

Comprueba que la lectura ha sido correcta usando el comando head para ver las primeras 6 líneas de la tabla de datos.

Solución:

# Así vemos el comienzo de la tabla.
head(practica1c01)
##   temperatura tipo
## 1         1.4    A
## 2         2.3    A
## 3         1.0    B
## 4         1.8    A
## 5         1.7    B
## 6         1.2    B

Apartado e.

Para acceder a los datos en las columnas del data.frame que has creado en el paso anterior puedes utilizar la notación practica01$temperatura y practica01$tipo. Compruébalo y usa esa notación para crear dos variables temperatura (de tipo numérico) y tipo (de tipo carácter) para cada una de las columnas de esa tabla (para no tener que usar esa notación $ en los siguientes apartados).

Solución:

# Creamos las variables tipo y temperatura
(tipo = practica1c01$tipo)
##  [1] "A" "A" "B" "A" "B" "B" "B" "B" "A" "A" "B" "B" "B" "A" "A" "A" "B"
## [18] "B" "A" "B" "A" "A" "B" "B" "A" "B" "B" "B" "B" "A"
(temperatura = practica1c01$temperatura)
##  [1] 1.4 2.3 1.0 1.8 1.7 1.2 1.4 1.6 2.1 1.6 1.7 1.6 1.5 1.4 1.9 1.8 1.1
## [18] 1.5 1.1 1.7 1.6 2.1 1.3 1.2 1.8 1.5 1.6 1.7 1.6 1.5

Apartado f.

La conversión de grados centígrados (Celsius) a grados Fahrenheit se hace con la fórmula \[T (^0F) = 1.8\cdot T(^0C)+ 32\] Cree la variable temperaturaF, que contenga las reducciones de temperatura de los treinta niños, en grados Fahrenheit, con una sola cifra decimal. Usa para eso la función round.

Solución:

# Conversión de Celsius a Fahrenheit
(temperaturaF = round(1.8 * temperatura + 32, digits = 1))
##  [1] 34.5 36.1 33.8 35.2 35.1 34.2 34.5 34.9 35.8 34.9 35.1 34.9 34.7 34.5
## [15] 35.4 35.2 34.0 34.7 34.0 35.1 34.9 35.8 34.3 34.2 35.2 34.7 34.9 35.1
## [29] 34.9 34.7

¿Cuál es la diferencia si en lugar de round usas signif, también con una cifra decimal?
¡Cuidado! No sobrescribas la variable temperaturaF al hacer esto.

Solución:

# Comparación de round con signif
signif(1.8 * temperatura + 32, digits = 1)
##  [1] 30 40 30 40 40 30 30 30 40 30 40 30 30 30 40 40 30 30 30 40 30 40 30
## [24] 30 40 30 30 40 30 30

Apartado g.

Cuál es la reducción de temperatura mínima que no ha sido superada por:

  • la cuarta parte de los niños.

Solución:

# En R los percentiles se calculan con quantile, indicando el pocentaje como tanto con uno mediate probs.
quantile(temperatura,probs = 0.25)
## 25% 
## 1.4
  • la mitad de los niños.

Solución:

# Podemos usar quantile como antes:
quantile(temperatura,probs = 0.5)
## 50% 
## 1.6
# o en el caso de la mediana:
median(temperatura)
## [1] 1.6
  • el 70% de los niños.

Solución:

# Como las anteriores:
quantile(temperatura,probs = 0.7, type=1)
## 70% 
## 1.7

Apartado h.

Determine el valor del coeficiente de asimetría

Solución:

# Para calcular esto se debe instalar una librería adicional en R.
if(!require(e1071)){install.packages("e1071", dependencies = TRUE)}
## Loading required package: e1071
# Ahora cargamos la librería
library(e1071)

# y usamos la función skewness de esa librería para calcular el coef. de asimetría.
(gamma1 = skewness(temperatura))
## [1] 0.2229468

¿Puede considerarse que la distribución de los datos de reducciones de temperaturas presenta asimetría a la izquierda?

Apartado i.

Divida el conjunto de datos en 5 intervalos de clase, entre 0.8 y 2.5. Solución:

# La forma general de agrupar por intervalos en R es con la función `cut`. Para usarla tenemos que definir los límites entre clases. Así que primero usamos `seq` para fabricar esos límites.
(limsClases = seq(0.8, 2.5, length.out = 6))
## [1] 0.80 1.14 1.48 1.82 2.16 2.50
# Ahora ya podemos usar cut así:
(temperaturaClases = cut(temperatura, breaks = limsClases, 
                         include.lowest = TRUE))
##  [1] (1.14,1.48] (2.16,2.5]  [0.8,1.14]  (1.48,1.82] (1.48,1.82]
##  [6] (1.14,1.48] (1.14,1.48] (1.48,1.82] (1.82,2.16] (1.48,1.82]
## [11] (1.48,1.82] (1.48,1.82] (1.48,1.82] (1.14,1.48] (1.82,2.16]
## [16] (1.48,1.82] [0.8,1.14]  (1.48,1.82] [0.8,1.14]  (1.48,1.82]
## [21] (1.48,1.82] (1.82,2.16] (1.14,1.48] (1.14,1.48] (1.48,1.82]
## [26] (1.48,1.82] (1.48,1.82] (1.48,1.82] (1.48,1.82] (1.48,1.82]
## Levels: [0.8,1.14] (1.14,1.48] (1.48,1.82] (1.82,2.16] (2.16,2.5]
# Como ves el resultado es una clasificación de los valores: para cada valor obtenemos el intervalo de clase al que pertenece. Podemos, por ejemplo, hacer una tabla de frecuencias de esas clases, usando:
table(temperaturaClases)
## temperaturaClases
##  [0.8,1.14] (1.14,1.48] (1.48,1.82] (1.82,2.16]  (2.16,2.5] 
##           3           6          17           3           1

Obtenga el histograma y decida a partir de su forma si la distribución es simétrica.

Solución:

# Al crear el histograma vamos a usar los límites de clase que hemos creado antes:
hist(temperatura, breaks = limsClases)

# Las alturas de estas columnas son las frecuencias de la tabla que acabamos de ver. 

# Añado la función de densidad, como complemento al histograma. Para eso hay que ajustar la escala vertical del histograma, usando 
hist(temperatura, breaks = limsClases, freq = FALSE)
lines(density(temperatura), col="red")

# Como comentario adicional, en conjuntos con tan pocos datos el histograma es muy sensible al número de clases y a la posición de las marcas de clase. 
# Prueba a ejecutar el comando
# hist(tempeartura, breaks = 6)
# cambiando el valor 6 por otros valores cercanos (por ejemplo de 5 a 10) para ver cómo cambia el histograma.

Apartado j.

Indique el valor de la marca de clase correspondiente al tercer intervalo de clase.

Solución:

# Esa marca de clase es la media entre los límites del intervalo.
mean(limsClases[3:4])
## [1] 1.65

¿Cuál es el porcentaje de datos que pertenecen a dicho intervalo?

Solución:

# La forma más fácil es hacer una tabla de frecuencias relativas de las clases:
table(temperaturaClases)/length(temperaturaClases)
## temperaturaClases
##  [0.8,1.14] (1.14,1.48] (1.48,1.82] (1.82,2.16]  (2.16,2.5] 
##  0.10000000  0.20000000  0.56666667  0.10000000  0.03333333

Apartado k.

Construya un diagrama de cajas para los datos.

Solución:

# Aparte de usar la función boxplot voy a asignar el resultado a una variable para poder acceder después a los componentes de la gráfica:
bxpltTemperatura =  boxplot(temperatura)

Apartado l.

¿Existen datos puntuales que puedan calificarse como atípicos? En caso afirmativo, indíquelos.

Solución:

# Usando la variable que hemos creado esto es muy sencillo:
bxpltTemperatura$out
## [1] 2.3

Apartado m.

Exclusivamente para las reducciones de temperatura debidas al antitérmico A, halle los valores de los parámetros estadísticos siguientes:

  • media:

Solución:

# Empezamos seleccionando los valores de la temperatura para los que el tipo es A.
(antiaA = temperatura[tipo == "A"])
##  [1] 1.4 2.3 1.8 2.1 1.6 1.4 1.9 1.8 1.1 1.6 2.1 1.8 1.5
# Y ahora calculamos su media:
mean(antiaA)
## [1] 1.723077
  • mediana:

Solución:

# Ahora obtenemos la mediana:
median(antiaA)
## [1] 1.8
  • moda:

Solución:

# Moda:
max(table(antiaA))
## [1] 3

Apartado n.

Para la variable tipo obtenga la tabla de frecuencias e indique qué porcentaje de niños han tomado el analgésico A.

Solución:

# La tabla de frecuencias se obtiene con:
table(tipo)
## tipo
##  A  B 
## 13 17
# Y si queremos convertirla en frecuencias relativas:
prop.table(table(tipo))
## tipo
##         A         B 
## 0.4333333 0.5666667

Apartado o.

Construya el diagrama de sectores correspondiente.

Solución:

# Recuerda que este tipo de gráficos **se desaconsejan** en general. Pero si no te queda más remedio, puedes usar:
pie(table(tipo), col = terrain.colors(3))

# Una alternativa mejor:
barplot(table(tipo), col = terrain.colors(3))

Cuestión 2.

Los datos siguientes corresponden a la presión arterial, en mm de Hg, de 20 participantes en un ensayo clínico.

\[ 90,\, 106,\, 117,\, 102,\, 109,\, 120,\, 115,\, 113,\, 113,\, 101,\, 99,\, 95,\, 126,\, 108,\, 107,\, 109,\, 100,\, 107,\, 106,\, 112 \]

Cree con ellos la variable presion.

Solución:

# Creamos el vector con los valores de la presión usando c (y copia/pega desde el enunciado de la práctica). 
presion = c(90, 106, 117, 102, 109, 120, 115, 113, 113, 101, 99, 95, 126, 108, 107, 109, 100, 107, 106, 112)

Apartado a.

¿Qué valor máximo de la presión arterial es superado por el 30% de las presiones?

Solución:

# El percentil 70 es el valor máximo de la presión arterial que es superado por el 30% de las presiones.
quantile(presion, 0.7)
##   70% 
## 112.3

Apartado b.

¿Qué valor mínimo de la presión arterial no es superado por el 15% de las presiones?

Solución:

# El percentil 15 es el valor mínimo de la presión que no es superado por el 15% de las presiones.
quantile(presion, 0.15)
##   15% 
## 99.85

Apartado c.

Obtenga el diagrama de tallo y hojas.

Solución:

# Recuerda que este tipo de diagramas **se desaconsejan**. Es mejor usar en su lugar diagramas de columnas, histogramas y boxplots.
stem(presion)
## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##    9 | 059
##   10 | 0126677899
##   11 | 23357
##   12 | 06
# Mejor usar:
hist(presion)

# Luego veremos una versión mejor del histograma. También sirven estos otros:
plot(density(presion))

boxplot(presion)

Apartado d.

¿Sugiere la forma del diagrama simetría en los datos?

Solución:
Ni el diagrama de densidad ni el boxplot indican asimetría.

Apartado e.

Obtenga el histograma correspondiente a una agrupación de los datos en seis intervalos de clase, con límite inferior 90 y límite superior 128.

Solución:

# Usamos `seq` para crear los límites de clase:
hist(presion, breaks = seq(90, 128, length.out = 7), col="orange", probability = TRUE)

# Y añadimos la densidad al mismo gráfico con lines (plot crearía un nuevo gráfico).
lines(density(presion),col="red", lwd=5)  

Práctica 2. REGRESIÓN LINEAL.

Objetivos:

  • Utilizar el modelo lineal para relacionar dos variables, obteniendo su ecuación, el grado de ajuste de la misma con los datos y su representación gráfica.
  • Estudiar otros modelos alternativos al lineal cuando éste no es adecuado y determinar cuál es el que mejor se ajusta a los datos.
  • Predecir valores de una variable en función de la otra y obtener residuos atípicos y puntos influyentes.

Cuestión 1.

Apartado a.

Los pesos (en kg) y estaturas (en cm) de una muestra de 10 estudiantes universitarios son:

\(\quad\)
Peso 82 75 70 68 44 63 80 70 54 54
Altura 185 185 180 178 159 170 190 172 162 165
\(\quad\)

Carga estos datos desde el siguiente fichero adjunto:

Recuerda guardarlo en la subcarpeta datos de tu directorio de trabajo. El resultado de la lectura del fichero deber ser un data.frame llamado practica02c01. Crea dos variables llamadas peso y altura que contengan, respectivamente, la primera y segunda columna de ese data.frame (para no tener que usar $).

Atención: aunque no se pida explícitamente, en los próximos ejercicios tendrás que repetir estas operaciones de cargar los datos y crear unas variables con las que trabajar cómodamente.

Solución:

# ¡¡Recuerda fijar el directorio de trabajo antes de ejecutar el código!!
# Acuérdate también de empezar explorando el fichero para ver si hay cabecera, cuál es el separador, etc.

# Ahora leemos el fichero con read.table.
practica02c01 = read.table(file = "./datos/practica02c01.csv", header = TRUE, sep=",")

# Usamos head para comprobar que ha ido bien:
head(practica02c01)
##   peso altura
## 1   82    185
## 2   75    185
## 3   70    180
## 4   68    178
## 5   44    159
## 6   63    170
# Y creamos las dos variables a partir de las columnas del data.frame.
peso = practica02c01$peso
altura = practica02c01$altura

Apartado b.

Dibuja el diagrama de dispersión correspondiente a los datos del apartado a (con peso en el eje x, altura en el eje y).

Solución:

# Basta con usar la función plot. El peso va en el eje horizontal x y la altura en el eje vertical y.
plot(peso, altura)

# El resultado como ves es un diagrama muy básico, pero a medida que aprendas más sobre plot lo podrás ir enriqueciendo tanto como quieras (ver Tutorial10).

Apartado c.

Obtén las rectas de regresión de la altura frente al peso y del peso frente a la altura. Dibuja ambas rectas en el diagrama de dispersion ¿Son la misma recta?

Solución:

# Empezamos pintando los puntos de la muestra, usando cruces como símbolos y aumentando su tamaño con cex.
plot(peso, altura, pch="+", cex=2)

# Ahora usamos la función lm para obtener la recta de regresión lineal de la altura frente al peso.
# Fíjate en la notación de fórmula que se usa en R.
(lmAP = lm(altura ~ peso))
## 
## Call:
## lm(formula = altura ~ peso)
## 
## Coefficients:
## (Intercept)         peso  
##    119.9911       0.8274
# Y ahora podemos añadir la recta de regresión a la gráfica con abline. Usamos color azul. Con lwd se controla la anchura del trazo.
abline(lmAP, col="blue", lwd=2)

# La ordenada en el origen y la pendiente de este modelo están en la componente coefficients del modelo lineal:
lmAP$coefficients
## (Intercept)        peso 
## 119.9911111   0.8274074
# Así que podemos usar esto para guardarlos en unas variables a1 y b1 (añadimos el 1 porque luego hay otra recta):
(a1 = lmAP$coefficients[1])
## (Intercept) 
##    119.9911
(b1 = lmAP$coefficients[2])
##      peso 
## 0.8274074
# También puedes usar estos valores para dibujar la recta así (el dibujo no cambiará, porque la recta es la misma):
abline(a = a1, b = b1, col="blue", lwd=2)

# Repetimos lo anterior cambiando los papeles de peso y altura. Fíjate en que cambiamos el nombre del modelo. 
(lmPA = lm(peso ~ altura))
## 
## Call:
## lm(formula = peso ~ altura)
## 
## Coefficients:
## (Intercept)       altura  
##    -125.881        1.099
# Vamos a llamar a2 y b2 a los coeficientes de este segundo modelo para distinguirlos del primero que hicimos.
(a2 = lmPA$coefficients[1])
## (Intercept) 
##   -125.8813
(b2 = lmPA$coefficients[2])
##   altura 
## 1.098977

Ahora tenemos dos rectas de regresión. La primera es: \[\text{altura} = a1 + b1 \cdot \text{peso}\] y la segunda es: \[\text{peso} = a2 + b2 \cdot \text{altura}\] Para ver si coinciden vamos a representarlas juntas, en la misma gráfica. Pero si ponemos la altura en el eje vertical y elpeso en el eje horizontal, como hemos hecho antes, la ecuación de la segunda recta no está en la forma adecuada. Tenemos que despejar la altura. Se obtiene: \[\text{altura} = -\dfrac{a2}{b2} + \dfrac{1}{b2}\cdot \text{altura}\] Así que debemos usar estos valores \[-\dfrac{a2}{b2},\qquad \dfrac{1}{b2}\] de ordenada en el origen y pendiente respectivamente para pintar la segunda recta en los mismos ejes que la primera.

# Por eso usamos las siguientes expresiones. La segunda recta se dibuja en color rojo.

plot(peso, altura, pch="+", cex=2) # Los puntos de la muestra.
abline(lmAP, col="blue", lwd=2) # Primera recta.
abline(-a2/b2, 1/b2, col="red", lwd=2) # Segunda recta. 

Apartado d.

¿Cuánto vale el coeficiente de correlación R para las variables peso y altura? ¿Qué valoración haces?

Solución:

# Usamos la función cor
cor(peso, altura)
## [1] 0.953573
#  Puedes comprobar que el orden no afecta al resultado:
cor(altura, peso)
## [1] 0.953573

Apartado e.

¿Cuál es la altura prevista para un estudiante que pesa 63kg? ¿Y cuál es el peso estimado para uno que tiene una estatura de 175cm?

Solución:

# Hay dos forma de predecir valores. Directamente, susituyendo en la ecuación de la recta de regresión.
peso0 = 63
(altura0 = lmAP$coefficients[1] + lmAP$coefficients[2] * peso0)
## (Intercept) 
##    172.1178
# O esta segunda forma usando predict, que puede parecer inicialmente más complicada,
# pero que es mucho más potente. Naturalmente, las respuestas coinciden.
predict(lmAP, newdata = data.frame(peso=63))
##        1 
## 172.1178

Cuestión 2.

Se realiza un estudio para establecer una ecuación mediante la cual se pueda utilizar la concentración de estrona en saliva, para predecir la concentración del esteroide en plasma libre. Se extrajeron, en pg/ml, los siguientes datos de 14 varones sanos:

\(\quad\)
saliva 7.4 7.5 8.5 9.0 9.0 11.0 13.0 14.0 14.5 16.0 17.0 18.0 20.0 23.0
plasma 30.0 25.0 31.5 27.5 39.5 38.0 43.0 49.0 55.0 48.5 51.0 64.5 63.0 68.0
\(\quad\)

Estos datos están disponibles en el siguiente fichero adjunto.

Solución:

# Lectura de datos y creación de variables:
practica02c02 = read.table("./datos/practica02c02.csv", header = TRUE, sep=",")

(saliva = practica02c02$saliva)
##  [1]  7.4  7.5  8.5  9.0  9.0 11.0 13.0 14.0 14.5 16.0 17.0 18.0 20.0 23.0
(plasma = practica02c02$plasma)
##  [1] 30.0 25.0 31.5 27.5 39.5 38.0 43.0 49.0 55.0 48.5 51.0 64.5 63.0 68.0

Apartado a.

Dibuja el correspondiente diagrama de dispersión para las variables:

  • concentración de estrona en saliva (eje x).
  • concentración del esteroide en plasma libre (eje y)

Solución:

# De nuevo, usamos plot. Vamos a añadir rótulos a los ejes:
plot(saliva, plasma, 
     xlab = "Concentración de estrona en saliva (pg/ml).", 
     ylab = "Concentración del esteroide en plasma (pg/ml).")

Apartado b.

Obtén la correspondiente recta de regresión para esas variables.

Solución:

# Recuerda, usamos lm para construir la recta y luego plot junto con abline para dibujarla:
(lm01 = lm(plasma ~ saliva))
## 
## Call:
## lm(formula = plasma ~ saliva)
## 
## Coefficients:
## (Intercept)       saliva  
##       8.645        2.727
plot(saliva, plasma)
abline(lm01, lwd=2, col="blue")

Apartado c.

Calcula el coeficiente de correlación. ¿Qué opinas?

Solución:

# Coeficiente de correlación:
cor(saliva, plasma)
## [1] 0.9544798
# El coeficiente de correlación no detecta ningún error evidente en el ajuste de la recta. Pero recuerda 
# que NUNCA debes usar ese valor como el único criterio para concluir que el ajuste es adecuado.

Apartado d.

¿Cuál sería la estimación de la concentración del esteroide en plasma libre que correspondería a una concentración de estrona en saliva de 17.7 pg/ml?

Solución:

# Otro  caso en el que podemos optar por sustituir en la recta
saliva0 = 17.7
lm01$coefficients[1] + lm01$coefficients[2] * saliva0
## (Intercept) 
##    56.91922
# o usar predict.
predict(lm01, newdata = data.frame(saliva=17.7))
##        1 
## 56.91922

¿Cuál sería la estimación de la concentración del esteroide en plasma libre que correspondería a una concentración de estrona en saliva de 26.3 pg/ml? ¿Cabría hacer alguna objeción a estas estimaciones?

Solución:
Es un caso claro de extrapolación, porque el valor excede el recorrido de valores de la concentración de estrona en saliva de la muestra. El valor que predice el modelo no es por lo tanto fiable.

Apartado e.

Obtén los valores predichos por el modelo para los puntos de la muestra.

Solución:

# El resultado de la función lm incluye los valores predichos para los puntos de la muestra bajo el nombre fitted.values:  
lm01$fitted.values
##        1        2        3        4        5        6        7        8 
## 28.82737 29.10011 31.82747 33.19116 33.19116 38.64588 44.10061 46.82797 
##        9       10       11       12       13       14 
## 48.19166 52.28270 55.01007 57.73743 63.19216 71.37425
# Vamos a ilustrar gráficamente la relación entre esos valores y los puntos originales de la muestra:
plot(saliva, plasma, col="black", pch="*", cex=2)
points(saliva, lm01$fitted.values, col="red", pch="*", cex=2)
abline(lm01, lwd=2, col="blue")
segments(saliva, plasma, saliva, lm01$fitted.values, lty=2)

Cuestión 3.

Se han tomado cinco muestras de la misma cantidad de glucógeno y se les ha aplicado una cantidad de glucogenasa, X, (en milimoles/litro) anotando en cada caso la velocidad de reacción, Y, (en micromoles/minuto). Se han obtenido los siguientes datos:

\(\quad\)
glucógeno 0.2 0.5 1.0 2.0 3.0
glucogenasa 8 10 18 35 60
\(\quad\)

Estos datos están disponibles en el siguiente fichero adjunto.

Solución:

# Lectura de datos y creación de variables:
practica02c03 = read.table("./datos/practica02c03.csv", header = TRUE, sep=",")

(glucogeno = practica02c03$glucogeno)
## [1] 0.2 0.5 1.0 2.0 3.0
(glucogenasa = practica02c03$glucogenasa)
## [1]  8 10 18 35 60

Apartado a.

Dibuja el diagrama de dispersión para esas variables: glucógeno en el eje x, glucogenasa en el eje y.

Solución:

# Como en los anteriores:
plot(glucogeno, glucogenasa)

Apartado b.

Obtén la recta de regresión para esas variables.

Solución:

# Usamos lm, plot, abline:
(lm01 = lm(glucogenasa ~ glucogeno))
## 
## Call:
## lm(formula = glucogenasa ~ glucogeno)
## 
## Coefficients:
## (Intercept)    glucogeno  
##       1.211       18.648
plot(glucogeno, glucogenasa)
abline(lm01)

Apartado c.

Si a una de las muestras le hubiésemos aplicado una concentración de 2.5 milimoles/litro de glucogenasa, ¿cuál habrá sido la velocidad de reacción? ¿Es fiable esta predicción?

Solución:
De nuevo, es extrapolación. Así que no, no es fiable y ni siquiera la haremos.

Cuestión 4.

Se sabe, por experiencia, que el aumento de peso de los embriones de pollo al transcurrir el tiempo sigue una ley de tipo exponencial. En un experimento se obtuvieron los pesos, en gramos, de un embrión desde el sexto día de su nacimiento hasta el decimosexto que aparecen a continuación:

\(\quad\)
dia 6 7 8 9 10 11 12 13 14 15 16
peso 0.029 0.052 0.079 0.125 0.181 0.261 0.425 0.738 1.130 1.882 2.812
\(\quad\)

Estos datos están disponibles en el siguiente fichero adjunto.

Solución:

# Lectura de datos y creación de variables:
practica02c04 = read.table("./datos/practica02c04.csv", header = TRUE, sep=",")

(dia = practica02c04$dia)
##  [1]  6  7  8  9 10 11 12 13 14 15 16
(peso = practica02c04$peso)
##  [1] 0.029 0.052 0.079 0.125 0.181 0.261 0.425 0.738 1.130 1.882 2.812

Apartado a.

Dibuja el diagrama de dispersión para esas variables: día en el eje x, peso en el eje y.

Solución:

# Fíjate en el aspecto del diagrama de dispersión. ¡Eso no es una recta!
plot(dia, peso)

# Aunque no se pide, vamos a calcular la correlación:
cor(dia, peso)
## [1] 0.8626562
# Fíjate en que no es especialmente alta, aunque la relación entre las dos variables es más que evidente.

Apartado b.

Construye un modelo exponencial para estos datos. Para ello, usa una ecuación de la forma: \[\text{peso} = a\cdot e^{b\cdot\text{ dia}}\] donde \(a\) y \(b\) son constantes. Indicación: toma logaritmos en esta ecuación.

Solución:
Tomando logaritmos (siempre neperianos, salvo aviso en contra) y usando sus propiedades: \[\log(\text{peso}) = \log a + \log(e^{b\cdot \text{ dia}}) =\log a + b\cdot\text{ dia}\] Hacemos un cambio de notación: \[A =\log a,\qquad \text{lpeso }= \log(\text{peso})\] Y entonces se tiene: \[\text{lpeso} = A + b \cdot \text{ dia}\] Esto es la ecuación de una recta en la que la variable \(x\) es dia, la variable \(y\) es lpeso, la ordenada en el origen es \(A\) y la pendiente es \(b\). Así que podemos usar R para ajustar una recta a los valores de \(\text{lpeso}\) frente a los de dia.

# Definimos la variable lpeso (recuerda que en R log es el logaritmo neperiano)
lpeso = log(peso)
  
# Esa recta se obtiene en R así:
(lm01 = lm(lpeso ~ dia))
## 
## Call:
## lm(formula = lpeso ~ dia)
## 
## Coefficients:
## (Intercept)          dia  
##      -6.192        0.451
# Vamos a representar el resultado en las variables originales. Empezamos con la muestra: 
plot(dia, peso)

# Ahora vamos a necesitar los coeficientes de la recta, que son:
A = lm01$coefficients[1] # ordenada en el origen
b  = lm01$coefficients[2] # pendiente

# Si miras la ecuación que hemos obtenido tomando logaritmos verás que:
a = exp(A)
# mientras que el valor de b no ha cambiado al tomar logaritmos. 

# Para representar la función exponencial vamos a darle valores en una sucesión de 30 puntos (nodos)
# que cubren todo el recorrido de dia en la muestra. 
diaNodes = seq(min(dia), max(dia), length.out = 30)

# Ahora usamos points para dibujar la exponencial. Al usar type="l" conseguimos que se pinte una curva 
# en lugar de puntos aislados.
points(diaNodes, a * exp(b * diaNodes ), type="l", lwd=2, col="blue")

# Como ves, el ajuste de la curva exponencial a los datos parece muy bueno.

Apartado c.

Calcula el coeficiente de correlación \(R\) y también \(R^2\). Interpreta sus valores.

Solución:

# Ya hemos visto antes la correlación de las variables originales.
cor(dia, peso)
## [1] 0.8626562
# Cuando pensamos en el modelo exponencial se obtiene un valor mucho más esperanzador, 
# acorde con el buen ajuste que vemos en la gráfica:
cor(dia, log(peso))
## [1] 0.9991654
# El cuadrado es:
cor(dia, log(peso))^2
## [1] 0.9983315

Apartado d.

Obtén estimaciones del peso para embriones de ocho días y cuarto, catorce días y medio y 20 días. Comenta la fiabilidad de estas estimaciones.

Solución:

# El valor para día = 20 sería extrapolación, asi que no lo calcularemos. 

# Para los otros dos usamos predict y luego tomamos la exponencial. La razón para hacer esto es que 
# nuestra recta predice valores del *logaritmo de peso*, no directamente del peso.
exp(predict(lm01, newdata = data.frame(dia = c(8.25, 14.5))))
##          1          2 
## 0.08449279 1.41604105

Apartado e.

Compara este modelo exponencial con una recta de regresión de peso frente a dia.

Solución:

# Aunque la conclusión parece evidente a partir del diagrama de dispersión, 
# vamos a añadir la recta de regresión de peso frente a dia a la gráfica.
(lm02 = lm(peso ~ dia))
## 
## Call:
## lm(formula = peso ~ dia)
## 
## Coefficients:
## (Intercept)          dia  
##     -1.8845       0.2351
plot(dia, peso)
diaNodes = seq(min(dia), max(dia), length.out = 30)
points(diaNodes, exp(lm01$coefficients[1] + lm01$coefficients[2] * diaNodes ), 
       type="l", lwd=2, col="blue")
abline(lm02, col="red", lwd=2)

# La conclusión es la esperada; el modelo exponencial es superior.

Cuestión 5.

Una compañía de electricidad quiere predecir el consumo mensual (en kw/h) de energía eléctrica de un hogar en función de la superficie (en \(m^2\)) de la casa y ha obtenido estos datos:

\(\quad\)
superficie 120 125 135 150 160 170 185 205 225 275
consumo 1182 1172 1264 1493 1571 1711 1804 1840 1956 1954
\(\quad\)

Estos datos están disponibles en el siguiente fichero adjunto.

Solución:

# Lectura de datos y creación de variables:
practica02c05 = read.table("./datos/practica02c05.csv", header = TRUE, sep=",")
(superficie = practica02c05$superficie)
##  [1] 120 125 135 150 160 170 185 205 225 275
(consumo = practica02c05$consumo)
##  [1] 1182 1172 1264 1493 1571 1711 1804 1840 1956 1954

Apartado a.

Dibuja el diagrama de dispersión para esas variables: superficie en el eje x, consumo en el eje y.

Solución:

# Como en ejercicios anteriores:
plot(superficie, consumo)

Apartado b.

Ajuste los datos por una curva exponencial de la forma: \[\text{consumo} = a\cdot \text{ superficie}^b\] donde de nuevo \(a\) y \(b\) son constantes y la idea vuelve a ser tomar logaritmos.

Solución:
Tomando logaritmos: \[\log(\text{consumo}) = \log a + b\cdot\log(\text{superficie})\] Así que podemos ajustra una recta a los valores de \(\log(\text{consumo})\) frente a los de \(\log(\text{superficie})\).

# Ajustamos esa recta
(lm01 = lm(log(consumo)  ~ log(superficie)))
## 
## Call:
## lm(formula = log(consumo) ~ log(superficie))
## 
## Coefficients:
##     (Intercept)  log(superficie)  
##          3.7242           0.7079
# Y representamos gráficamente el resultado:
(a = exp(lm01$coefficients[1]))
## (Intercept) 
##    41.43684
(b = lm01$coefficients[2])
## log(superficie) 
##       0.7078515
supfNodes = seq(min(superficie), max(superficie), length.out = 30)
(consumNodes = a * supfNodes^b)
##  [1] 1227.830 1266.294 1304.281 1341.817 1378.925 1415.625 1451.936
##  [8] 1487.876 1523.461 1558.707 1593.627 1628.233 1662.539 1696.555
## [15] 1730.291 1763.759 1796.966 1829.922 1862.634 1895.112 1927.361
## [22] 1959.389 1991.202 2022.807 2054.209 2085.415 2116.429 2147.256
## [29] 2177.902 2208.371
plot(c(superficie, supfNodes), c(consumo, consumNodes), type="n")
points(superficie, consumo)
points(supfNodes, consumNodes, type="l", lwd=2, col="blue")

# **Indicación:** En GeoGebra para obtener el mismo resultado puedes usar:
superficie = {120, 125, 135, 150, 160, 170, 185, 205, 225, 275}
consumo = {1182, 1172, 1264, 1493, 1571, 1711, 1804, 1840, 1956, 1954}
datos = Secuencia[(Elemento[superficie, i], Elemento[consumo, i]), i, 1, Longitud[consumo]]
AjustePotencia[datos]

Apartado b.

Calcule los valores del coeficiente de correlación \(R\) y \(R^2\). Interprételos.

Solución:

# Estos son R y R^2.
cor(log(superficie), log(consumo)) # R
## [1] 0.9393047
cor(log(superficie), log(consumo))^2 # R^2
## [1] 0.8822933
# Aunque los valores no son malos, el ajuste en la gráfica no es demasiado espectacular. ¿Hay puntos influyentes?

# También puede ser interesante aprender a fijarse en la información que proporciona 
# summary aplicado al resultado de lm.
(sumLm01 = summary(lm01))
## 
## Call:
## lm(formula = log(consumo) ~ log(superficie))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12238 -0.05026  0.02307  0.04152  0.08528 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.72417    0.46967   7.929 4.65e-05 ***
## log(superficie)  0.70785    0.09141   7.744 5.52e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07314 on 8 degrees of freedom
## Multiple R-squared:  0.8823, Adjusted R-squared:  0.8676 
## F-statistic: 59.97 on 1 and 8 DF,  p-value: 5.516e-05
# Por supuesto, los valores que se obtienen así coinciden con los de antes. Por ejemplo, R^2
sumLm01$r.squared 
## [1] 0.8822933

Apartado c.

Esime el consumo mensual de energía eléctrica de un hogar de 250 \(m^2\)

Solución:

# Como en otros ejercicios lo podemos hacer directamente o vía predict (y aplicando exp, porque el valor predicho es el logaritmo).
supf0 = 250
(pred0 = a * supf0^b)
## (Intercept) 
##    2064.297
exp(predict(lm01, newdata = data.frame(superficie = 250)))
##        1 
## 2064.297
# Como comentario adicional, este valor que estimamos está demasiado cerca del borde superior de la muestra.
# (y más aún teniendo en cuenta que parece haber un punto influyente).
# Es un buen momento para recordar que las predicciones del modelo son tanto peores cuanto más nos alejamos
# del centro de la muestra. Gráficamente si añadimos esa predicción al gráfico anterior es:

plot(c(superficie, supfNodes), c(consumo, consumNodes), type="n", asp=1)
points(supf0, pred0, col="orange", pch="*", cex=3)
points(superficie, consumo)
points(supfNodes, consumNodes, type="l", lwd=2, col="blue")

Cuestión 6.

Usa esta serie de datos para responder a las siguientes cuestiones.

\(\quad\)
x 1 3 6 8 11
y 11.9 13.1 15.2 18.5 29.1
\(\quad\)

Los datos están también disponibles en el siguiente fichero adjunto.

Solución:

# Lectura de datos y creación de variables:
practica02c06 = read.table("./datos/practica02c06.csv", header = TRUE, sep=",")
(x = practica02c06$x)
## [1]  1  3  6  8 11
(y = practica02c06$y)
## [1] 11.9 13.1 15.2 18.5 29.1

Apartado a.

Dibuja el diagrama de dispersión para esas variables.

Solución:

# Gráfica con plot:
plot(x, y)

Apartado b.

Usando el método de mínimos cuadrados ajuste a esos datos una parábola de la forma \[y = a + b\, x + c\,x^2\]

Solución:

# Una forma de ajustar usando R un polinomio es esta:  
(lm01 = lm(y ~ 1 + I(x) + I(x^2)))
## 
## Call:
## lm(formula = y ~ 1 + I(x) + I(x^2))
## 
## Coefficients:
## (Intercept)         I(x)       I(x^2)  
##     13.2003      -0.9918       0.2189
# A continuación podemos dibujar el resultado como en otros ejemplos, creando una lista de nodos,
# encontrando los valores que predice el modelo y uniéndolos en una curva:
xNodes =seq(min(x), max(x), length.out = 30)
yNodes = predict(lm01, newdata = data.frame(x=xNodes))
plot(x, y)
points(xNodes, yNodes, type = "l", col="blue")

Apartado c.

Estime el valor de \(y\) para \(x = 5.6\).

Solución:

# Usamos predict:
predict(lm01, newdata = data.frame(x = 5.6))
##        1 
## 14.51024

Cuestión 7.

Los datos de la tabla siguiente son el resultado de un estudio del efecto de la temperatura de cristalización primaria de una solución (medida en grados centígrados) sobre el contenido en fósforo (medido en gramos por litro):

\(\quad\)
temperatura -6 -3 0 3 6 9 12 15 20 25
fósforo 2.0 2.8 3.9 4.2 5.8 6.2 7.5 8.2 9.3 10.9
\(\quad\)

Estos datos están disponibles en el siguiente fichero adjunto.

Solución:

# Lectura de datos y creación de variables:
practica02c07 = read.table("./datos/practica02c07.csv", header = TRUE, sep=",")
(temperatura = practica02c07$temperatura)
##  [1] -6 -3  0  3  6  9 12 15 20 25
(fosforo = practica02c07$fosforo)
##  [1]  2.0  2.8  3.9  4.2  5.8  6.2  7.5  8.2  9.3 10.9

Apartado a.

Dibuja el diagrama de dispersión para esas variables (temperatura en el eje x, fósforo en el eje y).

Solución:

# Gráfica básica con plot:
plot(temperatura, fosforo, asp=1)

Apartado b.

Ajuste estos datos:

  • por una recta del tipo \(y = a + bx\).
  • por una función exponencial del tipo \(y = a\cdot e^{b\,x}\)
  • por una función parabólica del tipo \(y = a + b\, x + c\,x^2\).

Solución:

# Ajuste por recta:
lm01 = lm(fosforo ~ temperatura)
summary(lm01)
## 
## Call:
## lm(formula = fosforo ~ temperatura)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41088 -0.12507 -0.03329  0.14807  0.32493 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.746694   0.101162   37.04 3.10e-10 ***
## temperatura 0.288062   0.008087   35.62 4.22e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2438 on 8 degrees of freedom
## Multiple R-squared:  0.9937, Adjusted R-squared:  0.993 
## F-statistic:  1269 on 1 and 8 DF,  p-value: 4.223e-10
# Ajuste por exponencial:
lm02 = lm(log(fosforo) ~ temperatura)
summary(lm02)
## 
## Call:
## lm(formula = log(fosforo) ~ temperatura)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25094 -0.07668  0.03848  0.09991  0.18414 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.258902   0.062261   20.22 3.74e-08 ***
## temperatura 0.052469   0.004977   10.54 5.71e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.15 on 8 degrees of freedom
## Multiple R-squared:  0.9329, Adjusted R-squared:  0.9245 
## F-statistic: 111.1 on 1 and 8 DF,  p-value: 5.713e-06
# Ajuste por parábola:
lm03 = lm(fosforo ~ 1 + I(temperatura) + I(temperatura^2))
summary(lm03)
## 
## Call:
## lm(formula = fosforo ~ 1 + I(temperatura) + I(temperatura^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44260 -0.15743  0.05481  0.13622  0.27633 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.7508292  0.1052153   35.65 3.55e-09 ***
## I(temperatura)    0.2990413  0.0188584   15.86 9.62e-07 ***
## I(temperatura^2) -0.0005947  0.0009147   -0.65    0.536    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2531 on 7 degrees of freedom
## Multiple R-squared:  0.9941, Adjusted R-squared:  0.9924 
## F-statistic: 588.9 on 2 and 7 DF,  p-value: 1.585e-08
# Y además comparamos gráficamente los tres:

plot(temperatura, fosforo, asp=1)
tempNodes = seq(min(temperatura), max(temperatura), length.out = 30)
pred01 = predict(lm01, newdata = data.frame(temperatura=tempNodes))
points(tempNodes, pred01, type = "l", col="blue", lwd=5)
pred02 = exp(predict(lm02, newdata = data.frame(temperatura=tempNodes)))
points(tempNodes, pred02, type = "l", col="red", lwd=5)
pred03 = predict(lm03, newdata = data.frame(temperatura=tempNodes))
points(tempNodes, pred03, type = "l", col="orange", lwd=2)

Apartado b.

¿Cuál de los modelos de ajuste anteriores es el más adecuado para representar la relación entre las dos variables? ¿Por qué?

Solución:
La parábola y la recta son claramente mejores que el exponencial. Aunque el \(R^2\) de la parábola sea ligeramente mejor, el Principio de Parsimonia (a calidades similares, modelos simples mejor que modelos complicados) nos llevaría a elegir la recta.

Ejemplos adicionales para esta práctica.

## Este ejemplo ilustra porque no se puede usar sólo el coeficiente de correlación para juzgar la bondad del ajuste lineal

set.seed(2017) # Esta línea sirve para que tus resultados coincidan con los nuestros (reproducibilidad).

n = 150

x = sort(runif(n, min = -1, max = 1))

y = x +  sin(20 * x)/10 + rnorm(n)/50

plot(x, y)
abline(lm(y ~x))

# Como puede verse en la gráfica, la relación entre x e y es claramente no lineal,  a pesar de que:
cor(x, y)
## [1] 0.9919605
# Por eso, insistimos, el coeficiente de correlación NUNCA debe usarse como el único criterio para concluir que el ajuste por una recta es adecuado, sin acompañarlo de análisis exploratorios, gráficas, etc.

## Las cosas aún peden ser peores. Fíjate en este otro ejemplo:

set.seed(2017)

n = 150

x = sort(runif(n, min = -1, max = 1))

y = floor(3 * x) + rnorm(n)/20

plot(x, y)
abline(lm(y ~x))

# En este caso la gráfica muestra no sólo que la relación entre x e y es claramente no lineal, sino que además debería hacernos sospechar que la variable $y$ es discreta. Y en ese 
# caso usar un modelo de regresión lineal es una muy mala idea. Pero el coeficiente de correlación no detecta nada de esto.
cor(x, y)
## [1] 0.9863509

Práctica 3. DISTRIBUCIONES DE VARIABLES ALEATORIAS DISCRETAS.

Objetivos

  • Aprender a manejar las principales variables aleatorias discretas, empleándolas en los contextos adecuados.
  • Estudiar sus funciones de distribución, y calcular probabilidades, porcentajes y percentiles.
  • Generar y analizar números aleatorios asociados a una distribución dada.
  • Calcular e interpretar los puntos críticos de las distribuciones más importantes.

Cuestión 1.

Sea X una variable aleatoria con distribución binomial de parámetros \(n = 10\) y \(p = 0. 4\). Calcule las probabilidades siguientes:

  1. \(P(X < 4)\)
  2. \(P(X\leq 4)\)
  3. \(P(X \geq 8)\)
  4. \(P (2 \leq X < 7)\)
  5. \(P (4 < X \leq 9)\)
  6. \(P (3 \leq X \leq 6)\)

Solución:
En este ejercicio se trata de practicar el uso de las funciones dbinom y pbinom de R, que actúan como tablas para cualquier distribución binomial. En este caso usamos una binomial \(B(n = 10, p = 0.4)\).

# Empezamos definiendo los valores de n y p
n = 10
p = 0.4
# Apartado a
# P(X < 4) = P(X <= 3) = suma de probabilidades para X = 0, 1, 2, 3
pbinom(3, size = n, prob = p)
## [1] 0.3822806
sum(dbinom(0:3, size = n, prob = p))
## [1] 0.3822806
# Apartado b
# P(X <= 4)

pbinom(4, size = n, prob = p)
## [1] 0.6331033
sum(dbinom(0:4, size = n, prob = p))
## [1] 0.6331033
# Apartado c
# P(X >= 8) = 1 - P(X <= 7)

1 - pbinom(7, size = n, prob = p)
## [1] 0.01229455
sum(dbinom(8:n, size = n, prob = p))
## [1] 0.01229455
# Apartado d 
# P(2 <= X < 7) = P(X <= 6) - P(X <= 1)

pbinom(6, size = n, prob = p) - pbinom(1, size = n, prob = p)
## [1] 0.8988807
sum(dbinom(2:6, size = n, prob = p))
## [1] 0.8988807
# Apartado e 
# P(4 < X <= 9) = P(X <= 9) - P(X <= 4)

pbinom(9, size = n, prob = p) - pbinom(4, size = n, prob = p)
## [1] 0.3667919
# Apartado f
# P(3 <= X <= 6) = P(X <= 6) - P(X <= 2)

pbinom(6, size = n, prob = p) - pbinom(2, size = n, prob = p)
## [1] 0.7779484

Representación gráfica:
Creemos que es una buena idea que pienses cómo es la representación gráfica de la distribución binomial \(B(n = 10, p = 0.4)\).

barplot(dbinom(0:n, size = n, prob = p), names.arg = 0:n)

Cuestión 2.

La probabilidad de que un paciente se recupere de una rara enfermedad de la sangre es 0.4. Si se sabe que 100 personas contrajeron esa enfermedad, calcule la probabilidad de que sobrevivan:

  1. 30 pacientes o menos.
  2. Exactamente 45 pacientes.
  3. Entre 35 y 40, ambos valores incluidos.
  4. ¿Cuántos pacientes de los 100 enfermos se espera que sobrevivan a la enfermedad?

Solución:
Vamos a usar la variable \[ X = \text{número de pacientes que se recuperan} \] Si suponemos que la recuperación de cada paciente es independiente de la de los demás, entonces \(X\) es una variable binomial \(B(n = 100, p = 0.4)\).

n = 100
p = 0.4
# Apartado a
# P(X <= 30)

pbinom(30, size = n, prob = p)
## [1] 0.02478282
# Apartado b
# P(X = 45)

dbinom(45, size = n, prob = p)
## [1] 0.04781118
# Si te sorprende este resultado recuerda que se trata de *exactamente* 45, no 44 ni 46. La probabilidad de cada valor concreto de una binomial es cada vez más pequeña a medida que n aumenta.

# Apartado c
# P(35 <= X <= 40)

pbinom(40, size = n, prob = p) - pbinom(34, size = n, prob = p)
## [1] 0.412958
# Apartado d
# Se trata de calcular la media o valor esperado que en una binomial es simplemente:

n * p
## [1] 40

Representación gráfica:
De nuevo, es una buena idea representar gráficamente esta distribución binomial y compararla con la del ejercicio anterior.

barplot(dbinom(0:n, size = n, prob = p), names.arg = 0:n)

Fíjate en que a diferencia de lo que ocurría antes, la forma de la distribución empieza a recordar a la curva normal. Es muy importante que experimentes con esta gráfica cambiando el valor de \(p\). Prueba por ejemplo con \(p= 0.01\) y observa lo que sucede.

Cuestión 3.

La probabilidad de que el próximo nacimiento en un hospital sea de un niño es 0.52. En los próximos mil nacimientos que se produzcan en dicho hospital, calcule la probabilidad de que haya:

  1. Más de 480 niños.
  2. Exactamente 510 niñas.
  3. Más niños que niñas.

Solución:
Vamos a usar la variable \[ X = \text{número de niños (vs niñas) que nacen en el hospital} \] Si suponemos que el género de cada nacimiento es independiente del resto (y es un buen ejercicio que pienses cómo podría ocurrir lo contrario), entonces \(X\) es una variable binomial \(B(n = 1000, p = 0.52)\).

n = 1000
p = 0.52

# Apartado a
# P(X > 480) = 1 - P(X <= 480)

1 - pbinom(480, size = n, prob = p)
## [1] 0.9937742
# Apartado b
# P(X = 510 niñas) = P(X = n - 510 niños)

dbinom(n - 510, size = n, prob = p)
## [1] 0.004165995
# Como hemos comentado antes, cada valor concreto es pequeño. Incluso el valor más grande es pequeño, sólo un 2%:

dbinom(n * p, size = n, prob = p)
## [1] 0.02524521
# Apartado c
# P(mayoría de niños) = P(X >  500) = 1 - P(X <= 500)

1 - pbinom(500, size = n, prob = p)
## [1] 0.8914189
# Este último puede resultar chocante, así que el siguiente gráfico ilustra la razón de un valor tan alto de la probabilidad. Los valores por debajo de la barra azul corresponden a más niñas que niños:

bplot = barplot(dbinom(0:n, size = n, prob=p), plot = FALSE)
barplot(dbinom(0:n, size = n, prob=p), border = "grey", col="gray")
axis(1, at=seq(min(bplot), max(bplot), length.out = 11), labels = seq(0, n, 100))
abline(v = bplot[n*p +1], col="red", lwd=3)
abline(v = bplot[500], col="blue", lwd=3)

# Y como puede verse, casi todo el área corresponde a caso con más niños que niñas.

Cuestión 4.

La variable aleatoria X tiene una distribución de Poisson de parámetro 9.

  1. Halle la media de X.

  2. Calcule: \[P (X \leq 3),\qquad P (X \geq 6),\qquad P (4 \leq X \leq 12)\] Solución:
    Ahora vamos a practicar el uso de las funciones dpois y ppois que juegan en la distribución de Poisson el mismo papel que dbinom y pbinom en la binomial R. Sólo debemos introducir el valor del parámetro \(\lambda\).

lambda = 9

# Apartado a
# P(X <= 3)

ppois(3, lambda)
## [1] 0.02122649
sum(dpois(0:3, lambda))
## [1] 0.02122649
# P(X >= 6) = 1 - P(X <= 5

1 - ppois(5, lambda)
## [1] 0.8843095
# P(4 <= X <= 12) = P(X <= 12) - P(X <= 3)

ppois(12, lambda) - ppois(3, lambda)
## [1] 0.8545469

Cuestión 5.

Se sospecha que muchas muestras de agua, todas del mismo tamaño y tomadas de un río, han sido contaminadas por operarios irresponsables de una planta de tratamiento de aguas. Se contó el número de organismos coliformes en cada muestra. El número medio de organismos por muestra fue de 15. Suponiendo que el número de organismos se distribuye según una distribución de Poisson, calcule la probabilidad de que:

  1. Las dos próximas muestras contengan al menos 17 organismos.

Solución:
Si el número de bacterias en una muestra es una distribución de Poisson de parámetro 15, entonces el número de bacterias en dos muestras es una Poisson de parámetro \(\lambda = 30\).

# Apartado a
L = 15
n = 2

#P(X >= 17) = 1 - P(X <= 16)

1 - ppois(16, lambda = n * L)
## [1] 0.9961273

Atención:
Dentro de este mismo apartado también se puede pensar en que se pregunta por la probabilidad de que ambas muestras superen el número de 17. En ese caso tenemos que usar una Poisson de parámetro 15 y usando la independencia calcularíamos:

(1 - ppois(16, lambda =  L))^2
## [1] 0.1128132
  1. Las diez próximas muestras contengan al menos 100 organismos.

Atención:
En este caso el número total esperado de bacterias en 10 muestras es de 150 bacterias. Así que usamos una Poisson de parámetro 150 y calculamos:

n = 10

# P(X >= 100) = 1 - P(X <= 99)

1 - ppois(99, lambda = n * L)
## [1] 0.9999941

Cuestión 6.

En cierto cultivo, el número medio de células de Rickettsia typhi (células que causan el tifus) es de cinco por 20 \(\mu m^2\). a) ¿Cuántas células de ese tipo se esperan encontrar en un cultivo de 16 \(\mu m^2\)? b) ¿Qué probabilidad hay de que no se encuentre ninguna en un cultivo de 16 \(\mu m^2\)? c) ¿Y de que haya al menos nueve células en un cultivo de este tamaño?

Solución:
Si el número de células es 5 en 20 \(\mu m^2\) entonces cuando tenemos una superficie de 16 \(\mu m^2\) usaremos una distribución de Poisson con un parámetro \(\lambda\) calculado así:

L0 = 5/20
(L16 = 16 * L0)
## [1] 4

Solución:
En el apartado (a) se pide la media de la distribución de Poisson, que es precisamente el valor de \(\lambda\) que hemos calculado.

# Apartado b
# P(X = 0)

dpois(0, lambda = L16)
## [1] 0.01831564
ppois(0, lambda = L16)
## [1] 0.01831564
# Apartado c
# P(X >= 9) = 1 - P(X <= 8)

1 - ppois(8, lambda = L16)
## [1] 0.02136343

Representación gráfica:
Vamos a ver también en este caso la representación gráfica de la distribución, limitándola a los primeros valores porque la probabilidad \(P(X = k)\) de la distribución de Poisson decae rápidamente cuando \(k\) se hace grande. Juega con los valores de \(\lambda\) para ver cómo cambia la gráfica.

barplot(dpois(0:20, lambda = L16), names.arg = 0:20)

Cuestión 7.

Aproximación de una distribución binomial por una distribución de Poisson.

Se pretende comprobar que en las condiciones estudiadas en clase, la distribución de Poisson es una buena aproximación de la distribución binomial. Para ello:

  1. Represente gráficamente la variable binomial de parámetros n = 50 y p = 0. 1.
barplot(dbinom(0:50, size = 50, prob = 0.1))

  1. Represente gráficamente la variable de Poisson de parámetro \(\lambda = 50 · 0.1 = 5\). Haga que los valores de X varíen entre 0 y 50, de 10 en 10, y que los valores de Y varían entre 0 y 0.2, de 0.04 en 0.04.
barplot(dpois(0:50, lambda = 5))

  1. Compare ambas gráficas y haga un comentario.
barplot(rbind(dbinom(0:50, size = 50, prob = 0.1), dpois(0:50, lambda = 5)), col = c("red", "blue"), beside = TRUE)

Cuestión 8.

Un laboratorio afirma que una droga causa efectos secundarios en una proporción de 3 de cada 100 pacientes. Si se eligen 5 pacientes al azar a los que se les aplica la droga, calcula la probabilidad de que:

  1. Ningún paciente tenga efectos secundarios.
# Usando binomial
dbinom(0, size = 5, prob = 0.03)
## [1] 0.858734
# Usando Poisson
dpois(0, lambda = 5*0.03)
## [1] 0.860708
  1. Al menos dos tengan efectos secundarios.
# Usando binomial
1 - pbinom(1, size = 5, prob = 0.03)
## [1] 0.008472053
# Usando Poisson
1 - ppois(1, lambda = 5*0.03)
## [1] 0.01018583
  1. ¿Cuál es el número medio de pacientes que espera el laboratorio que sufran efectos secundarios si elige 250 pacientes al azar?
250*0.03
## [1] 7.5

Cuestión 9.

Algunas especies de paramecios producen y segregan partículas “letales” que causarían la muerte de un individuo susceptible, si este se pusiera en contacto con ellas. Todos los paramecios que no son capaces de producir partículas letales son susceptibles a ellas. El número medio de partículas letales emitidas por un paramecio de tales especies es de uno cada cinco horas.

  1. ¿Cuál es la probabilidad de que uno de estos paramecios no emita tales partículas en un periodo de dos horas y media?
L = 1/5
dpois(0, lambda = 2.5*L)
## [1] 0.6065307
  1. ¿Cuál es la probabilidad de que emita al menos una partícula letal?
1 - ppois(0, lambda = 2.5*L)
## [1] 0.3934693

Práctica 4. DISTRIBUCIONES DE VARIABLES ALEATORIAS CONTINUAS

Objetivos.

  • Aprender a manejar las principales variables aleatorias continuas, empleándolas en los contextos adecuados.
  • Estudiar sus funciones de densidad y calcular probabilidades, porcentajes y percentiles.
  • Generar y analizar números aleatorios asociados a una distribución dada.
  • Calcular e interpretar los puntos críticos de las distribuciones más importantes.

Cuestión 1.

Represente gráficamente de manera conjunta las siguientes distribuciones normales:

\[N (0, 1), \qquad N (0, 2), \qquad N (0, 3)\]

De la observación de las gráficas, ¿qué se puede decir acerca de la influencia que tienen en ellas los cambios de valores de la desviación típica?

Solución:
Este apartado y el siguiente los haremos usando el calculador de probabilidades de GeoGebra, al que se puede acceder instalando el programa en vuestros ordenadores o a través de la red en:

https://www.geogebra.org/m/YhMm8vgX

Cuestión 2.

Represente gráficamente de manera conjunta las siguientes distribuciones normales:

\[N (0, 1), \qquad N (1, 1), \qquad N (2, 1)\]

De la observación de las gráficas, ¿qué se puede decir acerca de la influencia que tienen en ellas los cambios de valores de la media?

Cuestión 3.

Calcule las siguientes probabilidades para una distribución N(0,1).

  1. \(\displaystyle P(Z < -1. 5)\)
  2. \(\displaystyle P(Z \geq 2)\)
  3. \(\displaystyle P(-0. 5 < Z \leq 1. 8)\)

Solución:

Usamos la función pnorm así:

Apartado a:

pnorm(-1.5)
## [1] 0.0668072

Apartado b:

1 - pnorm(2)
## [1] 0.02275013

Apartado c:

pnorm(1.8) - pnorm(-0.5)
## [1] 0.6555321

Cuestión 4.

Calcule las siguientes probabilidades para una distribución \(X \sim N (15, 2.5)\):

  1. \(\displaystyle P (X < 14)\)
  2. \(\displaystyle P (13 \leq X \leq 17)\)
  3. \(\displaystyle P (X \geq 18. 5)\)

Solución:

Se hace como el anterior con pnorm, pero indicando la media y desviación típica:

mu = 15
sigma = 2.5

Apartado a:

pnorm(14, mean = mu, sd = sigma)
## [1] 0.3445783

Apartado b:

pnorm(13, mean = mu, sd = sigma) - pnorm(17, mean = mu, sd = sigma)
## [1] -0.5762892

Apartado c:

1 - pnorm(18.5, mean = mu, sd = sigma)
## [1] 0.08075666

Cuestión 5.

Se admite que el peso de una población se distribuye según una \(N (67, 8)\). Si se toma una persona al azar, cuál es la probabilidad de que pese:

  1. Menos de 65
  2. Más de 70
  3. ¿Cuál es el peso que es superado solamente por el 10% de la población?
  4. ¿Cuál es el peso mínimo que no supera el 25% de la población?

Solución:

Trabajamos con una variable \(X\sim N(\mu = 67, \sigma = 8)\) y las funciones de R:

Apartado a:

mu = 67
sigma = 8

pnorm(65, mean = mu, sd = sigma)
## [1] 0.4012937

Apartado b:

1 - pnorm(70, mean = mu, sd = sigma)
## [1] 0.3538302

Apartado c:

Para el siguiente apartado usaremos qnorm, que permite averiguar cuál es el valor de X que deja en su cola izquierda una cantidad dada de probabilidad. En este caso se trata de calcular el valor \(x_0\) tal que: \[P(X > x_0) = 0.1\] Y esto es equivalente a pedir que sea: \[P(X < x_0) = 0.9\]

qnorm(0.9, mean = mu, sd = sigma)
## [1] 77.25241
1 - pnorm(70, mean = mu, sd = sigma)
## [1] 0.3538302

Cuestión 6.

El diámetro máximo de los hematíes de una persona con malaria por Plasmodium Vivas presenta las siguientes características: Si la célula está infectada dicha variable se distribuye de forma normal con media 7.6 micras y desviación típica 0.9 micras, y si la célula no está infectada dicha variable se distribuye de forma normal con media 9.6 micras y desviación típica 1.0 micras. Calcule la proporción de:

  1. Células no infectadas con un diámetro máximo mayor que 9.4 micras.
  2. Células no infectadas con un diámetro máximo inferior a 7 micras.
  3. Células infectadas con un diámetro máximo inferior a 9.4 micras

Solución:

En este ejercicio usamos una variable \(X_1\sim N(\mu_1 = 7.6, \sigma_1 = 0.9)\) para las celulas infectadas y otra \(X_2\sim N(\mu_2 = 9.6, \sigma_2 = 1.0)\) para las no infectadas:

mu1 = 7.6
sigma1 = 0.9

mu2 = 9.6
sigma2 = 1.0

Apartado a:

1 - pnorm(9.4, mean = mu2, sd = sigma2)
## [1] 0.5792597

Apartado b:

pnorm(7, mean = mu2, sd = sigma2)
## [1] 0.004661188

Apartado c:

pnorm(9.4, mean = mu1, sd = sigma1)
## [1] 0.9772499

Cuestión 7.

Siete mil corredores participan en una carrera local en la que se pueden clasificar para correr la maratón de la ciudad de New York si completan los 42 km en un tiempo inferior a 3 horas y 10 minutos. Solamente 6350 corredores han terminado la carrera. Si los tiempos empleados en completar los 42 km se distribuyen normalmente con una media de 3 horas 40 minutos y una desviación típica de 28 minutos, se pide:

  1. ¿Cuál es la probabilidad de que un corredor elegido al azar haya empleado menos de 3 horas en completar la carrera?
  2. ¿Cuántos corredores se han clasificado para la maratón de Nueva York?
  3. Si se clasifican más de 800 corredores, la organización aplicará como criterio de selección para acudir a la maratón de Nueva York haber conseguido un tiempo que esté incluido dentro del 5% de los mejores tiempos. ¿Qué tiempo máximo se debe emplear en la carrera local para ser seleccionado para la maratón de Nueva York con este criterio?

Solución:

Usamos una variable \(X\sim N(\mu = 3 + 40 / 60, \sigma = 28 / 60)\) para el tiempo empleado en finalizar la carrera:

mu = 3 + 40 / 60
sigma = 28 / 60

Apartado a:

pnorm(3, mean = mu, sd = sigma)
## [1] 0.07656373

Apartado b:
Empezamos averiguando que porcentaje de corredores emplean menos de 3 horas y 10 minutos, redondeándolo a entero:

y multiplicamos por el número que ha terminado para averiguar cuántos corresponden a ese porcentaje:

6350 * pnorm(3 + 10 / 60, mean = mu, sd = sigma)
## [1] 901.6263
round(6350 * pnorm(3 + 10 / 60, mean = mu, sd = sigma), 0)
## [1] 902

Apartado c:
El tiempo que representa al 5% más rápido de corredores es:

qnorm(0.05, mean = mu, sd = sigma)
## [1] 2.899068

Es decir, 2h y aproximadamante 54 minutos.

Cuestión 8.

Se considera que el nivel de plomo en sangre, en los niños, no es admisible para su salud si supera los 10 \(\mu\)g/dl. Se supone que la variable X: “cantidad de plomo presente en la sangre de un niño de una cierta población” tiene una distribución \(N(7, 2)\).

  1. Si se elige un niño al azar de esa población, ¿cuál es la probabilidad de que su nivel de plomo en sangre no sea admisible?
  2. Si se seleccionan al azar veinte niños de esa población, ¿cuál es la probabilidad de que haya entre 1 y 4 (ambos valores incluidos) con un nivel de plomo en sangre inadmisible?

Solución:

Usamos una variable \(X\sim N(\mu = 7, \sigma = 2)\) para el tiempo empleado en finalizar la carrera:

mu = 7
sigma = 2

Apartado a:

1 - pnorm(10, mean = mu, sd = sigma)
## [1] 0.0668072

Apartado b:
Para este apartado usamos una binomial \(B(n, p)\) con \(n = 20\) y \(p\) igual a la probabilidad que hemos calculado en el apartado a.

n = 20
p = 1 - pnorm(10, mean = mu, sd = sigma)
pbinom(4, size = n, prob = p) - pbinom(0, size = n, prob = p)
## [1] 0.7403082

Cuestión 9.

Unos laboratorios farmacéuticos almacenan en una nave cierto tipo de comprimidos que producen en dos plantas de fabricación. El 65% procede de la planta A y el 35% de la planta B. El peso, en miligramos, de los comprimidos procedentes de A tiene una distribución \(N (746, 15)\) y el peso de los comprimidos de B es \(N (754, 22)\).

  1. ¿Cuál es la probabilidad de que un comprimido de la planta A, elegido al azar, pese entre 751 y 760 mg.
  2. Si un comprimido elegido al azar en la nave de almacenamiento tiene un peso superior a 745 mg, ¿cuál es la probabilidad de que haya sido fabricado en la planta A?

Solución:

En este ejercicio usamos una variable \(X_1\sim N(\mu_1 = 746, \sigma_1 = 15)\) para los comprimidos de A y otra \(X_2\sim N(\mu_2 =754, \sigma_2 = 22)\) para los de B:

mu1 = 746
sigma1 = 15

mu2 = 754
sigma2 = 22

Apartado a:

pnorm(760, mean = mu1, sd = sigma1) - pnorm(751, mean = mu1, sd = sigma1)
## [1] 0.1941174

Apartado b:
Vamos a usar el teorema de Bayes así: \[ P(A | peso > 745) = \dfrac{P(peso > 745\, |\, A)\, P(A)}{P(peso > 745\, |\, A)\, P(A) + P(peso > 745\, |\, B)\, P(B)} \] Ahora las cuentas en R son fáciles:

(p745A = 1 - pnorm(745, mean = mu1, sd = sigma1))
## [1] 0.5265765
(p745B = 1 - pnorm(745, mean = mu2, sd = sigma2))
## [1] 0.6587635
pA = 0.65
(pB = 1 - pA)
## [1] 0.35
(pA745 = (p745A * pA) / (p745A * pA + p745B * pB)) 
## [1] 0.5975029

Cuestión 10.

Aproximación de una distribución binomial y una distribución de Poisson por una distribución normal.

  1. Represente gráficamente la variable binomial de parámetros \(n = 100\) y \(p = 0.3\).
  2. Represente gráficamente la variable \(N(100\cdot 0.3, \sqrt{100\cdot 0.3\cdot 0.7}) = N (30, 4.58)\).

Cuestión 11.

Aproximación de una distribución binomial y una distribución de Poisson por una distribución normal.

  1. Represente gráficamente la variable de Poisson de parámetros \(\lambda = 64\).
  2. Represente gráficamente la variable \(N (64, \sqrt{64}) = N (64, 8)\). Compare ambas gráficas.

Apartados a y b:
Ya hicimos ejercicios parecidos en la práctica anterior, así que aquí nos limitaremos a la representación gráfica. No te preocupes si no entiendes todo el código, son básicamente ajsutes de escala por razones técnicas.

bpPoisson = barplot(dpois(0:128, lambda =  64), ylim=c(0, 0.06))
(escala = (2 * 64)/diff(range(bpPoisson)))
## [1] 0.8333333
centro = bpPoisson[64]
densNorm = function(x)dnorm(x, mean = 64, sd = 8)
x = seq(from = 64 - 24, to = 64 + 24, length.out = 100)
y = densNorm(x)
points(x/escala, y, type = "l", col="red", lwd=7)

Cuestión 12.

  1. Calcule las siguientes probabilidades:
    \(P (T_{24} < 0.5)\)
    \(P (-0.3 < T_{13} \leq 1)\)
  2. Calcule las siguientes probabilidades:
    \(P (\chi^2_{35} \geq 12)\)
    \(P (10 < \chi^2_{11} < 30)\)

Solución:

Usando las funciones de R correspondientes:

# Apartado a
pt(0.5, df = 24)
## [1] 0.6891856
pt(1, df = 13) - pt(1, df = -0.3)
## Warning in pt(1, df = -0.3): NaNs produced
## [1] NaN
# Apartado b
1 - pchisq(12, df = 35)
## [1] 0.9998996
pchisq(30, df = 11) - pchisq(10, df = 11)
## [1] 0.5288026

Cuestión 13.

Calcule las siguientes probabilidades:
\(P (F_{12,17} < 1.6)\)
\(P (0.3 \leq F_{20,18} < 1.8)\)

Solución:

Otra pregunta que se responde con las funciones básicas de R:

# Apartado a
pf(1.6, df1 = 12, df2 = 17)
## [1] 0.8172578
pf(1.8, df1 = 20, df2 = 18) - pf(0.3, df1 = 20, df2 = 18)
## [1] 0.8869758

Cuestión 14.

Generación de una muestra aleatoria de una población normal. Análisis descriptivo.

  1. Genere una muestra aleatoria de tamaño 100 de una distribución \(N(170, 4)\) con el nombre de muestrap4.

  2. Obtenga el histograma y el gráfico de densidad suavizada y compruebe que estos sugieren normalidad en los datos.

  3. Confirme la normalidad con los valores de la asimetría y curtosis tipificadas.

Práctica 5. INTERVALOS DE CONFIANZA.

Objetivos.

  • Determinar las condiciones de aplicabilidad de los intervalos de confianza y el tipo de intervalo a utilizar en cada situación.
  • Construir con el ordenador los intervalos de confianza para la media, la varianza y la proporción de variables aleatorias con distribución normal.
  • Aprender a utilizar el ordenador para construir los intervalos de confianza para la diferencia de medias (datos pareados y no pareados), cociente de varianzas, y diferencia de proporciones de dos variables aleatorias normales e independientes.
  • Interpretar los resultados obtenidos.

Cuestión 1.

El nivel de creatina sérica en la sangre se considera un buen indicador de la presencia o ausencia de enfermedades de riñón. Las personas sanas tienen generalmente bajas concentraciones de creatina sérica (aproximadamente un valor de \(\mu\) = 1), mientras que las personas enfermas la tienen alta. Se desea conocer la relación entre el abuso de analgésicos y las enfermedades de riñón. Para ello, se estudia a 11 personas que trabajan en una fábrica y que se sabe que abusan de los analgésicos. Sus niveles de creatina han sido:

\[0.9,\, 0.7,\, 1.0,\, 1.1,\, 1.4,\, 2.0,\, 1.5,\, 2.2,\, 0.8,\, 0.8,\, 1.4\]

Puedes cargar estos datos desde el siguiente fichero adjunto:

Calcule un intervalo de confianza del 95% para el nivel medio de creatina en personas que abusan de los analgésicos. ¿Es razonable pensar con ese nivel de confianza que el consumo de analgésicos aumenta el nivel de creatina en la sangre? Justifíquelo.

Solución:

Como la muestra es pequeña usamos la t de Student. El intervalo de confianza para la media será: \[\mu = \bar X \pm t_{\alpha/2; n- 1}\dfrac{s}{\sqrt{n}}\] Para calcular el intervalo basta con ejecutar el siguiente codigo:

datos = read.table(file = "datos/practica05c01.csv", header = TRUE)
creatina = datos[,1]

# Valores muestrales
(xbar = mean(creatina))
## [1] 1.254545
(s = sd(creatina))
## [1] 0.4987256
(n = length(creatina))
## [1] 11
# Valor crítico tc de la t de Student.
nc = 0.95
(alfa = 1 - nc)
## [1] 0.05
(tc = qt(1 - alfa / 2, df = n - 1))
## [1] 2.228139
# Intervalo
(intervalo = xbar + c(-1, 1) * tc * s / sqrt(n))
## [1] 0.919497 1.589594
# Alternativamente se puede usar:
t.test(creatina, alternative = "two.sided", conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  creatina
## t = 8.343, df = 10, p-value = 8.137e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.919497 1.589594
## sample estimates:
## mean of x 
##  1.254545

Cuestión 2.

Los datos siguientes han sido obtenidos de una muestra aleatoria simple de tamaño 30, de la distribución X: “porcentaje de aumento del contenido de alcohol en la sangre de una persona, después de ingerir cuatro cervezas”. \[ \bar X = 41.2,\qquad s = 2.1 \]

  1. Calcule un intervalo de confianza del 90% para el porcentaje medio de aumento de alcohol en la sangre de una persona, después de tomar cuatro cervezas.

Solución:

Con una muestra de ese tamaño lo más adecuado es usar la t de Student. Como tenemos los descriptores de los datos, pasamos directamente al cálculo del intervalo:

# Valores muestrales
xbar = 41.2
s = 2.1
n = 30

# Valor crítico tc de la t de Student.
nc = 0.90
(alfa = 1 - nc)
## [1] 0.1
(tc = qt(1 - alfa / 2, df = n - 1))
## [1] 1.699127
# Intervalo
(intervalo = xbar + c(-1, 1) * tc * s / sqrt(n))
## [1] 40.54854 41.85146
  1. ¿Hay que hacer alguna suposición previa sobre la distribución de la variable X? En caso afirmativo, indíquela.

Solución:
Puesto que la muestra es pequeña, es necesario suponer que se trata de datos procedentes de una población normal.

Cuestión 3.

En un estudio de prevalencia de factores de riesgo en una cohorte de 412 mujeres mayores de quince años de una determinada población, resultó que el 17.6% eran hipertensas. Determine un intervalo de confianza del 90% para la proporción de mujeres hipertensas en dicha población.
¿Se podría afirmar a ese nivel de confianza que el porcentaje de mujeres hipertensas de esa población es del 21%?

Solución:

Para empezar debemos calcular un intervalo de confianza para la proporción. Como la muestra es grande podemos usar la expresión: \[p = {\hat p} \pm z_{\alpha/2}\sqrt{\dfrac{{\hat p}{\hat q}}{n}}\] Así que el código es:

# Valores muestrales
pHat = 0.176
qHat = 1 - pHat
n = 412

# Valor crítico zc de la normal estándar Z.
nc = 0.90
(alfa = 1 - nc)
## [1] 0.1
(zc = qnorm(1 - alfa / 2))
## [1] 1.644854
# Intervalo
(intervalo = pHat  + c(-1, 1) * zc * sqrt(pHat * qHat / n))
## [1] 0.1451398 0.2068602

Atención:

La última pregunta del enunciado plantea en realidad un contraste de hipótesis sobre la proporción. Desaconsejamos realizar dichos contrastes mediante el intervalo de confianza para la proporción muestral. Es más adecuado plantearlos por los métodos que estudiamos en el tema de contrastes.

Cuestión 4.

La variable aleatoria que representa la capacidad de las probetas producidas en una determinada empresa del vidrio tiene una distribución aproximadamente normal. Una muestra aleatoria de 7 de ellas dio como resultado una cuasivarianza de 62 mililitros\(^2\). Dé una estimación, mediante un intervalo de confianza del 95%, de la varianza de la capacidad de las probetas que fabrica dicha empresa.

Solución:

Sabemos que la expresión del intervalo de confianza para la varianza en una población (aproximadamente) normal es: \[\dfrac{(n-1)s^2}{\chi^2_{k,\alpha/2}}\leq\sigma^2\leq\dfrac{(n-1)s^2}{\chi^2_{k,1-\alpha/2}} , ext{ con }k = n - 1.\] Así que el código es:

# Valores muestrales
varianza = 62
(s = sqrt(varianza))
## [1] 7.874008
n = 7

# Valores críticos de chi cuadrado.
nc = 0.95
(alfa = 1 - nc)
## [1] 0.05
(chi1 = qchisq(1 - alfa / 2, df = n - 1))
## [1] 14.44938
(chi2 = qchisq(alfa/2, df = n - 1))
## [1] 1.237344
# Intervalo
(intervalo = (n - 1) * varianza / c(chi1, chi2))
## [1]  25.74506 300.64390

Cuestión 5.

Se tiene la sospecha de que en las aguas de un embalse las concentraciones de nitritos superan el umbral establecido para la existencia de vida piscícola de tipo ciprinícolas, que es de 0.03 mg NO2/l o menos. Para tratar de verificar esta sospecha, se midieron los niveles de nitritos en diez puntos aleatorios del embalse y se obtuvieron los siguientes valores:

\[ 0.02,\, 0.05,\, 0.03,\, 0.05,\, 0.04,\, 0.06,\, 0.07,\, 0.03,\, 0.04,\, 0.03 \]

Puedes cargar estos datos desde el siguiente fichero adjunto:

¿Se puede afirmar al nivel de confianza del 90% que las concentraciones de nitritos superan el umbral establecido para que sea factible la existencia de vida piscícola en el embalse? Justifíquelo.

Solución:

Es muy parecido al primer ejercicio de esta práctica, pero está formulado como un contraste de hipotesis. En lugar de eso vamos a calcular un intervalo de confianza al 90%.

datos = read.table(file = "datos/practica05c05.csv", header = TRUE)
nitritos = datos[,1]

# Valores muestrales
(xbar = mean(nitritos))
## [1] 0.042
(s = sd(nitritos))
## [1] 0.01549193
(n = length(nitritos))
## [1] 10
# Valor crítico tc de la t de Student.
nc = 0.9
(alfa = 1 - nc)
## [1] 0.1
(tc = qt(1 - alfa / 2, df = n - 1))
## [1] 1.833113
# Intervalo
(intervalo = xbar + c(-1, 1) * tc * s / sqrt(n))
## [1] 0.03301962 0.05098038
# Alternativamente se puede usar:
t.test(nitritos, alternative = "two.sided", conf.level = 0.9)
## 
##  One Sample t-test
## 
## data:  nitritos
## t = 8.5732, df = 9, p-value = 1.268e-05
## alternative hypothesis: true mean is not equal to 0
## 90 percent confidence interval:
##  0.03301962 0.05098038
## sample estimates:
## mean of x 
##     0.042

Cuestión 6.

Se ha aplicado un fármaco a una muestra de 200 enfermos y se ha observado una respuesta positiva en 140 de ellos. Estímese mediante un intervalo de confianza del 99%, la proporción de enfermos que responderían positivamente si este medicamento se aplicase a la población de la que se ha extraído la muestra.

Solución:

Se trata de un intervalo de confianza para la proporción, que ya hemos visto en un ejercicio anterior.

# Valores muestrales
pHat = 140 / 200
qHat = 1 - pHat
n = 200

# Valor crítico zc de la normal estándar Z.
nc = 0.99
(alfa = 1 - nc)
## [1] 0.01
(zc = qnorm(1 - alfa / 2))
## [1] 2.575829
# Intervalo
(intervalo = pHat  + c(-1, 1) * zc * sqrt(pHat * qHat / n))
## [1] 0.6165336 0.7834664
# Alternativamente
library(DescTools) # Debes instalar esta librería primero
## Warning: package 'DescTools' was built under R version 3.4.3
BinomCI(140, 200, conf.level = 0.99, method = "wald")
##      est    lwr.ci    upr.ci
## [1,] 0.7 0.6165336 0.7834664

Cuestión 7.

Se desea comparar los fármacos amlodipino y enalapril para el tratamiento de la hipertensión arterial. De dos poblaciones independientes se extraen dos muestras aleatorias simples de nueve y diez pacientes, respectivamente. A los primeros se les trata con amlodipino y a los segundos con enalapril. Posteriormente, se les mide el grado de disminución de la tensión arterial (en mg de mercurio), con los resultados siguientes:

\(\quad\)
Amlodipino 6 5 8 3 6 6 5 7 2
Enalapril 8 10 9 7 6 8 9 8 6 4
\(\quad\)
  1. Calcule un intervalo de confianza para la diferencia de las disminuciones medias de la presión arterial, del 90% de confianza. Se suponen distribuciones normales para las variables.
  2. ¿Se puede afirmar que el enalapril es más efectivo que el amlodipino en la disminución de la tensión arterial? Justifique la respuesta.

Solución:

Al tratarse de muestras pequeñas, la fórmula que se usa para el intervalo de confianza depende de si asumimos que las varianzas de ambas poblaciones son iguales o distintas. Vamos a usar la función t.test con la opción two.sided para obtener el intervalo. Veremos las dos posibilidades, asumiendo varianzas iguales y asumiendo varianzas distintas.

# Muestra 1.
muestra1 = c(6, 5 , 8, 3, 6, 6, 5, 7, 2)
muestra2 = c(8, 10, 9, 7, 6, 8, 9, 8, 6, 4)

t.test(muestra1, muestra2, alternative = "two.sided", conf.level = 0.90, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  muestra1 and muestra2
## t = -2.5866, df = 17, p-value = 0.0192
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
##  -3.6238238 -0.7095095
## sample estimates:
## mean of x mean of y 
##  5.333333  7.500000
t.test(muestra1, muestra2, alternative = "two.sided", conf.level = 0.90, var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  muestra1 and muestra2
## t = -2.5794, df = 16.568, p-value = 0.01977
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
##  -3.6300852 -0.7032481
## sample estimates:
## mean of x mean of y 
##  5.333333  7.500000

Como puedes ver, los dos intervalos son muy parecidos, pero el que asume varianzas distintas es más ancho. En general, si no dispones de ningún indicio sobre la igualdad de las varianzas es más prudente asumir que no son iguales; de esa forma el intervalo es más ancho y por tanto más fiable.
Alternativamente puedes usar los ficheros plantilla del Tutorial 9 para hacer un contraste de hipótesis sobre la igualdad de varianzas antes de calcular estos intervalos. También puedes hacerlo con var.test así:

var.test(muestra1, muestra2, alternative = "two.sided", conf.level = 0.95)
## 
##  F test to compare two variances
## 
## data:  muestra1 and muestra2
## F = 1.1053, num df = 8, denom df = 9, p-value = 0.8767
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2694479 4.8158892
## sample estimates:
## ratio of variances 
##           1.105263
sd(muestra1)
## [1] 1.870829
sd(muestra2)
## [1] 1.779513

El p-valor de este contraste nos lleva a rechazar la hipótesis alternativa, así que en este caso podríamos asumir que las varianzas son iguales (hemos incluido el cálculo de las cuasivarianzas muestrales para confirmar esto) y tomar el intervalo de confianza más pequeño.

Cuestión 8.

Se está realizando un estudio para considerar el efecto que tiene sobre el feto el hecho de que la madre fume. Para el estudio se tomó una muestra de 121 madres fumadoras y otra de 121 madres no fumadoras. La variable de interés es el peso en gramos del niño al nacer y se supone que esta variable se distribuye normalmente. Los datos obtenidos de las muestras son:

\(\quad\)
No fumadoras \(\bar X_1 = 348.0\)g \(s_1 = 8.68\)g \(n_1 = 121\)
Fumadoras \(\bar X_2 = 326.5\)g \(s_2 = 11.02\)g \(n_2 = 121\)
\(\quad\)
  1. Construya un intervalo de confianza del 90% para la diferencia de los pesos medios de los niños al nacer.
  2. ¿Qué puede concluirse acerca de la influencia que pueda tener en el peso del niño el que la madre sea fumadora?

Solución:

Al tratarse de muestras pequeñas, empezamos contrastando la igualdad de las varianzas. El fichero enlazado a continuación contiene la plantilla necesaria con los datos ya cumplimentados.

El p-valor que se obtiene es 0.0094. Por lo tanto (rechazamos la hipótesis nula de igualdad de varianzas) asumiremos que las varianzas son distintas.
Una vez hecho esto podemos usar el fichero plantilla correspondiente, enlazado aquí debajo con los datos, para calcular el intervalo de confianza deseado. Se obtiene el intervalo \((19.39, 23.61)\).

Corrección:

El ejercicio se ha resuelto como si las muestras fueran pequeñas, cuando en realidad no es necesario suponer eso. Con esos tamaños muestrales (>100) podemos usar la distribución \(Z\) y, mas importante, ahorrarnos el contraste de igualdad de varianzas. El siguiente fichero plantilla adjunto calcula el intervalo teniendo eso en cuenta.

Se obtiene el intervalo \((19.40, 23.60)\). Como puede verse, la diferencia con el uso de la \(t\) es muy pequeña: como sabemos, a medida que el tamaño muestral aumenta la \(t\) se parece cada vez más a la \(Z\) (recuerda que los intervalos que usan \(t\) son siempre más anchos que los que usan \(Z\)). En cualquier caso, la ventaja evidente es que nos ahorramos el contraste de igualdad de varianzas.

Cuestión 9.

Un laboratorio farmacéutico produce un mismo fármaco en dos plantas A y B. Pretende estudiar si la proporción de unidades fabricadas que mantienen su eficacia después de su fecha de caducidad es similar en ambas plantas. Para ello, se tomó una muestra aleatoria de 36 unidades de la planta A y se observó que 12 de ellas mantuvieron sus propiedades después de su fecha de caducidad. Otro tanto ocurrió con 16 unidades de otra muestra aleatoria de 36 unidades de la planta B.

  1. Calcule un intervalo de confianza del 99% para la diferencia de proporciones de fármacos caducados que mantienen sus propiedades.
  2. Justifique si la diferencia es o no significativa

Solución:
Se trata de calcular un intervalo de confianza para la diferencia de proporciones. Puede usarse la plantilla del Tutorial09 o bien hacerlo directamente así:

n1 = 36
n2 = 36
pHat1 = 12 / 36
qHat1 = 1 - pHat1
pHat2 = 16 / 36
qHat2 = 1 - pHat2
nc = 0.99
alfa = 1 - nc
za2 = qnorm(1 - alfa/2)
(intervalo = (pHat1 - pHat2) + c(-1, 1) * za2 * sqrt((pHat1 * qHat1 / n1) + (pHat2 * qHat2 / n2)))
## [1] -0.4051570  0.1829348

Alternativamente se puede usar la función prop.test así (atención a la opcion correct = FALSE):

prop.test(x = c(12, 16), n = c(36, 36), alternative = "two.sided", 
          conf.level = 0.99, correct = FALSE)
## 
##  2-sample test for equality of proportions without continuity
##  correction
## 
## data:  c(12, 16) out of c(36, 36)
## X-squared = 0.93506, df = 1, p-value = 0.3336
## alternative hypothesis: two.sided
## 99 percent confidence interval:
##  -0.4051570  0.1829348
## sample estimates:
##    prop 1    prop 2 
## 0.3333333 0.4444444

Cuestión 10.

Ocho personas seropositivas tomaron parte en un ensayo clínico con el objeto de determinar la efectividad de una droga para retardar la aparición del SIDA. Un indicador de esta aparición son las células CD4 del sistema inmunitario. Para cada paciente se contó el número de células CD4 antes y después del tratamiento con los siguientes resultados:

\(\quad\)
Antes 415 356 339 188 256 296 249 303
Después 474 329 555 282 423 323 256 431
\(\quad\)

Puedes cargar estos datos desde el siguiente fichero adjunto:

¿Se puede afirmar al nivel de confianza del 99% que el tratamiento retarda la aparición del SIDA? Justifíquelo.

Solución:
Comentario: Aunque aparece en la práctica 5 este ejercicio está planteado como un contraste de hipótesis. Las prácticas de la asignatura tradicionalmente mezclan los dos temas (los apartados b de muchos ejercicios de la práctica 5 son contrastes de hipótesis y en esta práctica hemos preferido ignorarlos). En el caso concreto de este ejercicio ni siquiera se pide un intervalo. En la solución que aparece aquí abajo hemos optado por calcular un intervalo de confianza al 99% para la diferencia de proporciones (emparejadas).

Cuidado: Se trata de calcular un intervalo de confianza para la diferencia de medias en muestras emparejadas. Lo hacemos, por ejemplo, con t.test de la manera más rápida. Fíjate en que para eso usamos la opción paired = TRUE y en que hemos ordenado las muestras de manera que el intervalo de confianza se refiere a la diferencia \[\text{media después} - \text{media antes}\] Es importante tenerlo en cuenta al interpretar el intervalo.

datos = read.table(file = "datos/practica05c10.csv", header = TRUE, sep=",")
antes = datos$antes
despues = datos$despues

t.test(despues, antes, alternative = "two.sided", paired = TRUE, conf.level = 0.99)
## 
##  Paired t-test
## 
## data:  despues and antes
## t = 2.8509, df = 7, p-value = 0.02466
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
##  -19.08281 186.83281
## sample estimates:
## mean of the differences 
##                  83.875

Otra opción igualmente válida es empezar calculando la diferencia, reduciendo así el problema a un problema de una muestra, y entonces calculando un intervalo de confianza para esas diferencias.

t.test(despues  - antes, alternative = "two.sided", conf.level = 0.99)
## 
##  One Sample t-test
## 
## data:  despues - antes
## t = 2.8509, df = 7, p-value = 0.02466
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
##  -19.08281 186.83281
## sample estimates:
## mean of x 
##    83.875

El resultado es evidentemente el mismo, pero la forma de usar t.test cambia ligeramente.

Cuestión 11.

Para estudiar la eficacia de un fármaco antidoloroso se eligieron aleatoriamente 29 mujeres que sufrían calambres uterinos después del parto. De ellas, 16 fueron tratadas con el medicamento y las 13 restantes con un placebo. La mejoría experimentada se midió en una escala entre 0 (ninguna mejoría) y 56 (mejoría completa durante 8 horas). Los resultados obtenidos fueron:

\(\quad\)
Fármaco 1 \(\bar X_2 = 31.96\) \(s_1 = 12.05\)g \(n_1 = 16\)
Fármaco 2 \(\bar X_2 = 25.32\) \(s_2 = 13.38\)g \(n_2 = 13\)
\(\quad\)
  1. Construya un intervalo de confianza del 98% para la diferencia de las mejorías medias entre las mujeres tratadas con el fármaco y las tratadas con placebo. (Se supone normalidad en los datos)
  2. ¿Se puede considerar que el fármaco es eficaz en la mejora del dolor producido por los calambres postparto? Justifíquelo

Solución:
Este ejercicio se puede hacer usando plantillas de forma análoga a como hemos hecho el Ejercicio 8. Pero vamos a aprovecharlo para presentar otra manera de proceder. Cuando no tenemos las muestras originales podemos simularlas. Para eso vamos a utilizar la función mvrnorm de la libraría MASS (viene instalada por defecto en R). Sería así:

n1 = 16
xbar1 = 31.96
s1 = 12.05

n2 = 13
xbar2 = 25.32
s2 = 13.38

library(MASS)
(muestra1 = mvrnorm(n = n1, mu = xbar1, Sigma = s1^2, empirical = TRUE)[,1])
##  [1] 21.60566 27.18233 35.80078 65.68955 34.56249 20.53763 31.55151
##  [8] 43.91405 43.94779 24.38885 23.64863 20.26590 28.40156 40.97587
## [15] 25.29708 23.59031
(muestra2 = mvrnorm(n = n2, mu = xbar2, Sigma = s2^2, empirical = TRUE)[,1])
##  [1]  8.17815 13.20956 17.18972 29.22121 37.48785 38.97519 10.57455
##  [8] 10.80379 19.66896 46.84818 42.84982 21.83558 32.31744

Fíjate en que Sigma se refiere a la cuasivarianza (hay que elevar al cuadrado) y que hemos usado empirical=TRUE. El detalle de [,1] se debe a que mvrnorm produce como resultado una matriz, y tenemos que tomr la primera (y única) columna como vector. Como ves, hemos obtenido dos muestras de datos. Ahora puedes comprobar que tienen los valores deseados. Por ejemplo:

length(muestra1)
## [1] 16
mean(muestra1)
## [1] 31.96
sd(muestra1)
## [1] 12.05

Análogamente:

length(muestra2)
## [1] 13
mean(muestra2)
## [1] 25.32
sd(muestra2)
## [1] 13.38

La ventaja de haber hecho esto es que ahora podemos usar var.test para decidir si las varianzas son iguales:

var.test(muestra1, muestra2, alternative = "two.sided", conf.level = 0.95)
## 
##  F test to compare two variances
## 
## data:  muestra1 and muestra2
## F = 0.81108, num df = 15, denom df = 12, p-value = 0.6915
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2552802 2.4034490
## sample estimates:
## ratio of variances 
##          0.8110766

Como ves, con ese p-valor asumiremos que son iguales. Así que ahora calculamos el intervalo para la diferencia de medias con t.test así:

t.test(muestra1, muestra2, alternative = "two.sided", conf.level = 0.98, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  muestra1 and muestra2
## t = 1.4048, df = 27, p-value = 0.1715
## alternative hypothesis: true difference in means is not equal to 0
## 98 percent confidence interval:
##  -5.047178 18.327178
## sample estimates:
## mean of x mean of y 
##     31.96     25.32

El intervalo es \((-5.047, 18.33)\)

Cuestión 12.

Se desea estudiar si un aumento de calcio en la dieta reduce la presión sanguínea. Para ello, se realizó un experimento con 20 personas. Diez de ellas (grupo 1) tomaron un suplemento de calcio durante doce semanas mientras que las diez restantes (grupo 2) tomaron un placebo. Tras las 12 semanas, se midió la disminución de la presión sistólica de la sangre de las 20 personas, con los siguientes resultados

\(\quad\)
Grupo 1 7 -4 18 17 -3 -5 1 10 11 -2
Grupo 2 -1 12 -1 -3 3 -5 5 2 -11 -3
\(\quad\)

Puedes cargar estos datos desde el siguiente fichero adjunto:

  1. Calcule un intervalo de confianza del 90% para la diferencia entre la disminución media la presión sanguínea de las personas que toman calcio y la de las personas que no toman calcio.
  2. ¿Es cierto que un aumento de calcio en la dieta reduce la presión sanguínea? Justifíquelo.

Solución corregida:

No se trata de muestras emparejadas, como se afirmaba erróneamente en la primera versión. Son dos muestras independientes y se debe calcular el intervalo de confianza para la diferencia de medias con dos muestras pequeñas independientes. Es decir, empezamos contrastando la igualdad de varianzas:

var.test(grupo1, grupo2, alternative = "two.sided", conf.level = 0.9)
## 
##  F test to compare two variances
## 
## data:  grupo1 and grupo2
## F = 1.9793, num df = 9, denom df = 9, p-value = 0.3237
## alternative hypothesis: true ratio of variances is not equal to 1
## 90 percent confidence interval:
##  0.6226339 6.2919403
## sample estimates:
## ratio of variances 
##           1.979287

Como el p-valor es grande, asumimos por tanto varianzas iguales. Ahora obtenemos el intervalo de diferencia de medias con:

t.test(grupo1, grupo2, alternative = "two.sided", conf.level = 0.9, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  grupo1 and grupo2
## t = 1.533, df = 18, p-value = 0.1427
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
##  -0.6822 11.0822
## sample estimates:
## mean of x mean of y 
##       5.0      -0.2

También puedes hacer estos cálculos usando ficheros plantilla. Se adjuntan las dos plantillas necesarias con los datos de este ejercicio ya incorporados.

Contraste de igualdad de varianzas:

Intervalo para la diferencia de medias:

Cuestión 13.

Para comprobar si un tratamiento con ácidos grasos es eficaz en pacientes con eczema atípico, se tomaron 10 pacientes con eczema y se les sometió durante tres semanas a un tratamiento ficticio (placebo) y durante las tres siguientes a un tratamiento con ácidos grasos. Tras cada periodo, un médico ajeno al proyecto evaluó la importancia del eczema en una escala de 0 (no eczema) a 10 (tamaño máximo de eczema).

\(\quad\)
Placebo 6 8 4 8 5 6 5 6 4 5
Tratamiento 5 6 4 5 3 6 6 2 2 6
\(\quad\)

Puedes cargar estos datos desde el siguiente fichero adjunto:

  1. Calcule un intervalo de confianza del 90% para la diferencia entre la importancia media del eczema atípico después del tratamiento ficticio (placebo) y después del tratamiento con ácidos grasos.
  2. Es cierto que el tratamiento del eczema atípico con ácidos grasos es eficaz? Justifíquelo

Solución corregida:
Este sí es el ejercicio de muestras emparejadas y el intervalo se obtiene así.

datos = read.table(file = "datos/practica05c13.csv", header = TRUE, sep=",")
placebo = datos$placebo
tratamiento = datos$tratamiento
t.test(placebo, tratamiento, alternative = "two.sided", paired = TRUE, conf.level = 0.95)
## 
##  Paired t-test
## 
## data:  placebo and tratamiento
## t = 2.25, df = 9, p-value = 0.051
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.00648382  2.40648382
## sample estimates:
## mean of the differences 
##                     1.2

Cuestión 14.

Cierta industria fabricante de lentes para gafas ha realizado un experimento para comparar dos tipos (1 y 2) de recubrimientos antirreflectantes y elegir aquel que sufra un menor desgaste. Ha seleccionado 12 personas al azar que requerían este tipo de lentes; ha preparado para las primeras 6 personas sendas gafas con un recubrimiento de tipo 1, y para las otras 6, sendas gafas con recubrimiento de tipo 2. Después de seis meses de uso se ha medido el desgaste sufrido en los cristales, con el resultado siguiente para cada una de las 12 personas.

\(\quad\)
Desgaste lente 4.1 4.4 4.8 4.5 4.3 5.2 5.0 5.5 4.0 4.6 5.9 7.0
Tipo de lente 1 1 2 1 2 1 2 2 1 2 1 2
\(\quad\)

Puedes cargar estos datos desde el siguiente fichero adjunto:

  1. Calcule un intervalo de confianza del 99% para la diferencia de desgastes medios en las lentes con los dos tipos de recubrimiento.
  2. ¿Se podría afirmar que con el de tipo 2 se obtiene un menor desgaste que con el de tipo 1? Justifíquelo.

Solución:
Vamos a leer los datos y a usar la variable tipo para separar las dos muestras. Si el siguiente código te resulta confuso recuerda que siempre te queda la opción de introducir los datos a mano.

datos = read.table(file = "datos/practica05c14.csv", header = TRUE, sep=",")
(tipo1 = datos$desgaste[datos$tipo == 1])
## [1] 4.1 4.4 4.5 5.2 4.0 5.9
(tipo2 = datos$desgaste[datos$tipo == 2])
## [1] 4.8 4.3 5.0 5.5 4.6 7.0

Solución:
A partir de aquí es muy parecido a los otros intervalos para diferencia de medias que hemos hecho usando primero var.test para decidir sobre la igualdad de varianzas y luego t.test: para calcular el intervalo.

var.test(tipo1, tipo2, alternative = "two.sided", conf.level = 0.95)
## 
##  F test to compare two variances
## 
## data:  tipo1 and tipo2
## F = 0.56773, num df = 5, denom df = 5, p-value = 0.5495
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.07944307 4.05721890
## sample estimates:
## ratio of variances 
##          0.5677305

Visto el p-valor asumiremos que las varianzas son iguales. Así que ahora calculamos el intervalo para la diferencia de medias con t.test así:

t.test(tipo1, tipo2, alternative = "two.sided", conf.level = 0.99, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  tipo1 and tipo2
## t = -1.0425, df = 10, p-value = 0.3217
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
##  -2.087332  1.053998
## sample estimates:
## mean of x mean of y 
##  4.683333  5.200000

Cuestión 15.

La siguiente tabla representa las presiones sanguíneas sistólicas de 10 individuos alcohólicos rehabilitados, antes y después de haber dejado la bebida.

\(\quad\)
Individuo 1 2 3 4 5 6 7 8 9 10
Antes 140 165 160 160 175 190 170 175 155 160
Después 145 150 150 160 170 175 160 165 145 170
\(\quad\)

Puedes cargar estos datos desde el siguiente fichero adjunto:

Estime mediante un intervalo de confianza del 95% el cambio de la presión sistólica que produce el abandono del alcohol

Solución:
Es muy parecido a cuestiones previas de muestras emparejadas.

datos = read.table(file = "datos/practica05c15.csv", header = TRUE, sep=",")
antes = datos$antes
despues = datos$despues

t.test(antes, despues, alternative = "two.sided", paired = TRUE, conf.level = 0.95)
## 
##  Paired t-test
## 
## data:  antes and despues
## t = 2.25, df = 9, p-value = 0.051
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0324191 12.0324191
## sample estimates:
## mean of the differences 
##                       6

Cuestión 16.

En un sondeo de opinión llevado a cabo en el año 2002 con 1000 individuos de una ciudad se encontró que el 55% de ellos eran favorables a la construcción de un nuevo aeropuerto. Dicho sondeo se repitió en el año 2003 tras una intensa campaña publicitaria, encontrándose esta vez un 60% de individuos favorables a esta opción, entre 1500 encuestados.

  1. Construya un intervalo de confianza al 95% para la diferencia de los porcentajes de individuos favorables al nuevo aeropuerto en los años 2003 y 2002.
  2. ¿Podría afirmarse al nivel de confianza del 99% que la campaña publicitaria en favor de la construcción del nuevo aeropuerto ha podido tener efecto?

Solución:
Es muy parecido a la Cuestión 9, un intervalo de confianza para la diferencia de proporciones. Como decíamos allí, puede usarse la plantilla del Tutorial09, que incluimos aquí con los datos ya cumplimentados:

Obtendrás el intervalo \(( -0.08957, -0.01043)\). Alternativamente puedes usar cálculo directo:

n1 = 1000
n2 = 1500
pHat1 = 0.55
qHat1 = 1 - pHat1
pHat2 = 0.60
qHat2 = 1 - pHat2
nc = 0.95
alfa = 1 - nc
za2 = qnorm(1 - alfa/2)
(intervalo = (pHat1 - pHat2) + c(-1, 1) * za2 * sqrt((pHat1 * qHat1 / n1) + (pHat2 * qHat2 / n2)))
## [1] -0.08956507 -0.01043493

y por supuesto se puede usar la función prop.test. Recuerda la opcion correct = FALSE y mira como hemos obtenido el número de éxitos de cada muestra a partr de las proporciones muestrales:

prop.test(x = c(0.55 * 1000, 0.60 * 1500), n = c(1000, 1500), alternative = "two.sided", 
          conf.level = 0.95, correct = FALSE)
## 
##  2-sample test for equality of proportions without continuity
##  correction
## 
## data:  c(0.55 * 1000, 0.6 * 1500) out of c(1000, 1500)
## X-squared = 6.1576, df = 1, p-value = 0.01308
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.08956507 -0.01043493
## sample estimates:
## prop 1 prop 2 
##   0.55   0.60

Práctica 6. CONTRASTES DE HIPÓTESIS PARAMÉTRICOS.

Objetivos.

  • Determinar las condiciones de aplicabilidad de los contrastes de hipótesis paramétricos y el tipo de contraste a utilizar en cada situación.
  • Realizar contrastes de hipótesis paramétricos para la media y la varianza de una población normal.
  • Realizar contrastes de hipótesis paramétricos para una proporción.
  • Realizar contrastes de hipótesis paramétricos para comparar dos medias de variables normales en muestras independientes y en muestras pareadas.
  • Realizar contrastes de hipótesis para comparar dos proporciones.
  • Interpretar los resultados obtenidos.

Cuestión 1.

A dieciocho esquiadores de fondo que han realizado un determinado recorrido se les mide la CPK en sangre (la cantidad de enzima CPK en sangre es una medida de estrés muscular). Los datos son los siguientes:

\[ \begin{array}{rrrrrrrrr} 520 & 300 & 480 & 1040 & 360 & 580 & 440 & 180 & 490\\ 520 & 380 & 640 & 1360 & 240 & 420 & 280 & 400 & 260 \end{array} \]

Puedes cargar estos datos desde el siguiente fichero adjunto:

  1. ¿Se puede afirmar al nivel de significación del 5% que la cantidad media de CPK en sangre de esquiadores de fondo que realicen ese mismo recorrido es 660?
  2. ¿Qué suposición debe hacerse acerca de la distribución de la variable aleatoria: “cantidad de CPK en la sangre de los esquiadores de fondo después del recorrido”?

Solución:

Apartado (a)
Haremos un contraste de hipótesis bilateral para la media con hipótesis nula \[H_0 = \{\mu = \mu_0\}\] donde \(\mu_0 = 660\). Empezamos leyendo los datos:

datos = read.table(file = "datos/practica06c01.csv", header = TRUE)
cpk = datos[ ,1]

Ahora vamos a calcular los descriptores muestrales:

(n = length(cpk))
## [1] 18
(xbar = mean(cpk))
## [1] 493.8889
(s = sd(cpk))
## [1] 289.8101

Con eso podemos calcular el valor del estadístico del contraste así:

mu0 = 660
estadistico = abs(xbar - mu0) / (s/sqrt(n))

Fíjate en que hemos tomado valor absoluto porque se trata de un contraste bilateral. Y puesto que la muestra es pequeña podemos usar la \(t\) de Student para calcular el p-valor mediante:

(pValor = 2 * (1 - pt(estadistico, df = n - 1)))
## [1] 0.02637288

¡Fíjate en que hemos multiplicado por dos, al ser un contraste bilateral! Como el p-valor es menor que 0.05, el contraste es significativo al 5% y rechazamos la hipótesis nula a ese nivel de significación (no lo haríamos, por ejemplo, al 1%).

Alternativamente puede hacerse mediante la función t.test:

t.test(cpk, alternative = "two.sided", mu = mu0, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  cpk
## t = -2.4318, df = 17, p-value = 0.02637
## alternative hypothesis: true mean is not equal to 660
## 95 percent confidence interval:
##  349.7697 638.0080
## sample estimates:
## mean of x 
##  493.8889

Apartado (b)

Para hacer el contraste así es necesario asumir que los datos proceden de una población normal.

Cuestión 2.

Un investigador afirma que en la prueba de esfuerzo para la práctica de un determinado deporte, tras un entrenamiento específico, el tiempo transcurrido hasta que se alcanza la frecuencia cardíaca límite, X, es mayor de 30 minutos. Extrae una muestra aleatoria de 9 atletas y obtiene los resultados siguientes: \[ \bar X = 35,\qquad s = 5. 42 \] Suponiendo que la variable X se distribuye de forma normal, compruebe tal hipótesis mediante un contraste, con un nivel de significación del 1%.

Solución:

Haremos un contraste de hipótesis unilateral para la media con hipótesis alternativa \[H_a = \{\mu > \mu_0\}\] donde \(\mu_0 = 30\). La hipótesis nula es, por supuesto, \[H_0 = \{\mu \leq \mu_0\}\] En este caso disponemos directamente de los descriptores muestrales,

n = 9
xbar = 35
s = 5.42

Ahora podemos calcular el valor del estadístico del contraste y el p-valor:

mu0 = 30
(estadistico = (xbar - mu0) / (s/sqrt(n)))
## [1] 2.767528
(pValor = 1 - pt(estadistico, df = n - 1))
## [1] 0.0121948

¡Fíjate en que usamos la cola derecha y no multiplicamos por dos!! Como el p-valor no es menor que 0.01 el contraste no es significativo y rechazamos la hipótesis alternativa.

En este caso no hemos usado t.test porque no disponemos de los datos en bruto. Pero podemos simularlos usando la función mvrnorm de la librería MASS así:

library(MASS)
(tiempos = mvrnorm(n = 9, mu = 35, Sigma = 5.42^2, empirical = TRUE))
##           [,1]
##  [1,] 36.30174
##  [2,] 36.48023
##  [3,] 30.97094
##  [4,] 35.53701
##  [5,] 40.43855
##  [6,] 41.65013
##  [7,] 37.54257
##  [8,] 23.71289
##  [9,] 32.36594

*¡Cuidado con la notación, que es confusa!** En esta función Sigma representa la cuasivarianza muestral. El resultado es un vector de datos que tiene precisamente los descriptores muestrales que necesitamos:

(n = length(tiempos))
## [1] 9
(xbar = mean(tiempos))
## [1] 35
(s = sd(tiempos))
## [1] 5.42

Así que podemos aplicarle la función t.test para hacer el contraste.

t.test(tiempos, alternative = "greater", mu = mu0, conf.level = 0.99)
## 
##  One Sample t-test
## 
## data:  tiempos
## t = 2.7675, df = 8, p-value = 0.01219
## alternative hypothesis: true mean is greater than 30
## 99 percent confidence interval:
##  29.76706      Inf
## sample estimates:
## mean of x 
##        35

El p-valor es naturalmente el mismo. Fíjate en el intervalo de confianza unilateral que se obtiene al usar esta función.

Cuestión 3.

En el ejercicio anterior, ¿podría afirmarse, al nivel de significación del 5%, que la varianza de los tiempos transcurridos hasta que se alcanza la frecuencia cardíaca límite es mayor de 5?

Solución:

Ahora debemos realizar un contraste de hipótesis unilateral para la varianza con hipótesis alternativa \[H_a = \{\sigma^2 > \sigma^2_0\}\] donde \(\sigma^2_0 = 5\). La hipótesis nula es \[H_0 = \{\sigma^2 \leq \sigma^2_0\}\] Vamos a calcular el estadistico del contraste. ¡Ten cuidado! El valor 5 es la cuasivarianza muestral):

sigma02 = 5

(estadistico = (n - 1) * s^2 / sigma02)
## [1] 47.00224

Y el p-valor se obtiene con la cola derecha de la correspondiente distribución \(\chi^2\).

(pValor = 1 - pchisq(estadistico, df = n - 1))
## [1] 1.531876e-07

Alternativamente podemos aprovechar los datos que hemos generado con mvrnorm y usar la librería TeachingDemos, como se explica en el Tutorial07:

library(TeachingDemos) # debes instalarla primero
sigma.test(tiempos, sigmasq = 5, alternative = "greater", conf.level=0.95)
## 
##  One sample Chi-squared test for variance
## 
## data:  tiempos
## X-squared = 47.002, df = 8, p-value = 1.532e-07
## alternative hypothesis: true variance is greater than 5
## 95 percent confidence interval:
##  15.15486      Inf
## sample estimates:
## var of tiempos 
##        29.3764

Cuestión 4.

En una ciudad, una muestra de 150 motoristas ha dado como resultado que 17 de ellos circulaban sin el casco reglamentario. Las autoridades locales consideran que debe iniciarse una campaña de concienciación sobre su uso cuando más del 10% de los motoristas no lo utilice. Se quiere determinar si hay evidencia suficiente, al nivel de significación \(\alpha\) = 0. 05, para poner en marcha la campaña.

Solución:

En este ejercicio vamos a realizar un contraste de hipótesis unilateral para la proporción con hipótesis alternativa \[H_a = \{p > p_0\}\] siendo \(p_0 = 0.1\). La hipótesis nula es \[H_0 = \{p \leq p_0\}\] El estadistico del contraste es:):

p0 = 0.1
q0 = 1 - p0
n = 150
x = 17
(pHat = x / n)
## [1] 0.1133333
(estadistico = (pHat - p0) / sqrt(p0 * q0 / n))
## [1] 0.5443311

y el p-valor es:

(pValor = 1 - pnorm(estadistico))
## [1] 0.2931068

Alternativamente:

prop.test(x = 17, n = 150, alternative = "greater", p = 0.1, conf.level = 0.95, correct = FALSE)
## 
##  1-sample proportions test without continuity correction
## 
## data:  17 out of 150, null probability 0.1
## X-squared = 0.2963, df = 1, p-value = 0.2931
## alternative hypothesis: true p is greater than 0.1
## 95 percent confidence interval:
##  0.07743674 1.00000000
## sample estimates:
##         p 
## 0.1133333

En cualquier caso, el p-valor que hemos obtenido nos lleva a rechazar la hipótesis alternativa.

Cuestión 5.

Se dispone de dos procedimientos A y B para la detección de glucosa en la sangre. Para ver cuál de ellos es más preciso, se han efectuado 10 mediciones con A y 11 con B, de la misma muestra de sangre, obteniéndose los resultados siguientes: \[ s_A^2 = 3.333,\qquad s_B^2= 3.818 \] ¿Se puede afirmar que ambos métodos no son igualmente precisos, al nivel de significación de 5%?

Solución:

En este ejercicio vamos a realizar un contraste de hipótesis bilateral para la varianza con hipótesis alternativa \[H_a = \{\sigma_1^2 \neq \sigma_2^2\}\] La hipótesis nula es \[H_0 = \{\sigma_1^2 = \sigma_2^2\}\] El estadistico del contraste es:):

n1 = 10
n2 = 11
var1 = 3.333
var2 = 3.818
(estadistico = var1 / var2)
## [1] 0.8729701

Fíjate en que puesto que hemos obtenido un estad´sitico menor que 1 debemos usar la cola izquierda de la distribución \(F\) de Fisher (alternativamente, cambiar var1 por var2 y usar la cola derecha). El p-valor es:

(pValor = 2 * pf(estadistico, df1 = n1 - 1, df2 = n2 - 1))
## [1] 0.847916
Alternativamente podemos usar la plantilla del Tutorial-09 que se incluye aquí con los datos ya cumplimentados:

Sea como sea, el p-valor que hemos obtenido nos lleva a rechazar la hipótesis alternativa. La conclusión es que los datos no permiten afirmar que ambos métodos no sean igualmente precisos, al nivel de significación del 95%.

Cuestión 6.

Para comprobar si la tolerancia a la glucosa en sujetos sanos tiende a decrecer con la edad se realizó un test oral de glucosa a dos muestras de pacientes sanos, unos jóvenes y otros adultos. El test consistió en medir el nivel de glucosa en sangre en el momento de la ingestión (nivel basal) de 100 g. de glucosa y a los 60 minutos de la toma. Los resultados fueron los siguientes:

Jovenes:

\(\quad\)
Basal 81 89 80 75 74 97 76 89 83 77
60 minutos 136 150 149 141 138 154 141 155 145 147
\(\quad\)

Adultos:

\(\quad\)
Basal 98 94 93 88 79 90 86 89 81 90
60 minutos 196 190 191 189 159 185 182 190 170 197
\(\quad\)

Al nivel de significación del 2%, conteste a las siguientes cuestiones:

  1. ¿Es mayor la concentración de glucosa en sangre a los 60 minutos, en adultos que en jóvenes?
  2. El contenido basal de glucosa en sangre, ¿es menor en jóvenes que en adultos?

Puedes cargar estos datos desde el siguiente fichero adjunto:

Cuestión 7.

Se ha realizado una prueba clínica con el propanolol para el tratamiento de infarto de miocardio. Se seleccionaron aleatoriamente dos grupos de pacientes afectados por esta enfermedad. De un grupo de 45 pacientes tratados con propanolol, 38 de ellos sobrevivieron los primeros 28 días a partir de su admisión en el hospital. Del segundo grupo, utilizado como grupo control, de 46 pacientes tratados con un placebo, 29 de ellos sobrevivieron a los 28 días.

¿Puede considerarse, al nivel de significación del 5%, que estos resultados constituyen una prueba suficiente de que el propanolol aumenta el índice de supervivencia a los 28 días?

Solución:

En este ejercicio vamos a realizar un contraste de hipótesis unilateral para la diferencia de proporciones con hipótesis alternativa \[H_a = \{p_1 > p_2\}\] siendo \(p_1\) la proporción de pacientes que sobreviven en el grupo tratado con propanolol y \(p_2\) la proporción que sobreviven en el grupo control. La hipótesis nula es \[H_0 = \{p_2 \leq p_1\}\] El estadistico del contraste se calcula así (recuerda que debemos usar la proporción ponderada):):

n1 = 45
pHat1 = 38 / n1
qHat1 = 1 - pHat1
n2 = 46
pHat2 = 29 / n2
qHat2 = 1 - pHat2

pHat = (n1 * pHat1 + n2 * pHat2) / (n1 + n2) # proporción ponderada
qHat = 1 - pHat

(estadistico = (pHat1 - pHat2) / sqrt(pHat * qHat * ((1 / n1) + (1 / n2))))
## [1] 2.316308

y el p-valor es:

(pValor = 1 - pnorm(estadistico))
## [1] 0.01027074

Alternativamente:

prop.test(x = c(38, 29), n = c(45, 46), alternative = "greater", 
          conf.level = 0.95, correct = FALSE)
## 
##  2-sample test for equality of proportions without continuity
##  correction
## 
## data:  c(38, 29) out of c(45, 46)
## X-squared = 5.3653, df = 1, p-value = 0.01027
## alternative hypothesis: greater
## 95 percent confidence interval:
##  0.0670367 1.0000000
## sample estimates:
##    prop 1    prop 2 
## 0.8444444 0.6304348

El p-valor que hemos obtenido nos lleva a rechazar la hipótesis nula. Los datos no permiten afirmar (al 95% de significación) que el porcentaje de supervivencia en el grupo tratado sea menor o igual que en el grupo control.

Cuestión 8.

A nueve jugadores de fútbol sala elegidos aleatoriamente se les ha medido la frecuencia cardíaca al final del primer tiempo y del segundo tiempo de un partido. Los datos obtenidos son:

\(\quad\)
Primer tiempo 177 164 168 173 163 171 170 168 160
Segundo tiempo 160 171 173 172 161 160 169 170 173
\(\quad\)
## [1] 168.2222
## [1] 167.6667

Puedes cargar estos datos desde el siguiente fichero adjunto:

¿Es razonable creer que el esfuerzo físico continuado aumenta la frecuencia cardíaca, con \(\alpha\) = 0. 01?

Solución:
Cuidado: Se trata de contraste de hipótesis para la diferencia de medias en muestras emparejadas. La hipótesis alternativa es \[H_a = \{\mu_2 > \mu_1\}\] siendo \(\mu_2\) la frecuencia cardiaca media al acabar el segundo tiempo y \(\mu_1\) la frecuencia media al acabar el primer tiempo. Por tanto la hipótesis nula es \[H_0 = \{\mu_2\leq\mu_1\}\] Lo hacemos, por ejemplo, con t.test. Como en los intervalos de confianza usamos la opción paired = TRUE Hemos ordenado las muestras de manera que la diferencia es \[\text{tiempo2} - \text{tiempo1}\]

datos = read.table(file = "datos/practica06c08.csv", header = TRUE, sep=",")
tiempo1 = datos$tiempo1
tiempo2 = datos$tiempo2

t.test(tiempo2, tiempo1, alternative = "two.sided", paired = TRUE, conf.level = 0.99)
## 
##  Paired t-test
## 
## data:  tiempo2 and tiempo1
## t = -0.18346, df = 8, p-value = 0.859
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
##  -10.716205   9.605094
## sample estimates:
## mean of the differences 
##              -0.5555556

La conclusión es que rechazamos la hipótesis alternativa. Los datos no permiten afirmar (al 99%) que la frecuencia cardiaca al final del segundo tiempo es mayor que al final del primero.

Corrección:
La función t.test debería haberse usado con greater en este caso.

t.test(tiempo2, tiempo1, alternative = "greater", paired = TRUE, conf.level = 0.99)
## 
##  Paired t-test
## 
## data:  tiempo2 and tiempo1
## t = -0.18346, df = 8, p-value = 0.5705
## alternative hypothesis: true difference in means is greater than 0
## 99 percent confidence interval:
##  -9.326498       Inf
## sample estimates:
## mean of the differences 
##              -0.5555556

Aunque el p-valor cambia, la conclusión es la misma.

Cuestión 9.

De un estudio sobre la incidencia de la hipertensión en una provincia española, se sabe que en la zona rural el porcentaje de hipertensos es del 27.7%. Tras una encuesta a 400 personas de una zona urbana, se obtuvo un 24% de hipertensos.

¿Se puede afirmar al nivel de significación del 5% que el porcentaje de hipertensos en la zona urbana es menor que en la zona rural?

Solución:
Se trata de un contraste de hipótesis para la proporción, con una única muestra. La hipótesis alternativa es \[H_a = \{p > p_0\}\] siendo \(p\) la proporción de hipertensos en la zona urbana y \(p_0 = 0.277\) la proporción de hipertensos en la zona rural. Por tanto la hipótesis nula es \[H_0 = \{p\leq p_0\}\] La manera más rápida de hacerlo es con prop.test (ver la cuestión 4 para métodos alternativos).

prop.test(x = 0.24 * 400, n = 400, alternative = "greater", p = 0.277, 
          conf.level = 0.95, correct = FALSE)
## 
##  1-sample proportions test without continuity correction
## 
## data:  0.24 * 400 out of 400, null probability 0.277
## X-squared = 2.7343, df = 1, p-value = 0.9509
## alternative hypothesis: true p is greater than 0.277
## 95 percent confidence interval:
##  0.206697 1.000000
## sample estimates:
##    p 
## 0.24

Y el p-valor obtenido confirma que los datos no permiten afirmar que la proporción de hipertensos en la zona urbana sea superior al de la zona rural (rechazamos \(H_a\)).

Comentario: Aunque la solución que aparece encima es correcta y el p-valor está bien calculado, la hipótesis alternativa natural para este ejercicio es la contraria. \[H_a = \{p < p_0\}\] En cuyo caso el p-valor se calcula con:

prop.test(x = 0.24 * 400, n = 400, alternative = "less", p = 0.277, 
          conf.level = 0.95, correct = FALSE)
## 
##  1-sample proportions test without continuity correction
## 
## data:  0.24 * 400 out of 400, null probability 0.277
## X-squared = 2.7343, df = 1, p-value = 0.04911
## alternative hypothesis: true p is less than 0.277
## 95 percent confidence interval:
##  0.0000000 0.2767966
## sample estimates:
##    p 
## 0.24

Fíjate en que este p-valor es 1 - (el p-valor de la otra versión).
Pero la conclusión no cambia: rechazamos la nula de este segundo contraste, que es esencialmente la alternativa del otro.

Cuestión 10.

Muchos autores afirman que los pacientes con depresión tienen una función cortical por debajo de lo normal debido a un riego sanguíneo cerebral por debajo de lo normal. A dos muestras de individuos, unos con depresión y otros normales, se les midió un índice que indica el flujo sanguíneo en la materia gris, dado en mg/(100g/min), obteniéndose:

\(\quad\)
Depresivos \(\bar X_1 = 47\) \(s_1 = 7.8\) \(n_1 = 19\)
Normales \(\bar X_2 = 53.8\) \(s_2 = 6.1\) \(n_2 = 22\)
\(\quad\)

¿Hay evidencia significativa a favor de la afirmación de los autores, con \(\alpha\) = 0. 01?

Solución:
Se trata de un contraste de hipótesis para la diferencia de medias, con muestras pequeñas. La hipótesis alternativa es \[H_a = \{\mu_1 < \mu_2\}\] siendo \(\mu_1\) el índice medio en los depresivos y \(\mu_1\) el índice medio en los no depresivos. Por tanto la hipótesis nula es \[H_0 = \{\mu_1\geq \mu_2\}\] Como las muestras son pequeñas para hacer este contraste debemos empezar por contarstar la igualdad de varianzas en ambas poblaciones. La forma más sencilla es usando la plantilla del Capítulo 9 que se adjunta con los datos del ejercicio. El p-Valor de este primer contraste es 0.28, así que rechazamos la alternativa y supondremos que las varianzas son iguales. Para terminar usamos la plantilla del Capítulo 9 adjunta con los datos del contraste de medias:

El p-Valor de este contraste es 0.0016, así que rechazamos la hipótesis nula. Los datos no permiten afirmar que el índice de la población depresiva sea igual o superior al de la no depresiva.

Cuestión 11.

De seis muestras de bilis hepática se ha obtenido por fistulización el pH, con los siguientes resultados: \[ 7.83,\quad 8.52,\quad 7.32,\quad 7.79,\quad 7.57,\quad 6.58 \] Se desea saber, al nivel de significación del 0.05, si la bilis hepática puede considerarse neutra.

Solución:
Se trata de un contraste de hipótesis para la media, con una muestra pequeña. La hipótesis alternativa es \[H_a = \{\mu\neq\mu_0\}\] siendo \(\mu_0 = 7\) el pH neutro y \(\mu\) el pH medio de la bilis hepática. Por tanto la hipótesis nula es \[H_0 = \{\mu =\mu_0\}\] La forma más rápida de hacer el ejercicio es usando t.test así.

muestra = c(7.83,  8.52,  7.32,  7.79,  7.57,  6.58)
t.test(muestra, alternative = "two.sided", mu = 7, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  muestra
## t = 2.2988, df = 5, p-value = 0.06988
## alternative hypothesis: true mean is not equal to 7
## 95 percent confidence interval:
##  6.928867 8.274466
## sample estimates:
## mean of x 
##  7.601667

Y como se deduce del p-valor, los datos no permiten rechazar la hipótesis de que la bilis hepática sea neutra (rechazamos la alternativa).

Cuestión 12.

La prueba de la d-xilosa permite la diferenciación entre una esteatorrea originada por una mala absorción intestinal y la debida a una insuficiencia pancreática, de modo que cifras inferiores a 4 g de d-xilosa, indican una mala absorción intestinal. Se realiza dicha prueba a 10 individuos, obteniéndose una media de 3.5 g y una cuasidesviación típica de 0.5 g ¿Se puede decir que esos pacientes padecen una mala absorción intestinal, al nivel de significación del 5%?

Solución:
Se trata de un contraste de hipótesis para la media, con una muestra pequeña. La hipótesis alternativa es \[H_a = \{\mu < \mu_0\}\] siendo \(\mu_0 = 4\) y \(\mu\) la absorción intestinal media de los pacientes. Por tanto la hipótesis nula es \[H_0 = \{\mu \geq\mu_0\}\] Vamos a aprovechar este ejercicio para practicar el método que consiste en simular los datos de lla muestra cuando no disponemos de ellos usando mvrnorm

library(MASS)
(muestra = mvrnorm(n = 10, mu = 3.5, Sigma = 0.5^2, empirical = TRUE)[,1])
##  [1] 3.323672 4.100371 3.795448 3.789525 3.642982 3.391972 3.800472
##  [8] 2.538307 3.841002 2.776249

Comprobamos los valores muestrales de estos datos simulados:

length(muestra)
## [1] 10
mean(muestra)
## [1] 3.5
sd(muestra)
## [1] 0.5

Como se ve coinciden con lo que queríamos, así que ya podemos hacer el contraste así:

t.test(muestra, alternative = "less", mu = 4, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  muestra
## t = -3.1623, df = 9, p-value = 0.005754
## alternative hypothesis: true mean is less than 4
## 95 percent confidence interval:
##      -Inf 3.789841
## sample estimates:
## mean of x 
##       3.5

El p-valor nos permite rechazar la hipótesis nula y por tanto concluir que los datos no permiten afirmar que la absorción intestinal media de los pacientes sea mayor o igual a 4.

Cuestión 13.

A los pacientes internados en un hospital traumatológico se les aplica un método de fisioterapia con el que el 70% de ellos requiere posteriormente algún tipo de intervención quirúrgica. Para determinar si un nuevo método de fisioterapia reduce el porcentaje de intervenciones posteriores, se aplica este a 30 pacientes de los cuales 17 necesitaron alguna intervención quirúrgica después. Compruebe que no hay razones suficientes para afirmar la eficacia del método con un nivel de significación del 5%

Solución:
Se trata de un contraste de hipótesis para la proporción. Para hacer este ejercicio vamos a tratar la muestra como grande, aunque ya hemos dicho que en general es mejor usar \(100\) como límite de tamaño para considerar las muestras grandes. La hipótesis alternativa es \[H_a = \{p < p_0\}\] siendo \(p_0 = 0.7\) la proporción de los que necesitan cirugía con el método previo y \(p\) la de los que la necesitan con el nuevo método. Por tanto la hipótesis nula es \[H_0 = \{p \geq p_0\}\] La forma más rápida de hacer el ejercicio es usando prop.test así.

prop.test(x = 17, n = 30, alternative = "less", p = 0.7, conf.level = 0.95)
## 
##  1-sample proportions test with continuity correction
## 
## data:  17 out of 30, null probability 0.7
## X-squared = 1.9444, df = 1, p-value = 0.08159
## alternative hypothesis: true p is less than 0.7
## 95 percent confidence interval:
##  0.0000000 0.7184049
## sample estimates:
##         p 
## 0.5666667

Solución:
Como se ve, rechazamos la alternativa: los datos no permiten afirmar que el nuevo tratamiento sea más eficaz que el anterior.

Cuestión 14.

Para estudiar si los problemas cardiovasculares son hereditarios se mide el nivel de colesterol en un grupo de 100 niños entre 2 y 14 años, hijos de hombres que sufren de colesterol alto, obteniéndose los valores siguientes: \[\bar X_1 = 207.3 mg/dL,\qquad s_1 = 2 mg/dL.\] Para un grupo de control de 74 niños elegidos de familias sin historial de problemas del corazón, \[\bar X_2 = 193.4 mg/dL, s_2 = 17.3 mg/dL.\] ¿Se puede considerar, con un nivel de significación del 5%, que los hijos de padres con nivel colesterol alto tienen mayor nivel de colesterol que el resto de la población?

Indicación:
Es muy parecido al 10.

Cuestión 15.

La tabla siguiente muestra los efectos de un placebo y de la hidroclorotiacida sobre la presión sanguínea sistólica de 11 pacientes.

\(\quad\)
Placebo 211 210 210 203 196 190 191 177 173 170 163
H - cloro 181 172 196 191 167 161 178 160 149 119 156
\(\quad\)

¿Se puede afirmar, al nivel del 5%, que existe diferencia en la presión sistólica media durante la utilización de estos dos fármacos?

Puedes cargar estos datos desde el siguiente fichero adjunto:

Indicación:
Es parecido al 8, se trata de un contraste de medias en muestras emparejadas.