Ejercicio_RLS_Ratones
Hugo Salazar Cedillo
13/11/2020
Verificacción de supuestos
Veremos si nuestro modelo cumple normalidad y misma varianza, no es necesario verificar independencia ya que el problema lo supone y haremos transformaciones de ser necesario.
Comenzaremos por importar los datos y las librerias que vamos a usar.
library(readr)
library(leaps)
library(SuppDists)
library(car)
library(carData)
library(broom)
library(lmtest)
library(nortest)
library(normtest)
Datos <- read_table2("C:/Users/chiva/Downloads/Iron.txt")
A <- Datos$Fe3
B <- Datos$Fe2 Ajustamos el modelo de regresión lineal simple
fit1 <- lm(B~A, data=Datos)
summary(fit1)##
## Call:
## lm(formula = B ~ A, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7407 -0.2778 -0.0545 0.7512 4.6269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.05019 0.92823 0.054 0.958
## A 1.16799 0.09507 12.286 1.46e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.135 on 16 degrees of freedom
## Multiple R-squared: 0.9042, Adjusted R-squared: 0.8982
## F-statistic: 150.9 on 1 and 16 DF, p-value: 1.462e-09Verificamos
Homocedasticidad
Una forma para verificar es graficar \(B\) v.s. \(A\)
plot(A, B)
Si no nos queda claro tambien podemos graficar \(e_s\) v.s.\(\hat{y}\) o \(e_s\) v.s. \(x\)
Datosfit <- augment(fit1)
plot(A, sqrt(abs(Datosfit$.std.resid)), xlab='x', ylab='e_s')
plot(fit1,3)
Se pude ver un patrón por lo que tal vez no cumple, veremos algunas pruebas que nos ayudarán a verificar el supuesto.
Primero tenemos que proponer las hipotesis:
\(H_0: p-value > 0.05\) v.s. \(H_a: p-value \leq 0.05\)
Y se busca no rechazar \(H_0\)
lmtest::bptest(fit1)##
## studentized Breusch-Pagan test
##
## data: fit1
## BP = 8.5786, df = 1, p-value = 0.003401car::ncvTest(fit1)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 14.45858, Df = 1, p = 0.00014328De las dos pruebas se rechaza \(H_0\) por lo que la varianza no es constante.
Normalidad
Como \(n\) es pequeña entonces tenemos que hacer pruebas y gráficas
plot(fit1,2)
qqPlot(Datosfit$.std.resid, dist="norm")
## [1] 18 15De las gráficas parece ser que no hay normalidad, procederemos con las prubas contrastanto las hipotesis:
\(H_0: p-value > 0.05\) v.s. \(H_a: p-value \leq 0.05\)
Buscamos no rechazar \(H_0\)
shapiro.test(Datosfit$.std.resid)##
## Shapiro-Wilk normality test
##
## data: Datosfit$.std.resid
## W = 0.88303, p-value = 0.02937nortest::lillie.test(Datosfit$.std.resid)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Datosfit$.std.resid
## D = 0.21886, p-value = 0.02259De las pruebas se rachaza \(H_0\) por lo que no cumple normalidad y se procede a hacer una transformación
Para la transformación usamos la transformación Box-Cox para la variable \(B\)
summary(powerTransform(fit1))## bcPower Transformation to Normality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1 -0.0407 0 -0.329 0.2476
##
## Likelihood ratio test that transformation parameter is equal to 0
## (log transformation)
## LRT df pval
## LR test, lambda = (0) 0.07761778 1 0.78055
##
## Likelihood ratio test that no transformation is needed
## LRT df pval
## LR test, lambda = (1) 31.86793 1 1.6502e-08Se rechaza la hipótesis de que no hay que transformar los datos, pero no se rechaza la hipótesis de la log transformación, además nuestra \(\lambda =0\) por lo que transformamos \(B\) a \(ln(B)\)
Datos$lnB <- log(B)Graficamos la transformación v.s. \(A\)
plot(A, Datos$lnB)
lmtest::bptest(fit2)##
## studentized Breusch-Pagan test
##
## data: fit2
## BP = 2.3584, df = 1, p-value = 0.1246car::ncvTest(fit2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.55471, Df = 1, p = 0.10997shapiro.test(Datosfit2$.std.resid)##
## Shapiro-Wilk normality test
##
## data: Datosfit2$.std.resid
## W = 0.93601, p-value = 0.2473nortest::lillie.test(Datosfit2$.std.resid)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Datosfit2$.std.resid
## D = 0.19756, p-value = 0.06151Vemos que con la transformación ya se ha corregido la homocedasticidad y la normalidad, pero vamos a intentar las transformaciones en \(A\).
Usamos la prueba Box-Tidwell para ver si necesitamos transformaciones tipo potencia en la variable \(A\)
boxTidwell(lnB~A, data = Datos)## MLE of lambda Score Statistic (z) Pr(>|z|)
## 0.77849 -0.9647 0.3347
##
## iterations = 3No se rechaza \(H_0:p-value > 0.05\) por lo que no hay evidencia estadística para considerar hacer una transformación tipo potencia sobre la variable \(A\).
Si hacemos una transformacion exponencial sobre \(A\)
Datos$eA=exp(A)
plot(Datos$eA,Datos$lnB)
Y podemos ver que no hay mejoría.
Si hacemos una transformacion logartimica
Datos$lnA <- log(A)
par(mfrow=c(1,2))
plot(A, Datos$lnB)
plot(Datos$lnA, Datos$lnB)
No podemos apreciar una mejoria
fitp <- lm(lnB~lnA, data=Datos)
par(mfrow=c(1,2))
plot(fitp,1)
plot(fit2,1)
Vemos que no transformar \(A\) es más factible para la modelación
Entonces nuestro modelo de regresión lineal queda: \[\hat{B}^{*}=\hat{\beta_0}+\hat{\beta_1}A\]
donde
\[\hat{B}^{*}=\ln(B)\]
Prueba de Hipotesis sobre la igualdad de medias en los dos grupos
Se busca ver que la varianza entre grupos sea la misma, contrastando:
\(H_0: \sigma_A^{2}=\sigma_B^{2}\) v.s. \(H_a:\sigma_A^{2} \neq \sigma_B^{2}\)
donde buscamos no rechazar \(H_0\)
para ello usamos la prueba de Bartlett y la de Levene
Datos1 = data.frame(cam = c(A,Datos$lnB),
Metodo = as.factor(c(rep("A",18),
rep("B",18))))
bartlett.test(cam ~ Metodo, data=Datos1)##
## Bartlett test of homogeneity of variances
##
## data: cam by Metodo
## Bartlett's K-squared = 51.814, df = 1, p-value = 6.103e-13car::leveneTest(cam ~ Metodo, data=Datos1)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 13.629 0.0007757 ***
## 34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se rechaza \(H_0\) por lo que no se supone varianzas iguales
para la alatoriedad vemos test para cada grupo
plot(fit2, 3)
lmtest::bptest(fit2)##
## studentized Breusch-Pagan test
##
## data: fit2
## BP = 2.3584, df = 1, p-value = 0.1246car::ncvTest(fit2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.55471, Df = 1, p = 0.10997la cumple, Veamos ahora la normalidad
Datosfit2 = broom::augment(fit2)
head(Datosfit2)## # A tibble: 6 x 8
## lnB A .fitted .resid .std.resid .hat .sigma .cooksd
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.40 2.2 1.47 -0.0724 -0.710 0.127 0.111 0.0366
## 2 1.43 2.93 1.54 -0.119 -1.15 0.111 0.108 0.0827
## 3 1.49 3.08 1.56 -0.0736 -0.713 0.108 0.111 0.0307
## 4 1.60 3.49 1.60 -0.00684 -0.0660 0.0996 0.113 0.000241
## 5 1.70 4.11 1.67 0.0366 0.351 0.0888 0.112 0.00600
## 6 1.75 4.95 1.75 -0.000608 -0.00579 0.0765 0.113 0.00000139plot(fit2, 2)
shapiro.test(Datosfit2$.std.resid)##
## Shapiro-Wilk normality test
##
## data: Datosfit2$.std.resid
## W = 0.93601, p-value = 0.2473nortest::lillie.test(Datosfit2$.std.resid)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Datosfit2$.std.resid
## D = 0.19756, p-value = 0.06151normtest::jb.norm.test(Datosfit2$.std.resid)##
## Jarque-Bera test for normality
##
## data: Datosfit2$.std.resid
## JB = 2.2813, p-value = 0.1025Con esto vemos que hay normalidad
Como no hay homocedasticidad pero si normalidad usamos \(t.test\)
fitfinal=t.test(A,Datos$lnB, var.equal = F )
fitfinal##
## Welch Two Sample t-test
##
## data: A and Datos$lnB
## t = 4.7355, df = 17.377, p-value = 0.0001807
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 3.394439 8.833250
## sample estimates:
## mean of x mean of y
## 8.203889 2.090045fit7 = oneway.test(cam ~ Metodo, data=Datos1,var.equal = F)
fit7##
## One-way analysis of means (not assuming equal variances)
##
## data: cam and Metodo
## F = 22.425, num df = 1.000, denom df = 17.377, p-value = 0.0001807Se recahaza \(H_0:\mu_A=\mu_B\) pues hay evidencia de que las medias difieren.