Caso en el que x es una variable categórica

Teoría

Recordemos que, cuando suponemos que $x$ es una variable categórica de dos niveles, el objetivo es comparar la media de dos poblaciones. En este caso consideraremos la paremetrización

\[ z=\mathbb{I}(x= cat \ \ 1) \]

de donde la linealidad se cumple por la forma de parametrización de las esperanzas

\[ \mathbb{E}[y; x=cat\ \ 1 ]=\mathbb{E}[z=1]=\beta_{0}+\beta_{1} \ \ \textrm{y}\ \ \mathbb{V}[y; x=cat\ \ 1 ]=\sigma^{2} \] y también

\[ \mathbb{E}[y; x=cat\ \ 2 ]=\mathbb{E}[z=0]=\beta_{0} \ \ \textrm{y}\ \ \mathbb{V}[y; x=cat\ \ 2 ]=\sigma^{2} \]

Dado lo anterior restaría revisar

Homocedasticidad.
Normalidad.
Independencia de los errores/observaciones.

Luego, en vez de revisar un scatterplot inicial como se hacía en el caso de $x$ continua, cuando $x$ es categórica se recomienda usar un boxplot de $y$ para cada categoría de $x$. Así como complementar con lo siguiente

Homocedasticidad.

Los boxplots de $y$ en cada categoría se espera que presenten un rango intercuantílico similar.
Se pueden realizar pruebas específicas sobre igualdad de varianzas de $y$ en los grupos, es decir,

\[ H_{0}: \sigma^{2}_{1}=\sigma^{2}_{2} \ \ vs \ \ H_{a}: \sigma^{2}_{1}\neq \sigma^{2}_{2} \]

por ejemplo: Barlett, Levene o Fligner-Killeen. Se espera no rechazar $H_{0}$.

Normalidad e independencia de las observaciones

Aplicar mismas pruebas que en el caso continuo, pero en cada grupo y directamente sobre $y$.

En muchos problemas de este tipo, por diseño se cumple la aleatoriedad, pero se sugiere hacer la revisión del proceso de generación de los datos.

Cuando hay evidencia en contra de la normalidad o de la homocedasticidad se sugiere lo siguiente:

Realizar transformaciones a $y$, por ejemplo utilizando las de tipo Box-Cox.
No hay homocedasticidad pero si normalidad en cada grupo. Procedemos a
- Usar una prueba $t$ para comparar las medias y que permita diferente varianza:
  - t.test(...,var.equal=FALSE)
  - oneway.test(...,var.equal=FALSE)
- Usar un modelo de regresión con la matriz de varianza que permita variar los valores de la diagonal por grupos.
  - gls(weights=varIdent(form= ~ 1 | group)) en nlme
No hay normalidad.
- Si hay homocedasticidad, utilizar el método bootstrap.
- Pruebas no paramétricas (Mann-Whitney test o Kruskal-Wallis test)

Ejemplo

Comenzamos por cargar los datos

# AADT: Average Annual Daily Traffic Data
# the response is taken as AADT, which is the annual average of the number
# of vehicles that pass through a road each day

# Sólo noninterstate. Rural vs Urban  

library(AID)
data(AADT)
Datos=AADT

Notamos a partir de

str(Datos)

## 'data.frame':    121 obs. of  8 variables:
##  $ aadt   : int  1616 1329 3933 3786 465 794 618 1150 1538 1769 ...
##  $ ctypop : int  13404 52314 30982 25207 20594 11507 9379 24991 30982 24991 ...
##  $ lanes  : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ width  : int  52 60 57 64 40 44 43 42 59 41 ...
##  $ control: Factor w/ 2 levels "access control",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ class  : Factor w/ 4 levels "rural interstate",..: 2 2 4 4 2 2 2 2 2 2 ...
##  $ truck  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ locale : Factor w/ 3 levels "rural","urban, population <= 50,000",..: 1 1 2 2 1 1 1 1 1 1 ...

que nos interesa trabajar sobre `.classque contiene cuatro niveles

levels(Datos$class)

## [1] "rural interstate"    "rural noninterstate" "urban interstate"   
## [4] "urban noninterstate"

donde dos de ellos son los que nos interesan (rural noninterstate y urban noninterstate). De tal manera procedemos a filtrar los datos

# Cargamos la librería necesaria
library(tidyverse)

# Filtramos del dataframe Datos, los elementos de la columna class para 
# rural noninterstate y urban noninterstate

Datosfil=Datos %>% filter (class %in% c( "rural noninterstate","urban noninterstate"))

# Después de realizar un filtro a una variable categórica se recomienda
# ver los niveles nuevamente 
levels(Datosfil$class)

## [1] "rural interstate"    "rural noninterstate" "urban interstate"   
## [4] "urban noninterstate"

De tal manera, a pesar de tener en el dataframe Datosfil los niveles rural noninterstate y urban noninterstate, de acuerdo a la función levels() continuamos teniendo cuatro niveles. Para eliminar los niveles que no nos interesan utilizaremos el código Datosfil$class=droplevels(Datosfil$class) con el cual conseguimos que los niveles que no están en el dataframe Datosfil sean eliminados

Datosfil$class=droplevels(Datosfil$class)

# Corroboramos
levels(Datosfil$class)

## [1] "rural noninterstate" "urban noninterstate"

Luego, por comodidad cambiamos el nombre de dichos niveles (pues están muy largos)

levels(Datosfil$class) <- list(RNI = "rural noninterstate", UNI = "urban noninterstate")
levels(Datosfil$class)

## [1] "RNI" "UNI"

Boxplot

Considerando boxplot(aadt ~ class, data = Datosfil, col = "white", outline=FALSE) nos permite contrastar a la variable $y$ (aadt) contra los niveles que tenga la variable class (de tipo factor); la parte de col="white" se utiliza para indicar a R que no pinte los fondos de los boxplots. Con lo anterior conseguimos

boxplot(aadt ~ class, data = Datosfil, col = "white", outline=FALSE)

pero, como en el caso de los scatterplot, nos gustaría ver los puntos referentes a los datos. Para ello utilizamos stripchart() como sigue

boxplot(aadt ~ class, data = Datosfil, col = "white", outline=FALSE)
stripchart(aadt ~ class, data = Datosfil,
           method = "jitter",
           pch = 19,
           col = 2:4,
           vertical = TRUE,
           add = TRUE)

donde utilizando method="jitter" agregando cierto ruido aleatorio a los datos, con lo cual conseguimos tener mejor visualización de los puntos. Sin este método los puntos graficados se verían como

boxplot(aadt ~ class, data = Datosfil, col = "white", outline=FALSE)
stripchart(aadt ~ class, data = Datosfil,
           pch = 19,
           col = 2:4,
           vertical = TRUE,
           add = TRUE)

Retomando el primer boxplot con los puntos, podemos ver que el rango intercuantílico del boxplot del nivel UNI es mayor al del nivel RNI; parece que enUNI hay mucha más variabilidad.

Observación: El rango intercuartílico nos indicará la variabilidad de los datos, mientras más grande más variabilidad e inversamente. También, buscamos una nube de puntos en los boxplots de la misma longitud respecto a la variable $y$ para que se pueda cumplir el supuesto de la homocedasticidad. Buscamos lo anterior pues recordemos que en los boxplots de $y$ en cada categoría se espera que presenten un rango intercuantílico similar así como una misma distribución de los puntos.

Así, parece verse que el supuesto de la homocedasticidad no se está cumpliendo

Luego, buscamos simétria en los boxplots, es decir que la mediana se encuentre cercana al centro del rango intercuantílico, para ver si el supuesto de la normalidad se está cumpliendo. Podemos ver que en ambos boxplots no hay simetría (aunque en el segundo boxplot tal vez si pueda darse) por lo que el supuesto de la normalidad parece no cumplirse.

Ajuste del modelo

De manera muy similar al caso en que $x$ es continua, para realizar el ajuste del modelo usamos

fit <- lm(aadt ~ class, data=Datosfil)

# y damos un rápido vistazo
summary(fit)

## 
## Call:
## lm(formula = aadt ~ class, data = Datosfil)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15955  -2939  -1849   1953  61122 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     3465       1320   2.625   0.0101 *  
## classUNI       13756       2060   6.679 1.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9875 on 93 degrees of freedom
## Multiple R-squared:  0.3242, Adjusted R-squared:  0.3169 
## F-statistic: 44.61 on 1 and 93 DF,  p-value: 1.732e-09

donde classUNI es la variable que R crea en automático para referirse a la variable $z$ que vimos en la parte teórica, donde nos explica que $z=1$ cuando $x=UNI$.

Cálculo “a mano”

Para ello

# Realizamos un filtro en Datosfil sólo de las columnas aadt y class
Datosfil=Datosfil %>% select(aadt, class)

# Agregamos una columna en donde asignamos el valor de 1 a aquellos datos que tengan
# "UNI" en la columna class
Datosfil$zclassUNI= (Datosfil$class=="UNI")*1

str(Datosfil)

## 'data.frame':    95 obs. of  3 variables:
##  $ aadt     : int  1616 1329 3933 3786 465 794 618 1150 1538 1769 ...
##  $ class    : Factor w/ 2 levels "RNI","UNI": 1 1 2 2 1 1 1 1 1 1 ...
##  $ zclassUNI: num  0 0 1 1 0 0 0 0 0 0 ...

cabe resaltar que Datosfil$class=="UNI" nos devolverá una columna con TRUE y FALSE, así, al multiplicar por 1(Datosfil$class=="UNI")*1 convertimos los TRUE en 1 y los FALSE en 0. De tal manera conseguimos así construir la variable $z$. Además, debemos de corroborar que la variable zclassUNI sólo tenga valores 0 y 1 para posteriormente realizar el ajuste del modelo. Luego realizamos el ajuste del modelos

fitB <- lm(aadt ~ zclassUNI, data=Datosfil)

Ahora, deseamos comparar los valores obtenidos en ambos modelos por lo que utilizamos

library(jtools)
export_summs(fit, fitB, scale = FALSE)

	Model 1	Model 2
(Intercept)	3464.55 *	3464.55 *
	(1319.59)	(1319.59)
classUNI	13756.47 ***
	(2059.53)
zclassUNI		13756.47 ***
		(2059.53)
N	95	95
R2	0.32	0.32
* p < 0.001; p < 0.01; * p < 0.05.

en los cuales podemos observar exactamente los mismos resultados.

Después, tenemos que $\hat{\beta_{0}}=3464.55$ y $\hat{\beta_{1}}=13756.47$ por lo que

$\mathbb{E}[aadt; class=RNI]=\hat{\beta_{0}}=3464.55$ pues en este caso tenemos el equivalente a $\mathbb{E}[y=aadt; x=RNI]=\mathbb{E}[z=0]=\mu_{1}$.
También

\[ \mathbb{E}[aadt; class=UNI]=\hat{\beta_{0}}+\hat{\beta_{0}}=3464.55+13756.47=17221.02 \]

pues en este caso tenemos el equivalente a $\mathbb{E}[y=aadt; x=UNI]=\mathbb{E}[z=1]=\mu_{2}$.

Otra forma de ajustar

Alternativamente podemos realizar el ajuste del modelo definiendo inicialmente la variable binaria, lo cual conseguiremos escribiendo I(class=="RNI") que se interpreta como

\[ \mathbb{I}(class=RNI)=\left\{\begin{array}{cc}1 & \textrm{si}\ \ x=RNI \\ 0 & \textrm{si}\ \ x=UNI\end{array}\right. \]

donde las $x$ son las entradas correspondientes a la columna class en el dataframe Datosfil. Por ende

# Realizamos el ajuste
fitC <- lm(aadt ~ I(class=="RNI"), data=Datosfil)

#  Número de decimales permitidos en las salidas
options(digits=7)

# Veamos la infomación
summary(fitC)

## 
## Call:
## lm(formula = aadt ~ I(class == "RNI"), data = Datosfil)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15955  -2939  -1849   1953  61122 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              17221       1581  10.891  < 2e-16 ***
## I(class == "RNI")TRUE   -13756       2060  -6.679 1.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9875 on 93 degrees of freedom
## Multiple R-squared:  0.3242, Adjusted R-squared:  0.3169 
## F-statistic: 44.61 on 1 and 93 DF,  p-value: 1.732e-09

coef(fitC)

##           (Intercept) I(class == "RNI")TRUE 
##              17221.03             -13756.47

en lo que se puede ver, a primera instancia, que el valor de los betas ha cambiado. Lo cual no es cierto pues al ver la interpretación

\[\begin{align*} &\mathbb{E}[aadt|class=RNI]=\hat{\beta_{0}}+\hat{\beta_{1}}= 17221.03 -13756.47 =3464.56=3464.56=\mu_{1}\\ &\mathbb{E}[aadt|class=RNI]=\hat{\beta_{0}}=17221.03=\mu_{2} \end{align*}\]

concluimos, en realidad, qeu los valores de las esperanzas son los mismos.

Análisis de los supuestos

Homocedasticidad

Para analizar este supuesto utilizaremos las mismas herramientas que en el caso en que $x$ era una variable continua:

# usamos la gráfica que R trae por defecto
plot(fit, 3)

en la cual buscamos que las nubes de puntos de ambas líneas coincidan o sean muy parecidas (al parecer en los puntos de la derecha hay 4 outliers). Pasamos a ver las pruebas de hipótesis

# Prueba basada en los errores studentizados
library(lmtest)
lmtest::bptest(fit)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 5.5846, df = 1, p-value = 0.01812

library(car)
car::ncvTest(fit)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 47.25489, Df = 1, p = 6.233e-12

en ambos casos se rechaza $H_{0}$, es decir, hay evidencia en contra de la homocedasticidad.

Normalidad

# Calculamos los errores
library(broom)
Datosfit=augment(fit)
# Usamos la gráfica que R trae por defecto (Q-Qplot)
plot(fit, 2)

donde vemos indicios de evidencia en contra de la normalidad. Pasamos a ver las pruebas de hipótesis

shapiro.test(Datosfit$.std.resid)

## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfit$.std.resid
## W = 0.7507, p-value = 2.112e-11

library(nortest)
nortest::lillie.test(Datosfit$.std.resid)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfit$.std.resid
## D = 0.1924, p-value = 2.82e-09

library(normtest)
normtest::jb.norm.test(Datosfit$.std.resid)

## 
##  Jarque-Bera test for normality
## 
## data:  Datosfit$.std.resid
## JB = 1013.3, p-value < 2.2e-16

en las que, en todas, se está rechazando $H_{0}$ y por ende también tenemos problemas con el supuesto de la normalidad.

Pruebas por grupos

Podemos complementar el análisis anterior realizando pruebas de homocedasticidad para grupos (donde $H_{0}$ es la hipótesis referente a que las varianzas de los grupos es la misma):

bartlett.test(aadt ~ class, data=Datosfil)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  aadt by class
## Bartlett's K-squared = 67.591, df = 1, p-value < 2.2e-16

rechazándose $H_{0}$. Otra prueba más robusta

library(car)
leveneTest(aadt ~ class, data=Datosfil)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value    Pr(>F)    
## group  1   20.62 1.676e-05 ***
##       93                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

nuevamente encontramos evidencia en contra de la igualdad entre varianzas.

# Prueba no parámetrica
fligner.test(aadt ~ class, data=Datosfil)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  aadt by class
## Fligner-Killeen:med chi-squared = 24.656, df = 1, p-value = 6.852e-07

De tal manera hallamos más evidencia en contra de la homocedasticidad.

Después, para analizar el supuesto de la normalidad por grupos utilizaremos las mismas pruebas, por ejemplo shapiro.test( Datosfil$aadt[Datosfil$class=="RNI"]) donde estamos considerando únicamente a la variable $y$ (aadt) con los datos que son de la clase RNI. Luego consideramos únicamente a la variable $y$ con los datos de la clase UNI y así sucesivamente en cada prueba

#Pruebas para la normalidad por grupos
# RNI:
shapiro.test( Datosfil$aadt[Datosfil$class=="RNI"])

## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfil$aadt[Datosfil$class == "RNI"]
## W = 0.73748, p-value = 1.146e-08

nortest::lillie.test(Datosfil$aadt[Datosfil$class=="RNI"])

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfil$aadt[Datosfil$class == "RNI"]
## D = 0.2533, p-value = 1.141e-09

normtest::jb.norm.test(Datosfil$aadt[Datosfil$class=="RNI"])

## 
##  Jarque-Bera test for normality
## 
## data:  Datosfil$aadt[Datosfil$class == "RNI"]
## JB = 31.822, p-value = 5e-04

#UNI
shapiro.test( Datosfil$aadt[Datosfil$class=="UNI"])

## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfil$aadt[Datosfil$class == "UNI"]
## W = 0.80763, p-value = 1.212e-05

nortest::lillie.test(Datosfil$aadt[Datosfil$class=="UNI"])

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfil$aadt[Datosfil$class == "UNI"]
## D = 0.15971, p-value = 0.01352

normtest::jb.norm.test(Datosfil$aadt[Datosfil$class=="UNI"])

## 
##  Jarque-Bera test for normality
## 
## data:  Datosfil$aadt[Datosfil$class == "UNI"]
## JB = 87.095, p-value < 2.2e-16

donde se rechaza en las primeras tres pruebas para el grupo RNI, más aún, en las tres pruebas para el grupo UNI también vemos p-values pequeños por lo que también se está rechazando $h_{0}$.

Fuerte evidencia contra la normalidad y homocedasticidad. Opción: Usar transformaciones para variable $y$.

Transformación de la variable $y$

Procedemos a utilizar transformaciones del tipo Box-Cox. A favor de lo anterior escribimos

summary(powerTransform(fit))

## bcPower Transformation to Normality 
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1    0.0657           0      -0.0714       0.2028
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                             LRT df    pval
## LR test, lambda = (0) 0.8808064  1 0.34798
## 
## Likelihood ratio test that no transformation is needed
##                            LRT df       pval
## LR test, lambda = (1) 163.5176  1 < 2.22e-16

no rechazándose $H_{0}$ por lo que es necesario realizar una transformación a la $y$. Luego, por la prueba realizada para $\lambda=0$ tenemos que la transformación que debemos realizar es una logarítmica. De tal modo que

Datosfil$Ytransf=bcPower(Datosfil$aadt, lambda=0)

y volvemos a ver los boxplots pero ahora con los datos transformados

boxplot(Ytransf  ~ class, data = Datosfil, col = "white", outline=FALSE)
stripchart(Ytransf  ~ class, data = Datosfil,
           method = "jitter",
           pch = 19,
           col = 2:4,
           vertical = TRUE,
           add = TRUE)

viéndose una nube de puntos similar (no igual, no perfecta pero si similar) en ambos boxplots y un rango intercuartílico parecido. Por otro lado, para el caso de la normalidad vemos ya cierta simetría en el boxplot de RNI y aunque la mediana del segundo boxplot no está a la mitad, vemos que en el resto de los puntos se presenta cierta simetría. Vemos que, gracias a al transformación, ya no es contundente la evidencia en contra de la homocedasticidad ni en contra de la normalidad.

Procedemos a realizar nuevamente el análisis de los spuestos

# Ajuste del modelo con los datos transformados
fitln <- lm(Ytransf  ~ class, data=Datosfil)

# Gráfica para revisar el supuesto de la homocedasticidad
plot(fitln, 3)

con la nube de puntos de ambas clases muy similar. Completamos con las pruebas de hipótesis

lmtest::bptest(fitln)

## 
##  studentized Breusch-Pagan test
## 
## data:  fitln
## BP = 2.4214, df = 1, p-value = 0.1197

car::ncvTest(fitln)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.709613, Df = 1, p = 0.19104

donde no se rechaza $H_{0}$ en ambas pruebas, lo cual es justo lo que buscamos.

Continuamos con el caso de la normalidad

plot(fitln, 2)

que presenta un mejor comportamiento que la Q-Qplot de los datos no transformados pero aún no presenta el comportamiento que quisieramos. De manera análoga, completamos con las pruebas de hipótesis

# Cálculo de los errores
Datosfitln=augment(fitln)

# Prueba 
shapiro.test(Datosfitln$.std.resid)

## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfitln$.std.resid
## W = 0.98477, p-value = 0.3402

nortest::lillie.test(Datosfitln$.std.resid)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfitln$.std.resid
## D = 0.061982, p-value = 0.4929

normtest::jb.norm.test(Datosfitln$.std.resid)

## 
##  Jarque-Bera test for normality
## 
## data:  Datosfitln$.std.resid
## JB = 1.4352, p-value = 0.412

no rechazándose $H_{0}$ en las tres pruebas. Parece que no hay fuerte evidencia en contra de la normalidad. Dado que aún puede presentarse la duda realizamos una prueba con los residuales studentizados y contra una $t$ con $n-3$ grados de libertad

library(MASS)
ResStudent=studres(fitln)
ks.test(ResStudent, "pt", length(ResStudent)-3)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  ResStudent
## D = 0.061617, p-value = 0.8635
## alternative hypothesis: two-sided

con la que tampoco se rechaza $H_{0}$.

Pruebas por grupos

comenzamos con la prueba para homocedasticidad

# H_o: las varianzas de los grupos es la misma vs H_a: al menos un grupo tiene una varianza diferente
bartlett.test(Ytransf  ~ class, data = Datosfil)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ytransf by class
## Bartlett's K-squared = 1.6986, df = 1, p-value = 0.1925

no rechazándose $H_{0}$.

#Otra prueba más robusta 
leveneTest(Ytransf  ~ class, data = Datosfil)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.8586 0.1761
##       93

#No paramétrica 
fligner.test(Ytransf  ~ class, data=Datosfil)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Ytransf by class
## Fligner-Killeen:med chi-squared = 1.5045, df = 1, p-value = 0.22

en ambos casos no se rechaza $H_{0}$. Posteriormente pasamos a analizar las pruebas de normalidad por grupos

#Para normalidad
shapiro.test( Datosfil$Ytransf[Datosfil$class=="RNI"])

## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfil$Ytransf[Datosfil$class == "RNI"]
## W = 0.96632, p-value = 0.1191

nortest::lillie.test(Datosfil$Ytransf[Datosfil$class=="RNI"])

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfil$Ytransf[Datosfil$class == "RNI"]
## D = 0.10262, p-value = 0.15

normtest::jb.norm.test(Datosfil$Ytransf[Datosfil$class=="RNI"])

## 
##  Jarque-Bera test for normality
## 
## data:  Datosfil$Ytransf[Datosfil$class == "RNI"]
## JB = 1.9897, p-value = 0.239

shapiro.test( Datosfil$Ytransf[Datosfil$class=="UNI"])

## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfil$Ytransf[Datosfil$class == "UNI"]
## W = 0.95676, p-value = 0.1388

nortest::lillie.test(Datosfil$Ytransf[Datosfil$class=="UNI"])

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfil$Ytransf[Datosfil$class == "UNI"]
## D = 0.1641, p-value = 0.009707

normtest::jb.norm.test(Datosfil$Ytransf[Datosfil$class=="UNI"])

## 
##  Jarque-Bera test for normality
## 
## data:  Datosfil$Ytransf[Datosfil$class == "UNI"]
## JB = 1.3727, p-value = 0.3435

en cualquier caso no se rechaza, salvo en la prueba Lilliefors. Si somos muy estrictos, al rechazarse $H_{0}$ en esa prueba podemos argumentar la existencia de evidencia en contra d ela normalidad en la clase UNI. Por otro lado, también podemos argumentar que dos pruebas de hipótesis para este grupo y dada la cantidad de datos no hay evidencia concreta en contra de la normalidad.

Observando

levels(Datosfil$class)

## [1] "RNI" "UNI"

R asigna un orden por defecto y al primer nivel (en este caso RNI) se le conoce como nivel de referencia ó valor cero en la variable $z$, la cual R nuca incluye en el ajuste del modelo. Pues por ejemplo, de

fit

## 
## Call:
## lm(formula = aadt ~ class, data = Datosfil)
## 
## Coefficients:
## (Intercept)     classUNI  
##        3465        13756

puede verse como la variable binaria creada automáticamente por R es referente al nivel UNI.

Luego, podemos cambiar el nivel de referencia. Por ejemplo, en vez de que el nivel de referencia sea RNI podemos considerar a UNI. Antes, veamos como

summary(fitln)

## 
## Call:
## lm(formula = Ytransf ~ class, data = Datosfil)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.26656 -0.82325  0.01484  0.68300  2.15391 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.5613     0.1365  55.381  < 2e-16 ***
## classUNI      1.8489     0.2131   8.677 1.29e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.022 on 93 degrees of freedom
## Multiple R-squared:  0.4474, Adjusted R-squared:  0.4414 
## F-statistic: 75.28 on 1 and 93 DF,  p-value: 1.286e-13

en el ajuste del modelo con los datos transformados la variable binaria es nuevamente UNI. Ahora

#Cambio de nivel de referencia y ajuste nuevamente del modelo

fitln3 <- lm(Ytransf  ~ relevel(class, ref = "UNI"), data=Datosfil)
summary(fitln3)

## 
## Call:
## lm(formula = Ytransf ~ relevel(class, ref = "UNI"), data = Datosfil)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.26656 -0.82325  0.01484  0.68300  2.15391 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      9.4102     0.1636  57.518  < 2e-16 ***
## relevel(class, ref = "UNI")RNI  -1.8489     0.2131  -8.677 1.29e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.022 on 93 degrees of freedom
## Multiple R-squared:  0.4474, Adjusted R-squared:  0.4414 
## F-statistic: 75.28 on 1 and 93 DF,  p-value: 1.286e-13

del cual podemos ver como “han cambiado” los valores de las betas y que la variable binaria creada ahora se basa en el nivel RNI.

Tenemos pues que

$\mathbb{E}[log(aadt)|class=RNI]=\hat{\beta_{0}}=\mu_{1}$
$\mathbb{E}[log(aadt)|class=UNI]=\hat{\beta_{0}}+\hat{\beta_{1}}=\mu_{2}$

Supongamos que el investigador tenía la hipótesis de que en las carreteras “urban noninterstate” hay más tráfico. Aunque la $y$ este transformada podemos utilizar los valores de las mu para analizar el supuesto del investigador. Lo planteamos como

\[ H_0: \mu_1\geq\mu_2 \ \ vs \ \ H_a: \mu_1<\mu_2 \]

o equivalentemente

\[ H_0: 0\geq \hat{\beta_{1}} \ \ vs \ \ H_a: 0<\hat{\beta_{1}} \]

Para lo anterior

# Utilizamos el siguiente paquete
library(multcomp)

#theta = zo*bo + z1*b1
MatZ0Z1=matrix(c(0,1), ncol=2, nrow=1)

# El valor con el que comparamos a b1:
c = 0

como en la prueba de hipótesis sólo tenemos involucrada a la $\hat{\beta_{1}}$ entonces $\hat{\theta}=0\cdot \hat{\beta_{0}} +1\cdot \hat{\beta_{1}}$ y por ello definimos en la matriz el vector c(0,1). De tal manera usamos el siguiente código, donde agregaremos alternative="greater" pues en la hipótesis alternativa tenemos que $0<\hat{\beta_{1}}$, así

prueba1=glht(fitln, linfct=MatZ0Z1, rhs=c, alternative ="greater")
summary(prueba1)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Ytransf ~ class, data = Datosfil)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value   Pr(>t)    
## 1 <= 0   1.8489     0.2131   8.677 6.43e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

en la que se rechaza $H_{0}$, es decir, es plausible asumir que $\mu_1<\mu_2$.

Interpretación

En otras palabras hay evidencia con una significancia de .05 de que en promedio el “Average Annual Daily Traffic Data” ($y$) en escala logarítmica es mayor en las carreteras “urban noninterstate”.

En la escala original podemos decir que hay evidencia con una significancia de .05 de que la mediana del “Average Annual Daily Traffic Data” es mayor en las carreteras “urban noninterstate”.

Por otro lado, en realidad a los investigadores no les interesa el resultado en términos de la mediana o de la esperanza, lo que le interesa a ellos es el resultado en términos de su hipótesis, dicho lo anterior podemos interpretar el resultado como

Hay evidencia con una significancia de .05 de que hay mayor probabilidad de observar mayores valores de “Average Annual Daily Traffic Data” (y) en las carreteras “urban noninterstate”

Ejemplo. Comparación de medias cuando se tienen dos grupos con diferente varianza, pero se puede asumir normalidad (en este ejemplo no es así, pero se ejemplifica)

# oneway.test (diferente varianzas, posiblemente más de dos niveles)
# Sólo comparación H_0: mu_1=mu_2 vs H_a: mu_1!=mu_2
oneway.test(aadt ~ class, data = Datosfil, var.equal = FALSE)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  aadt and class
## F = 32.597, num df = 1.000, denom df = 42.035, p-value = 1.037e-06

# usamos var.equal = FALSE asumiendo que las varianzas pueden ser distintas

en la que se rechaza $H_{0}$, es decir, es plausible asumir que hay diferencia entre las medias. Notemos que estamos usando cuando la variable $y$ no está transformada, pues lo que estamos realizando es sólo para ejemplificar el caso en el que se está trabajando. Por otro lado, escribiendo var.equal=TRUE obtendríamos un resultado equivalente a la información de summary(fit)

oneway.test(aadt ~ class, data = Datosfil, var.equal = FALSE)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  aadt and class
## F = 32.597, num df = 1.000, denom df = 42.035, p-value = 1.037e-06

summary(fit)

## 
## Call:
## lm(formula = aadt ~ class, data = Datosfil)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15955  -2939  -1849   1953  61122 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     3465       1320   2.625   0.0101 *  
## classUNI       13756       2060   6.679 1.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9875 on 93 degrees of freedom
## Multiple R-squared:  0.3242, Adjusted R-squared:  0.3169 
## F-statistic: 44.61 on 1 and 93 DF,  p-value: 1.732e-09

Alternativamente a la prueba anterior podemos utilizar, sólo para el caso de dos niveles, el código

#Otra alternativa, sólo dos niveles
levels(Datosfil$class)

## [1] "RNI" "UNI"

#H_0: mu_1=mu_2 vs H_a: mu_1!=mu_2
# Forma 1
t.test(aadt ~ class, data = Datosfil, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  aadt by class
## t = -5.7093, df = 42.035, p-value = 1.037e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -18618.853  -8894.091
## sample estimates:
## mean in group RNI mean in group UNI 
##          3464.554         17221.026

# Forma 2
t.test(Datosfil$aadt[Datosfil$class=="RNI"], Datosfil$aadt[Datosfil$class=="UNI"], var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  Datosfil$aadt[Datosfil$class == "RNI"] and Datosfil$aadt[Datosfil$class == "UNI"]
## t = -5.7093, df = 42.035, p-value = 1.037e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -18618.853  -8894.091
## sample estimates:
## mean of x mean of y 
##  3464.554 17221.026

en los cuales se obtiene un p-value igual al obtenido en la prueba oneway.test. Después, podemos realizar pruebas de la forma

\[ H_0: \mu_1\geq \mu_2 \ \ vs \ \ H_a: \mu_1<\mu_2 \] con la función t.test puesto que ahora la oneway.test ya no nos permite hacer pruebas de este estilo. Luego

# H_0: mu_1>=mu_2 vs H_a: mu_1<mu_2
#Considera como mu_1 a la media asociada al primer nivel de la variable.
# Forma 1
t.test(aadt ~ class, data = Datosfil, alternative ="less", var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  aadt by class
## t = -5.7093, df = 42.035, p-value = 5.183e-07
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf -9703.94
## sample estimates:
## mean in group RNI mean in group UNI 
##          3464.554         17221.026

# Forma 2
t.test(Datosfil$aadt[Datosfil$class=="RNI"], Datosfil$aadt[Datosfil$class=="UNI"], alternative ="less", var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  Datosfil$aadt[Datosfil$class == "RNI"] and Datosfil$aadt[Datosfil$class == "UNI"]
## t = -5.7093, df = 42.035, p-value = 5.183e-07
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf -9703.94
## sample estimates:
## mean of x mean of y 
##  3464.554 17221.026

donde es importante saber a qué clase le corresponde cada valor de $\mu$ para poder interpretar de manera correcta la hipótesis alternativa. Para

levels(Datosfil$class)

## [1] "RNI" "UNI"

a RNI siempre le corresponde el valor de la $\mu_{1}$ y a INI el valor de la $mu_{2}$.

Caso en el que x es una variable categórica

Luis Fernando Apáez Álvarez

17/11/2021

Teoría

Homocedasticidad.

Normalidad e independencia de las observaciones

Ejemplo

Boxplot

Ajuste del modelo

Cálculo “a mano”

Otra forma de ajustar

Análisis de los supuestos

Homocedasticidad

Normalidad

Pruebas por grupos

Transformación de la variable \(y\)

Pruebas por grupos

Interpretación