A continuación se presenta una serie de ejercicios que permitiran interiorizar los conocimientos aprendidos en relación a pruebas de hipótesis desde R. El documento contiene la solución de 6 ejercicios que hacen énfasis en:
El desarrollo del taller hace parte del contenido académico elaborado por el profesor Aquiles Enrique Darghan Contreras para la materia de de Métodos Multivariados.
Previo a la solución de los ejercicios es necesario instalar y cargar las siguientes librerías y funciones diseñadas para el análisis.
library(dplyr)
library(kableExtra)
library(ggplot2)
library(ICSNP)
library(biotools)
## ---
## biotools version 3.1
library(rrcov)
library(ergm)
pajust <- function(hipotesis,alfa){
pvalor = 1-((1-alfa)^(1/hipotesis))
return(pvalor)
}
Una línea base de longitud calibrada en 955m se midió 10 veces. Cada medida es independiente y se hizo con la misma precisión.
Los datos de las mediciones son:
M1 | M2 | M3 | M4 | M5 | M6 | M7 | M8 | M9 | M10 |
---|---|---|---|---|---|---|---|---|---|
0 | 955.1 | 954.8 | 956.2 | 956.4 | 957.2 | 958.1 | 960.4 | 953.2 | 954.8 |
med=c(965.1,955.1,954.8,956.2,956.4,957.2,958.1,960.4,953.2,954.8)
pruebat <- t.test(med,alternative = "t",mu = 0,var.equal = F, conf.level = 0.90);
pruebat #La media está en la cola inferior
##
## One Sample t-test
##
## data: med
## t = 877.89, df = 9, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 90 percent confidence interval:
## 955.1314 959.1286
## sample estimates:
## mean of x
## 957.13
set.seed(2020)
B = 100 #repeticiones Bootstrap
df=data.frame(med)
remuestreo <- with(df,matrix(sample(med,size = 10,replace = TRUE), B, 10))
dim(remuestreo)
## [1] 100 10
####### las medias de las submuestras
medias <- apply(remuestreo, 1, mean)
pruebat <- t.test(medias,alternative = "t",mu = 955,var.equal = F, conf.level = 0.90);
pruebat
##
## One Sample t-test
##
## data: medias
## t = 13.232, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 955
## 90 percent confidence interval:
## 959.3376 960.5824
## sample estimates:
## mean of x
## 959.96
Conclusión: Los resultados consideran que la media no es 955 sino que está dentro en el rango 959.3-960.5
Si los datos son mediciones, la primera prueba no tendría sentido porque solo es una media. En este caso la primera prueba estaría asumiendo que cada medición es una valor promedio diferente y no una medición. Para que el error estándar Bootstrap tenga sentido deberíamos tener más repeticiones y valores promedios de diferentes conjuntos de datos.
Se midió un ángulo en seis series en una condición atmosférica particular utilizando un instrumento especializado, pero usando dos operarios simultáneamente. Los datos de las medidas angulares fueron:
Operario X | 112°47’34’’ | 112°44’39’’ | 113°01’34’’ | 112°57’30’’ | 113°00’14’’ | 111°59’58’’ | 112°09’58’’ |
Operario Y | 113°47’24’’ | 112°14’39’’ | 112°01’34’’ | 112°17’32’’ | 113°10’04’’ | 112°19’18’’ | 112°59’48’ |
#Crear el dataframe con la información de operarios
set.seed(2020);options(digits = 9)
opx<-c(112.7927778, 112.7441667, 113.0261111, 112.9583333, 113.0038889, 111.9994444, 112.1661111)
opy<-c(113.79,112.2441667,112.0261111,112.2922222,113.1677778,112.3216667,112.9966667)
data=c(opx,opy)
operario <- gl(2, 7, 14, labels = c("operariox", "operarioy"))
df_oper <- data.frame(data,operario)
Definición de la hipótesis nula:
\[h_o: \bar{x}_{operario} = \bar{y}_{operario}\]
####### estadisticas descriptivas por grupo
grupos <- group_by(df_oper, operario)
summarise(grupos, media = mean(data, na.rm = T), desv = sd(data,na.rm=T),cv = desv*100/media, muestra = length(data))
## `summarise()` ungrouping output (override with `.groups` argument)
#Crear Datos por operario
dataX <- c(filter(df_oper, operario == "operariox"));
dataY <- c(filter(df_oper, operario == "operarioy"));
n.datos <- with(df_oper, summary(operario)) #cada tamano de muestra
B = 1000 #repeticiones Bootstrap
#las submuestras de cada metodo
met.x <- with(df_oper,matrix(sample(opx,size = n.datos[1]*B,replace = TRUE), B, n.datos[1]))
met.y <- with(df_oper,matrix(sample(opy,size = n.datos[2]*B,replace = TRUE), B, n.datos[2]))
dim(met.x)
## [1] 1000 7
dim(met.y)
## [1] 1000 7
####### las medias de las submuestras
medias.x <- apply(met.x, 1, mean)
medias.y <- apply(met.y, 1, mean)
######## el vector de diferencia de medias
stat.boot <- medias.x - medias.y
length(stat.boot)
## [1] 1000
ggplot(data.frame(x = stat.boot), aes(x = x)) + geom_density()+
labs(x = "diferencia de medias")+labs(y = "densidad")+
labs(title = "Distribución de la diferencia de medias",
subtitle = "Comparación de los dos operarios")+
labs(caption = "(datos suministrados en clase)")+
geom_vline(xintercept =0,col="azure4")
eebdm <- sd(stat.boot); eebdm # error bootstrap
## [1] 0.255560199
pruebat <- t.test(medias.x,medias.y, alternative = "two.sided", mu=0, conf.level = 0.95);
eec <- c((mean(medias.x)-mean(medias.y))/pruebat$statistic);eec
## t
## 0.00830908668
veces <- eebdm/eec;veces
## t
## 30.7567135
# Cálculo del estadastico t bootstrap
tboot <- abs((mean(medias.x)-mean(medias.y))/eebdm)
# Obtener el p-valor del estadístico
pvalor_boot <- pt(tboot,n.datos[1]+n.datos[2]-2,lower.tail = F)
pvalor_boot
## [1] 0.479290842
pvalor_usual <- pruebat$p.value
pvalor_usual
## [1] 0.10307391
curve(dt(x, 46), from = -5, to = 5, col = "orange",xlab = "cuantil t", ylab = "densidad", lwd = 2)
legend("topleft", legend = paste0("DF = ", 46),col = "sky blue",lty = 1, lwd = 2)
tt0025 <- qt(0.025,46,lower.tail = FALSE)
tt0975 <- qt(0.025,46,lower.tail = TRUE)
abline(h = 0, col = "darkgoldenrod2")
segments(tt0025, 0, tt0025, dt(tt0025,df = 46))
segments(-tt0025, 0, -tt0025, dt(tt0025,df = 46))
qbonf <- qt(pajust(2, 0.05)/2, 46, lower.tail = F)
segments(qbonf, 0, qbonf, dt(qbonf, df = 46), col = "darkblue")
segments(-qbonf, 0, -qbonf, dt(qbonf, df = 46), col = "darkblue")
Suponga del problema anterior que se seleccionó como mejor operario aquel que generó el menor coeficiente de variación de las mediciones angulares y se utilizó en el siguiente estudio. Ahora se generaron las mismas 6 medidas, pero una a las 10 de la mañana y otra a las 12 del mediodía, cuando la temperatura se había incrementado 5° C. Los datos se muestran a continuación:
10:00 am | 110°17’14’’ | 110°41’19’’ | 110°11’14’’ | 111°00’30’’ | 110°08’24’’ | 110°19’38’’ | 110°09’58’’ |
12:00 m | 111°17’24’’ | 110°44’39’’ | 110°21’11’’ | 111°17’32’’ | 113010’14’’ | 110°29’28’’ | 110°20’08’’ |
Determinar al 95% de nivel de confianza si se incrementó la medida angular en las dos horas registradas. Haga una representación gráfica radial para ilustrar el comportamiento de ambas medidas.
am=c(110.2872222,110.6886111,110.1872222,111.0083333,110.14,110.3272222,110.1661111)
m=c(111.29,110.7441667,110.3530556,111.2922222,113.1705556,110.4911111,110.3355556)
data=c(am,m)
hora <- gl(2, 7, 14, labels = c("10:00 am", "12:00 m"))
temp <- gl(2, 7, 14, labels = c(10, 15))
df_hora <- data.frame(data,temp,hora)
n.datos <- with(df_hora, summary(hora)) #cada tamano de muestra
B = 1000 #repeticiones Bootstrap
met.10 <- with(df_hora,matrix(sample(am,size = n.datos[1]*B,replace = TRUE), B, n.datos[1]))
met.12 <- with(df_oper,matrix(sample(m,size = n.datos[2]*B,replace = TRUE), B, n.datos[2]))
dim(met.10)
## [1] 1000 7
dim(met.12)
## [1] 1000 7
####### las medias de las submuestras
medias.10 <- apply(met.10, 1, mean)
medias.12 <- apply(met.12, 1, mean)
######## el vector de diferencia de medias
stat.boot <- medias.12 - medias.10
length(stat.boot)
## [1] 1000
ggplot(data.frame(x = stat.boot), aes(x = x)) + geom_density()+ coord_polar() + labs(x = "diferencia de medias")+labs(y = "densidad") + labs(title = "Distribucion de la diferencia de medias",subtitle = "Comparacion de los dos horarios")+labs(caption = "(datos suministrados en clase)")
pruebat<-t.test(medias.12,medias.10, alternative = "greater", mu=0, conf.level = 0.95, pared=TRUE)
pruebat
##
## Welch Two Sample t-test
##
## data: medias.12 and medias.10
## t = 61.54576, df = 1194.691, p-value < 2.22e-16
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 0.687272497 Inf
## sample estimates:
## mean of x mean of y
## 111.100554 110.394395
Construya una función basada en la distribución Uniforme[0,1] para que genere medidas angulares entre 0°0’0’’ y 360°0’0’’. Construya un histograma circular generando una muestra de 500 medidas angulares.
angulalet <- function(nm){
x<-runif(nm,0,360)
return(x)
}
sm=angulalet(500)
ggplot(data.frame(x = sm), aes(x = x)) + geom_density()+ coord_polar() + labs(x = "Medidas angulares aleatorias")+labs(y = "Densidad")+ labs(title = "Distribucion de la medidas angulares aleatorias")
Una línea de base se observa repetidamente usando un instrumento EDM durante un período de tiempo. Cada día se hacen 10 observaciones y se promedian. Las varianzas de las observaciones se enumeran a continuación.
Dia | 1 | 2 | 3 | 4 | 5 | 0 | 0 |
Varianza | 50 | 61 | 51 | 53 | 50 | 0 | 0 |
A un nivel de confianza del 95%, ¿los resultados del día 2 son estadísticamente diferentes de los del día 5?
var_d2=61
var_d5=50
coc=var_d2/var_d5;
#Prueba F
ftest=qf(0.05,9,9,lower.tail = F)
if(coc > ftest){
result="Se rechaza la hipotesis nula de igualdad de varianzas";result
} else {
result="No se puede rechaza la hipotesis nula de igualdad de varianzas";result
}
## [1] "No se puede rechaza la hipotesis nula de igualdad de varianzas"
Se están probando dos dispositivos que permiten medir el ángulo recto y la distancia a una serie de puntos en un transecto. La figura muestra la situación:
Toma de datos
La siguiente tabla muestra los datos de toma:
Datos
Por lo anterior se procede a cargar los datos que se encuentran en la tabla mostrada, haciendo acotación de que los valores angulares deben ser convertidos a valores decimales:
distancia=c(100.02,200.12,300.08,399.96,419.94,519.99,620.04,720.08,100.06,199.45,298.08,398.96,420.15,520.02,621.01,721.11)
angulo=c(90.004444444, 90.002222222, 89.996666667, 90.017222222, 89.996666667, 90.005555556, 89.993333333, 90.001666667, 90.035833333, 90.021111111, 90.988333333, 89.017222222, 89.835555556, 89.038888889, 89.165555556, 89.029444444)
disp <- gl(2, 8, 16, labels = c("A", "B"))
datos5=data.frame(distancia,angulo,disp)
datos5a=datos5[1:8,]
datos5b=datos5[9:16,]
A partir de los datos presentados, se proponen la realización de los siguientes ejercicios:
pda1=t.test(datos5a$distancia,datos5b$distancia, alternative = "two.sided", var.equal = T, mu=0, conf.level = 0.95)
pda1
##
## Two Sample t-test
##
## data: datos5a$distancia and datos5b$distancia
## t = 0.001662933, df = 14, p-value = 0.998697
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -223.922263 224.269763
## sample estimates:
## mean of x mean of y
## 410.02875 409.85500
pda2=t.test(datos5a$angulo,datos5b$angulo, alternative = "two.sided", var.equal = T, mu=0, conf.level = 0.95)
pda2
##
## Two Sample t-test
##
## data: datos5a$angulo and datos5b$angulo
## t = 1.441763, df = 14, p-value = 0.171364
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.175896373 0.897354707
## sample estimates:
## mean of x mean of y
## 90.0022222 89.6414931
Para los dos casos no se puede rechazar la hipotesis nula de igualdad de medias, es decir que las medias de los dos grupos pueden asumirse como iguales.
pajust <- function(hipotesis,alfa){
pvalor = 1-((1-alfa)^(1/hipotesis))
return(pvalor)
}
crrb=pajust(2, 0.05) ## Cálculo del nuevo nivel de significancia por Bonferroni
crrb
## [1] 0.0253205655
crrb=1-crrb ## Nivel de Confianza con Bonferroni
crrb
## [1] 0.974679434
A pesar de que es más exigente por corrección de Bonferroni (nivel de confianza de 97.5%), aún no se puede rechazar las hípotesis nulas debido a que los p-valor de las pruebas son mayores a 0.10, es decir mucho más que el nivel de significancia.
cor(datos5$distancia,datos5$angulo)
## [1] -0.448013061
Según el índice de correlación existe una leve relación lineal negativa entre las dos variables.
n.datos <- with(datos5, summary(disp)) #cada tamaño de muestra
B = 1000 #repeticiones Bootstrap
#las submuestras de cada método
met.a.dist <- with(datos5a,
matrix(sample(datos5a$distancia,
size = n.datos[1]*B,
replace = TRUE), B, n.datos[1]))
met.a.ang <- with(datos5a,
matrix(sample(datos5a$angulo,
size = n.datos[1]*B,
replace = TRUE), B, n.datos[1]))
met.b.dist <- with(datos5b,
matrix(sample(datos5b$distancia,
size = n.datos[1]*B,
replace = TRUE), B, n.datos[2]))
met.b.ang <- with(datos5b,
matrix(sample(datos5b$angulo,
size = n.datos[1]*B,
replace = TRUE), B, n.datos[2]))
medias.a.dist<- apply(met.a.dist, 1, mean)
medias.a.ang<- apply(met.a.ang, 1, mean)
medias.b.dist <- apply(met.b.dist, 1, mean)
medias.b.ang<- apply(met.b.ang, 1, mean)
pdistd<-t.test(medias.a.dist,medias.b.dist,alternative = "two.sided",0, conf.level = 0.95, var.equal = F)
pdistd
##
## Welch Two Sample t-test
##
## data: medias.a.dist and medias.b.dist
## t = -0.07360255, df = 1997.76, p-value = 0.941334
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.1262375 5.6830325
## sample estimates:
## mean of x mean of y
## 411.796494 412.018096
pdang<-t.test(medias.a.ang,medias.b.ang,alternative = "two.sided", mu=0, conf.level = 0.95)
pdang
##
## Welch Two Sample t-test
##
## data: medias.a.ang and medias.b.ang
## t = 46.54637, df = 999.2003, p-value < 2.22e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.336422318 0.366037196
## sample estimates:
## mean of x mean of y
## 90.0021763 89.6509466
Con Bootstrap se observa que se rechaza la hipótesis de igualdad de medias para la variable de ángulos, mientras que en distancia persiste la misma conclusión.
Por corrección por Bonferroni el nivel de significancia para cada prueba es de 2.5%, sin embargo por los p-valor persiste la conclusión manifestada previamente.
\(H_0 : [\mu_{d_A}, \mu_{a_A}] = [\mu_{d_B}, \mu_{a_B}]\)
Z=as.matrix(datos5[,1:2])
g=as.matrix(datos5[,3])
PruebaH2 <- HotellingsT2(Z~g, mu = c(0,0)) ##
PruebaH2
##
## Hotelling's two sample T2-test
##
## data: Z by g
## T.2 = 1.255882, df1 = 2, df2 = 13, p-value = 0.317201
## alternative hypothesis: true location difference is not equal to c(0,0)
De acuerdo a la prueba T\(^2\) de Hotelling no se puede rechazar la hípotesis de igualdad de medias de los dos grupos para las dos variables de interés desde un enfoque bivariado.
## Matriz de varianzas y covarianzas metodo A
data.a <- as.matrix((filter(datos5[, 1:2], disp == "A")), ncol = 2)
sa <- var(data.a)
sa
## distancia angulo
## distancia 43540.489698214 -2.79522189e-01
## angulo -0.279522189 5.45855368e-05
## Matriz de varianzas y covarianzas metodo B
data.b <- as.matrix((filter(datos5[, 1:2], disp == "B")), ncol = 2)
sb <- var(data.b)
sb
## distancia angulo
## distancia 43794.854429 -100.166801766
## angulo -100.166802 0.500745567
## Prueba M de Box
boxM(datos5[,1:2], datos5[,3])
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: datos5[, 1:2]
## Chi-Sq (approx.) = 46.49637, df = 3, p-value = 4.44783e-10
De acuerdo a la prueba M de Box no se puede asumir que las matrices de varianzas y covarianzas son iguales, por lo anterior se buscaron otras funciones de R que de las librerias ergm y rrcov que permiten realizar la prueba T\(^2\) de Hotelling indicando este aspecto.
approx.hotelling.diff.test(data.a[,-3], data.b[,-3], mu0 = 0, assume.indep = T, var.equal = F)
##
## Hotelling's Two -Sample T^2-Test
##
## data:
## T^2 = 2.704977, param = 2.00000, df = 12.65632, p-value = 0.323394
## alternative hypothesis: two.sided
## null values:
## distancia angulo
## 0 0
## sample estimates:
## distancia angulo
## 0.173750000 0.360729167
T2.test(data.a[,-3], data.b[,-3], method = "c")
##
## Two-sample Hotelling test
##
## data: data.a[, -3] and data.b[, -3]
## T2 = 2.704977, F = 1.255882, df1 = 2, df2 = 13, p-value = 0.317201
## alternative hypothesis: true difference in mean vectors is not equal to (0,0)
## sample estimates:
## distancia angulo
## mean x-vector 410.02875 90.0022222
## mean y-vector 409.85500 89.6414931
Obsérvese que cuando se indica que la matriz de covarianza no es igual, aumenta el p-valor, no obstante se mantiene la misma conclusión sobre la existencia de igualdad del vector de medias entre los dos grupos de datos.
ggplot(data=datos5, aes(x=angulo, y=distancia, col=disp, group=disp)) + geom_point() + coord_polar()
Debido a que no se pudo rechazar la \(H_o\) de igualdad de medias, se realiza un modelo lineal entre la distancia y el ángulo:
model1=lm(datos5$distancia ~ datos5$angulo)
summary(model1)
##
## Call:
## lm(formula = datos5$distancia ~ datos5$angulo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -278.0460 -157.5174 18.3621 107.8981 341.5290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16090.9629 8363.2881 1.92400 0.074927 .
## datos5$angulo -174.5791 93.1083 -1.87501 0.081805 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 186.823 on 14 degrees of freedom
## Multiple R-squared: 0.200716, Adjusted R-squared: 0.143624
## F-statistic: 3.51567 on 1 and 14 DF, p-value: 0.0818049
De acuerdo de la prueba t para cada variable y la prueba F del modelo de regresión, el modelo posee significancia para un nivel de confianza del 95%, pero se acerca bastante al valor permitido, además los valores de R\(^2\) indican que el modelo no posee algún grado de significado, por la proporción de la varianza explicado.
Por lo anterior, se optó realizar el segundo modelo propuesto (inverso de la distancia)
model2=lm(1/datos5$distancia ~ datos5$angulo)
summary(model2)
##
## Call:
## lm(formula = 1/datos5$distancia ~ datos5$angulo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.002420437 -0.001537940 -0.000851183 0.000508107 0.006184345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.14146822 0.12201305 -1.15945 0.26567
## datos5$angulo 0.00161416 0.00135837 1.18831 0.25448
##
## Residual standard error: 0.00272559 on 14 degrees of freedom
## Multiple R-squared: 0.0916218, Adjusted R-squared: 0.0267376
## F-statistic: 1.41208 on 1 and 14 DF, p-value: 0.254479
El modelo mejora considerablemente, pero aún no es del todo adecuado.
Se calculan los valores ponderados por el inverso de la distancia:
datos5a$inversdis=1/datos5a$distancia
datos5a$angulopond=datos5a$inversdis*datos5a$angulo
datos5b$inversdis=1/datos5b$distancia
datos5b$angulopond=datos5b$inversdis*datos5b$angulo
cor(datos5b$angulopond,datos5a$angulopond)
## [1] 0.999956602
mean(datos5a$angulopond)
## [1] 0.316513643
mean(datos5a$angulo)
## [1] 90.0022222
mean(datos5b$angulopond)
## [1] 0.316452432
mean(datos5b$angulo)
## [1] 89.6414931
Se observa que los valores ponderados por el inverso de la distancia tienen una correlación lineal cercana a 1. Además, las medias ponderadas para los dos grupos son casi idénticas, siendo esta conclusión un poco más confusa si se realizará a partir de las medias usuales de solamente el ángulo, puesto que la ponderada posee implicitamente la contribución de la distancia.
Richard A. Johnson & Dean W. Wichern (2014). Inferences about the mean vector. In: Richard A. Johnson & Dean W. Applied Multivariate Statistical Analysis. Pearson Prentice Hall. 775 pp.
Douglas C. Montgomery (2004). Design and Analysis of Experiments. Editorial Limusa S. A. 692 pp.