# Del libro R-By-Example by Jim Albert and Maria Rizzo
# Capitulo 7 Regresion
setwd("C:\\Users\\Luis\\Documents\\R-by-Example")
plot(cars)
names(cars)
## [1] "speed" "dist"
# Para ajustar un modelo de linea recta
# dist = ß0+ß1 speed+e,
# necesitamos estimados del y-intercept ß0 y
# la pendiente ß1 de la linea.
# La funcion lm
model <- lm(cars$dist~cars$speed)
model
##
## Call:
## lm(formula = cars$dist ~ cars$speed)
##
## Coefficients:
## (Intercept) cars$speed
## -17.579 3.932
coef(model)
## (Intercept) cars$speed
## -17.579095 3.932409
summary(model)
##
## Call:
## lm(formula = cars$dist ~ cars$speed)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.069 -9.525 -2.272 9.215 43.201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.5791 6.7584 -2.601 0.0123 *
## cars$speed 3.9324 0.4155 9.464 1.49e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438
## F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12
plot(residuals(model))
abline(h=-2)
abline(model)
plot(cars,main="dist=2.909 x speed",
xlim=c(0,25))
# linea con intercepto en cero,pendiente=2.909
abline(0,2.909)
plot(model,which = 1,add.smooth = FALSE)
plot(model,which = 1,add.smooth = TRUE)
# un modelo cuadratico podria ser considerado
# para esta data
# 7.3 Analisis de regresion para la data
# con dos variables independientes
setwd("C:\\Users\\Luis\\Documents\\R-by-Example\\R-by-Example\\Rx-data")
arboles <- read.table("C:\\Users\\Luis\\Documents\\R-by-Example\\R-by-Example\\Rx-data\\cherry.txt",
header=TRUE)
names(arboles)
## [1] "Diam" "Height" "Volume"
pairs(arboles)
cor(arboles)
## Diam Height Volume
## Diam 1.0000000 0.5192801 0.9671194
## Height 0.5192801 1.0000000 0.5982497
## Volume 0.9671194 0.5982497 1.0000000
# La correlacion entre diametro y volumen
# es 0.97, indicando fuerte y positiva
# relacion lineal entre
# Diam and Volume. La correlacion entre
# height and volume is 0.60, indica
# una asociacion moderada entre Height and Volume.
# primero, usaremos un modelo de regresion
# lineal con diameter como variable independiente
# Y = ß0+ß1X1+e,
modelo.arbol <- lm(formula = arboles$Volume~arboles$Diam,
data = arboles)
plot(arboles$Diam,arboles$Volume)
abline(modelo.arbol$coefficients)
modelo.arbol
##
## Call:
## lm(formula = arboles$Volume ~ arboles$Diam, data = arboles)
##
## Coefficients:
## (Intercept) arboles$Diam
## -36.943 5.066
summary(modelo.arbol)
##
## Call:
## lm(formula = arboles$Volume ~ arboles$Diam, data = arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.065 -3.107 0.152 3.495 9.587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.9435 3.3651 -10.98 7.62e-12 ***
## arboles$Diam 5.0659 0.2474 20.48 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.9331
## F-statistic: 419.4 on 1 and 29 DF, p-value: < 2.2e-16
plot(modelo.arbol)
# Para predecir el volumen de arboles nuevos,
# podemos usar el metodo predict. Guardamos
# # el valor del diametro para el nuevo arbol(es)
# en una dataframe llamada Diam del
# # la formula del modelo original especificado
# en la llamada de funcion lm . Por ejemplo, la prediccion
# # del volumen para un nuevo arbol con un diametero de 16
# pulgadas es obtenido:
new <- data.frame(Diam=16)# da error
predict(modelo.arbol,new)# da error
## Warning: 'newdata' had 1 row but variables found have 31 rows
## 1 2 3 4 5 6 7
## 5.103149 6.622906 7.636077 16.248033 17.261205 17.767790 18.780962
## 8 9 10 11 12 13 14
## 18.780962 19.287547 19.794133 20.300718 20.807304 20.807304 22.327061
## 15 16 17 18 19 20 21
## 23.846818 28.406089 28.406089 30.432431 32.458774 32.965360 33.978531
## 22 23 24 25 26 27 28
## 34.991702 36.511459 44.110244 45.630001 50.695857 51.709028 53.735371
## 29 30 31
## 54.241956 54.241956 67.413183
# 7.3.2 Modelo de regresion multiple
# Y = ß0+ß1X1+ß2X2+e,
# dos variables independientes
modelo.arbol1 <- lm(Volume~Diam+Height,data=arboles)
modelo.arbol1
##
## Call:
## lm(formula = Volume ~ Diam + Height, data = arboles)
##
## Coefficients:
## (Intercept) Diam Height
## -57.9877 4.7082 0.3393
summary(modelo.arbol1)
##
## Call:
## lm(formula = Volume ~ Diam + Height, data = arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4065 -2.6493 -0.2876 2.2003 8.4847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
## Diam 4.7082 0.2643 17.816 < 2e-16 ***
## Height 0.3393 0.1302 2.607 0.0145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
## F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16
plot(modelo.arbol1)
# modelo.arbol1 tambien muestra problemas
# cuando revisamos los graficos de los residuales
# Finalmente, usaremos un modelo
# que incluye tambien el cuadrado del diameter
# como variable independiente
modelo.arbol2 <- lm(arboles$Volume~arboles$Diam+
I(arboles$Diam^2)+arboles$Height,data = arboles)
modelo.arbol2
##
## Call:
## lm(formula = arboles$Volume ~ arboles$Diam + I(arboles$Diam^2) +
## arboles$Height, data = arboles)
##
## Coefficients:
## (Intercept) arboles$Diam I(arboles$Diam^2)
## -9.9204 -2.8851 0.2686
## arboles$Height
## 0.3764
plot(modelo.arbol2)
# El Modelo (modelo.arbol2), el grafico de residuals vs fits
# no tiene la forma de “U”
# aparente en los modelos anteriores
# Los residuales estan aproximadamente
# centrados en 0 con varianza constante.