Cargar las librerías requeridas
library(survival)
library(KMsurv)
library(survMisc)
library(survminer)
## Loading required package: ggplot2
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:survMisc':
##
## autoplot
## Loading required package: ggpubr
## Loading required package: magrittr
library(flexsurv)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#AJUSTAR EL DIRECTORIO DE TRABAJO
setwd("G:/UAAAN/MATERIAS/2019/REGRESION APLICADA/LECCION 7")
#LEER LOS DATOS
mis_datos<-read.csv("Intercepcion_de_lluvia.csv", header = T)
attach(mis_datos)
print(mis_datos)
## X Bueno Malo
## 1 0.7 0.2 0.2
## 2 1.1 0.4 0.4
## 3 11.0 2.5 2.5
## 4 1.6 0.4 0.4
## 5 4.0 0.8 0.8
## 6 7.1 1.6 1.6
## 7 14.4 2.9 2.9
## 8 40.0 8.0 2.0
Resumen de los datos
summary(mis_datos)
## X Bueno Malo
## Min. : 0.700 Min. :0.2 Min. :0.200
## 1st Qu.: 1.475 1st Qu.:0.4 1st Qu.:0.400
## Median : 5.550 Median :1.2 Median :1.200
## Mean : 9.988 Mean :2.1 Mean :1.350
## 3rd Qu.:11.850 3rd Qu.:2.6 3rd Qu.:2.125
## Max. :40.000 Max. :8.0 Max. :2.900
Gráficos de dispersión del bueno y del malo
x <- c(1, 2)
m <- matrix(x, ncol = 2)
layout(m)
nf <- layout(m)
#layout.show(nf)
#HACIENDO LA FIGURA DE DISPERSION DEL BUENO
plot(X, Bueno, xlab="Precipitacion total (cm)", col.lab=3, ylab="Intercepcion (cm)", pch=1, col="green", lwd = 5, col.axis=3)
abline(lm(Bueno~X), col="green", col.axis=1)
#HACIENDO LA FIGURA DE DISPERSION DEL MALO
plot(X, Malo, xlab="Precipitacion total (cm)", col.lab=2, ylab="Intercepcion (cm)", pch=1, col="red", lwd = 5, col.axis=2)
abline(lm(Malo~X), col="red")
Los Gráficos de box plot del bueno y del malo
x0 <- c(1, 2) # dos boxplots juntos
m0 <- matrix(x0, ncol = 2)
layout(m0)
nf0 <- layout(m0)
#layout.show(nf0)
#La dispersion del bueno
boxplot(Bueno, main="El Bueno")
abline(h=median(Bueno)); abline(h=mean(Bueno), col="green")
#La dispersion del malo
boxplot(Malo, main="El Malo")
abline(h=median(Malo)); abline(h=mean(Malo), col="red")
Gráficos qqplot del bueno y del malo
par(new=T)
## Warning in par(new = T): llamada par(new=TRUE) sin gráfico
ggqqplot(Bueno, ylab = "B", conf.int = TRUE, conf.int.level = 0.95, color = "green")
ggqqplot(Malo, ylab = "M", conf.int = TRUE, conf.int.level = 0.95, color = "red")
Prueba de hipótesis de la correlación del bueno y del malo
# Prueba de hipæ¼ã¸³tesis H0:rho=0 (Prueba de hipotesis) Bueno
cor.test(X, Bueno, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: X and Bueno
## t = 61.696, df = 6, p-value = 1.219e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9954645 0.9998636
## sample estimates:
## cor
## 0.9992128
# Prueba de hipæ¼ã¸³tesis H0:rho=0 (Prueba de hipotesis) Malo
cor.test(X, Malo, alternative="two.sided", method="pearson")
##
## Pearson's product-moment correlation
##
## data: X and Malo
## t = 1.8222, df = 6, p-value = 0.1183
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1860527 0.9161947
## sample estimates:
## cor
## 0.5968742
La regresión del bueno
n<-length(X)
plot(X, Bueno, ylab="Intercepcion (mm", xlab="Precipitacion total (mm)", main="", col.lab=3)
# Estimaciæ¼ã¸³n de paræ¼ã¸±metros y prueba de hipæ¼ã¸³tesis
fit_Bueno<-lm(Bueno~X); summary(fit_Bueno)
##
## Call:
## lm(formula = Bueno ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11724 -0.06693 -0.03589 0.05932 0.19999
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.127079 0.050613 2.511 0.0459 *
## X 0.197539 0.003202 61.696 1.22e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.111 on 6 degrees of freedom
## Multiple R-squared: 0.9984, Adjusted R-squared: 0.9982
## F-statistic: 3806 on 1 and 6 DF, p-value: 1.219e-09
abline(fit_Bueno,col="green")
# Podemos obtener el valor predicho de Y para cada valor de X
fitted_bueno<-fit_Bueno$ fitted.value; fitted_bueno
## 1 2 3 4 5 6 7
## 0.2653567 0.3443723 2.3000082 0.4431418 0.9172353 1.5296062 2.9716408
## 8
## 8.0286388
# Con un bucle for y la funciæ¼ã¸³n lines() podemos representar los resæ¼ã¹¤duos
for(i in 1:n)
{
lines( c(X[i],X[i]), c(Bueno[i], fitted_bueno[i]), col="green")
}
La regresión del malo
n<-length(X)
# Grafica
plot(X, Malo, ylab="Intercepcion (mm", xlab="Precipitacion total (mm)", main="", col.lab=2)
# Estimaciæ¼ã¸³n de paræ¼ã¸±metros y prueba de hipæ¼ã¸³tesis
fit_Malo<-lm(Malo~X); summary(fit_Malo)
##
## Call:
## lm(formula = Malo ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7788 -0.5900 -0.3959 0.5660 1.3399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.87454 0.41297 2.118 0.0785 .
## X 0.04761 0.02612 1.822 0.1183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9054 on 6 degrees of freedom
## Multiple R-squared: 0.3563, Adjusted R-squared: 0.249
## F-statistic: 3.321 on 1 and 6 DF, p-value: 0.1183
abline(fit_Malo,col="red")
# Podemos obtener el valor predicho de Y para cada valor de X
fitted_malo<-fit_Malo$ fitted.value; fitted_malo
## 1 2 3 4 5 6 7
## 0.9078684 0.9269104 1.3982001 0.9507129 1.0649649 1.2125405 1.5600572
## 8
## 2.7787457
# Con un bucle for y la funciæ¼ã¸³n lines() podemos representar los resæ¼ã¹¤duos
for(i in 1:n)
{
lines( c(X[i],X[i]), c(Malo[i], fitted_malo[i]), col="red")
}
Análisis de varianza bueno y del malo
# Anæ¼ã¸±lisis de varianza del bueno
confint(fit_Bueno); anova(fit_Bueno)
## 2.5 % 97.5 %
## (Intercept) 0.003233408 0.2509254
## X 0.189704377 0.2053736
## Analysis of Variance Table
##
## Response: Bueno
## Df Sum Sq Mean Sq F value Pr(>F)
## X 1 46.866 46.866 3806.3 1.219e-09 ***
## Residuals 6 0.074 0.012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Anæ¼ã¸±lisis de varianza del Malo
confint(fit_Malo); anova(fit_Malo)
## 2.5 % 97.5 %
## (Intercept) -0.13594743 1.8850371
## X -0.01631962 0.1115297
## Analysis of Variance Table
##
## Response: Malo
## Df Sum Sq Mean Sq F value Pr(>F)
## X 1 2.7218 2.7218 3.3205 0.1183
## Residuals 6 4.9182 0.8197
Graficos de residuales ordinarios del bueno y del malo
x <- c(1, 2) # dos græ¼ã¸±ficos diferentes
m <- matrix(x, ncol = 2)
layout(m)
nf <- layout(m)
#layout.show(nf)
# Residuos ordinarios del Bueno
fit_Bueno$residuals
## 1 2 3 4 5 6
## -0.06535668 0.05562773 0.19999178 -0.04314176 -0.11723533 0.07039382
## 7 8
## -0.07164077 -0.02863878
boxplot(fit_Bueno$residuals)
# Residuos ordinarios del Malo
fit_Malo$residuals
## 1 2 3 4 5 6
## -0.7078684 -0.5269104 1.1017999 -0.5507129 -0.2649649 0.3874595
## 7 8
## 1.3399428 -0.7787457
boxplot(fit_Malo$residuals)
Pruebas de normalidad del bueno y del malo
shapiro.test(fit_Bueno$residuals) #Bueno
##
## Shapiro-Wilk normality test
##
## data: fit_Bueno$residuals
## W = 0.90713, p-value = 0.3343
shapiro.test(fit_Malo$residuals) #Malo
##
## Shapiro-Wilk normality test
##
## data: fit_Malo$residuals
## W = 0.84071, p-value = 0.07661
# Cuantiles del Bueno
qqnorm(fit_Bueno$residuals)
qqline(fit_Bueno$residuals, col="green")
# Cuantiles del Malo
qqnorm(fit_Malo$residuals)
qqline(fit_Malo$residuals, col="red")
Analisis de autocorrelacion de residuales del bueno y del malo
# Prueba de Durbin -Watson del Bueno
library(lmtest, pos=4)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
#Correlograma de resduales del bueno
acf(fit_Bueno$ residuals)
#Correlograma de resduales del malo
acf(fit_Malo$ residuals)
Mas analisis de residuos
# Residuos frente a puntuaciones ajustadas Bueno
oldpar <- par(oma=c(0,0,3,0), mfrow=c(1,2))
plot(fit_Bueno)
par(oldpar)
# Residuos frente a puntuaciones ajustadas Malo
oldpar <- par(oma=c(0,0,3,0), mfrow=c(1,2))
plot(fit_Malo)
par(oldpar)
Los residuos estandarizados del bueno y del malo
# Residuos estandarizados Bueno
ei=fit_Bueno$residuals/summary(fit_Bueno)$sigma
boxplot(ei)
# Residuos estandarizados Malo
ei=fit_Malo$residuals/summary(fit_Malo)$sigma
boxplot(ei)
Los residuos estudentizados
# Residuos estudentizados del Bueno
ri=rstudent(fit_Bueno)
boxplot(ri)
data.frame(ei,ri)
## ei ri
## 1 -0.7818542 -0.6227885
## 2 -0.5819826 0.5224315
## 3 1.2169592 2.8522792
## 4 -0.6082730 -0.3990065
## 5 -0.2926589 -1.1879941
## 6 0.4279565 0.6470784
## 7 1.4799927 -0.6633896
## 8 -0.8601396 -0.6980660
par(new=T)
abline(h=3, col="red");abline(h=-3, col="red")
# Residuos estudentizados del Malo
ri=rstudent(fit_Malo)
boxplot(ri)
data.frame(ei,ri)
## ei ri
## 1 -0.7818542 -0.8522834
## 2 -0.5819826 -0.6123287
## 3 1.2169592 1.4026284
## 4 -0.6082730 -0.6391514
## 5 -0.2926589 -0.2930914
## 6 0.4279565 0.4268821
## 7 1.4799927 1.9227644
## 8 -0.8601396 -18.9818113
par(new=T)
abline(h=3, col="red");abline(h=-3, col="red")
Diagnóstico de datos potenciales e influyentes
# Estadæ¼ã¹¤sticos y graficas de diagnæ¼ã¸³stico del Bueno
inflm.SR=influence.measures(fit_Bueno); summary(inflm.SR)
## Potentially influential observations of
## lm(formula = Bueno ~ X) :
##
## dfb.1_ dfb.X dffit cov.r cook.d hat
## 8 0.54 -1.71_* -1.85_* 9.56_* 1.86_* 0.87_*
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
outlierTest(fit_Bueno)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 3 2.852279 0.035728 0.28583
influencePlot(fit_Bueno)
## StudRes Hat CookD
## 1 -0.6227885 0.1968198 0.05292264
## 3 2.8522792 0.1258536 0.26751009
## 5 -1.1879941 0.1548495 0.12099784
## 8 -0.6980660 0.8749822 1.86458915
# Estadæ¼ã¹¤sticos y graficas de diagnæ¼ã¸³stico del malo
inflm.SR=influence.measures(fit_Malo); summary(inflm.SR)
## Potentially influential observations of
## lm(formula = Malo ~ X) :
##
## dfb.1_ dfb.X dffit cov.r cook.d hat
## 8 14.66_* -46.49_* -50.22_* 0.00 20.71_* 0.87_*
library(car)
outlierTest(fit_Malo)
## rstudent unadjusted p-value Bonferroni p
## 8 -18.98181 7.478e-06 5.9824e-05
influencePlot(fit_Malo)
## StudRes Hat CookD
## 1 -0.8522834 0.1968198 0.09325334
## 7 1.9227644 0.1412112 0.20969383
## 8 -18.9818113 0.8749822 20.70919747
Predicciones del bueno
# Banda de confianza para E(Y/X0) del Bueno
new <- data.frame(X=seq(from=min(X), to=max(X), length=X)); new
## Warning in seq.default(from = min(X), to = max(X), length = X): first
## element used of 'length.out' argument
## X
## 1 0.7
CI95 <- predict( fit_Bueno, newdata=new, se.fit_Bueno=TRUE,interval="confidence", level=0.95)
CI95
## fit lwr upr
## 1 0.2653567 0.1449007 0.3858127
# Græ¼ã¸±fica Bueno
X0<-seq(min(X),max(X),length=length(Bueno))
pred.m<-predict(fit_Bueno,data.frame(X=X0),interval="confidence",se.fit=T)
par(mfrow=c(1,1))
matplot(X0,cbind(pred.m$fit),lty=c(1,2,2,3,3),
col=c("black","red","red","blue","blue"),type="l",xlab="Precipitacion total (cm)",
ylab="Intercepcion (cm)")
points(X,Bueno)
Predicciones del bueno
# Banda de confianza para E(Y/X0) Malo
new <- data.frame(X=seq(from=min(X), to=max(X), length=X)); new
## Warning in seq.default(from = min(X), to = max(X), length = X): first
## element used of 'length.out' argument
## X
## 1 0.7
CI95 <- predict( fit_Malo, newdata=new, se.fit_Malo=TRUE,interval="confidence", level=0.95)
CI95
## fit lwr upr
## 1 0.9078684 -0.07496398 1.890701
# Græ¼ã¸±fica Malo
X0<-seq(min(X),max(X),length=length(Malo))
pred.m<-predict(fit_Malo,data.frame(X=X0),interval="confidence",se.fit=T)
par(mfrow=c(1,1))
matplot(X0,cbind(pred.m$fit),lty=c(1,2,2,3,3),
col=c("black","red","red","blue","blue"),type="l",xlab="Precipitacion total (cm)",
ylab="Intercepcion (cm)")
points(X,Malo)
