Estudiando Agricultura en Taiwán: Regresión Lineal

Las regresiones ayudan a crear ecuaciones para explicar cómo variables independientes afectan o cambian a otras variables dependientes.

# Agregando ecuación de regresión y coeficiente de determinación (R2) a gráfico
library(ggplot2)
library(ggpmisc)

# Región sombreada representa el 95% intervalo de confianza

my.formula <- y ~ x
p <- ggplot(data         = hoja1, aes(x = NFERTI, y = NHOJA))                          +
  geom_smooth(method     = "lm", se=TRUE, color="blue", formula = my.formula)          +
  stat_poly_eq(formula   = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")), 
               parse     = TRUE) + geom_point(); p

El Capítulo 15 de Experimentación en Agricultura, inicia la descripción de modelos estadísticos lineales. Colectamos códigos para los ejemplos que ofrece el Capítulo.

Códigos en R

#***********************************************************
# EXPERIMENTACIÓN EN AGRICULTURA                              
# Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010       
# CAPÍTULO 15: REGRESIÓN LINEAL 
#************************************************************

# Importar Tabla 15.1 & 15.2, Hoja de cálculo 1
library(readxl)

hoja1 <- read_excel("C:/Users/Administrator/Desktop/Tabla 15.1 & 15.2.xls", sheet = 1)
attach(hoja1); knitr::kable(hoja1)

NFERTI	NHOJA
0.00	1.41
0.12	1.62
0.25	1.74
0.50	1.78
1.00	1.93

# Matrix de dispersión y correlación
library(GGally)

ggpairs(hoja1)

# Diagramas de dispersión básicos con líneas de regresión
library(ggplot2)

ggplot(hoja1, aes(x=NFERTI, y=NHOJA)) + 
  geom_point(shape=1)       # Con círculos vacíos

ggplot(hoja1, aes(x=NFERTI, y=NHOJA)) + 
  geom_point(shape=1) +     # Con círculos vacíos
  geom_smooth(method=lm,    # Agregar línea de regresión
              se=FALSE)     # DSuprimir región sombreada de intervalo de confianza

ggplot(hoja1, aes(x=NFERTI, y=NHOJA)) + 
  geom_point(shape=1) +     # Con círculos vacíos
  geom_smooth(method=lm)    # Agregar línea de regresión

                            # Agrega de forma predeterminada la región de 95% de intervalo de confianza 


# Agregando ecuación de regresión y coeficiente de determinación (R2) a gráfico
library(ggplot2)
library(ggpmisc)

my.formula <- y ~ x
p <- ggplot(data         = hoja1, aes(x = NFERTI, y = NHOJA))                          +
  geom_smooth(method     = "lm", se=FALSE, color="blue", formula = my.formula)         +
  stat_poly_eq(formula   = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")), 
               parse     = TRUE) + geom_point(); p

# Región sombreada representa el 95% intervalo de confianza
my.formula <- y ~ x
p <- ggplot(data         = hoja1, aes(x = NFERTI, y = NHOJA))                          +
  geom_smooth(method     = "lm", se=TRUE, color="blue", formula = my.formula)          +
  stat_poly_eq(formula   = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")), 
               parse     = TRUE) + geom_point(); p

# Análisis de varianza para modelo de regresión (página 204)
reg1 <- lm(NHOJA ~ NFERTI, hoja1)
summary(reg1)

## Call:
## lm(formula = NHOJA ~ NFERTI, data = hoja1)
## 
## Residuals:
##        1        2        3        4        5 
## -0.12038  0.03648  0.09891  0.02820 -0.04321 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.53038    0.06313  24.243 0.000154 ***
## NFERTI       0.44282    0.12254   3.614 0.036408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09707 on 3 degrees of freedom
## Multiple R-squared:  0.8132, Adjusted R-squared:  0.7509 
## F-statistic: 13.06 on 1 and 3 DF,  p-value: 0.03641

# Intervalos de Confianza (95%) e Intervalos de Predicción
var_temp <- predict(reg1, interval="prediction") # # Aparecerá un 'Warning'-Advertencia

nuevo    <- cbind(hoja1, var_temp)

ggplot(nuevo, aes(NFERTI, NHOJA))+ geom_point()               +
  geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  geom_line(aes(y = upr), color = "red", linetype = "dashed") +
  geom_smooth(method =lm, se    = TRUE)

# Valores ajustados, intervalos, residuos y Cook´s D
library(MASS) # Para obtener Normalización y Estudentization: stdres(), studres()

fit       <- lm(NHOJA ~ NFERTI)
ajustado  <- predict(fit, interval = "prediction")

tablaDiag <- data.frame(hoja1, ajustado, reg1$resid, stdres(reg1), studres(reg1), cooks.distance(reg1))
knitr::kable(tablaDiag)

NFERTI	NHOJA	fit	lwr	upr	reg1.resid	stdres.reg1.	studres.reg1.	cooks.distance.reg1.
0.00	1.41	1.530384	1.161883	1.898886	-0.1203844	-1.6325072	-3.9893242	0.9764992
0.12	1.62	1.583523	1.230915	1.936131	0.0364769	0.4500414	0.3805268	0.0439841
0.25	1.74	1.641090	1.299243	1.982937	0.0989100	1.1570697	1.2695941	0.1937897
0.50	1.78	1.751796	1.409838	2.093754	0.0282044	0.3301101	0.2745666	0.0158458
1.00	1.93	1.973207	1.555933	2.390481	-0.0432069	-1.0624388	-1.0983886	2.6511994

# Gráficos de diagnósticos de residuos 
par(mfrow=c(2,2)) # Cambiar distribución de panel a 2 x 2
plot(reg1)

# Ajuste de recta de regresión por el origen
regZero <- lm(NHOJA ~ NFERTI - 1, hoja1)
summary(regZero)

## Call:
## lm(formula = NHOJA ~ NFERTI - 1, data = hoja1)
## 
## Residuals:
##       1       2       3       4       5 
##  1.4100  1.3080  1.0901  0.4802 -0.6696 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)  
## NFERTI    2.600      1.024   2.538   0.0641 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.18 on 4 degrees of freedom
## Multiple R-squared:  0.617,  Adjusted R-squared:  0.5212 
## F-statistic: 6.444 on 1 and 4 DF,  p-value: 0.06409

# Gráfico con recta por el origen
ggplot(hoja1, aes(x  = NFERTI, y = NHOJA)) + geom_point(shape = 1) + 
  geom_smooth(method = "lm",  se = FALSE, formula = y ~ x - 1)

# Gráficos de diagnósticos de regresión por el origen
plot(regZero)

# Importar Tabla 15.1 & 15.2, Hoja de cálculo 2
library(readxl)
hoja2 <- read_excel("C:/Users/Administrator/Desktop/Tabla 15.1 & 15.2.xls", sheet = 2)
attach(hoja2)

# Cambiar la variable 'MUESTRA' a factor y etiquetar valores
MUESTRA <- factor(MUESTRA, levels = c(1, 2), labels = c("Año 1", "Año 2"))
hoja2   <- data.frame(MUESTRA, NFERTI, NHOJA)
knitr::kable(hoja2)

MUESTRA	NFERTI	NHOJA
Año 1	0.00	1.41
Año 1	0.12	1.62
Año 1	0.25	1.74
Año 1	0.50	1.78
Año 1	1.00	1.93
Año 2	0.00	1.45
Año 2	0.12	1.60
Año 2	0.25	1.56
Año 2	0.50	1.70
Año 2	1.00	1.75

# Gráfico de ambas muestras
ggplot(hoja2, aes(x = NFERTI, y = NHOJA, color = MUESTRA)) + geom_point(shape = 1) # MUESTRA fijará color

# Línea de regresión para cada año
ggplot(hoja2, aes(x     = NFERTI, y = NHOJA, color = MUESTRA)) +
  geom_point(shape      = 1) + scale_colour_hue(l = 50)        +  # Usar una paleta relativamente oscura que la normal
  geom_smooth(method    = lm,                                     # Agrega líneas de regresión
              se        = FALSE,                                  # No agrega la región de confianza 
              fullrange = TRUE)                                   # Extender líneas de regresión

# Regressión simple (página 213)
reg2 <- lm(NHOJA ~ NFERTI, hoja2)
summary(reg2)

## Call:
## lm(formula = NHOJA ~ NFERTI, data = hoja2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12618 -0.06594  0.01871  0.05557  0.13001 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.52126    0.04200  36.219  3.7e-10 ***
## NFERTI       0.35492    0.08153   4.353  0.00244 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09134 on 8 degrees of freedom
## Multiple R-squared:  0.7031, Adjusted R-squared:  0.666 
## F-statistic: 18.95 on 1 and 8 DF,  p-value: 0.002435

# Intervalos de Confianza (95%) e Intervalos de Predicción
var_temp <- predict(reg2, interval = "prediction") # # Aparecerá un 'Warning'-Advertencia

nuevo    <- cbind(hoja2,  var_temp)

ggplot(nuevo, aes(NFERTI, NHOJA)) + geom_point()                 +
  geom_line(aes(y    = lwr), color = "red", linetype = "dashed") +
  geom_line(aes(y    = upr), color = "red", linetype = "dashed") +
  geom_smooth(method = lm, se = TRUE)

# Valores ajustados, intervalos, residuos y Cook´s D
library(MASS) # Para obtener Normalización y Estudentization: stdres(), studres()
fit        <- lm(NHOJA ~ NFERTI)
ajustado   <- predict(fit, interval = "prediction")

tablaDiag2 <- data.frame(hoja2, ajustado, reg2$resid, stdres(reg2), studres(reg2), cooks.distance(reg2))
knitr::kable(tablaDiag2)

MUESTRA	NFERTI	NHOJA	fit	lwr	upr	reg2.resid	stdres.reg2.	studres.reg2.	cooks.distance.reg2.
Año 1	0.00	1.41	1.521260	1.289427	1.753092	-0.1112596	-1.3717089	-1.4672075	0.2522761
Año 1	0.12	1.62	1.563850	1.337836	1.789864	0.0561499	0.6673258	0.6423608	0.0397271
Año 1	0.25	1.74	1.609990	1.387852	1.832128	0.1300102	1.5106758	1.6714893	0.1442827
Año 1	0.50	1.78	1.698720	1.476542	1.920898	0.0812800	0.9446584	0.9374793	0.0566441
Año 1	1.00	1.93	1.876181	1.625872	2.126489	0.0538195	0.7685622	0.7470318	0.2071483
Año 2	0.00	1.45	1.521260	1.289427	1.753092	-0.0712596	-0.8785526	-0.8645764	0.1034875
Año 2	0.12	1.60	1.563850	1.337836	1.789864	0.0361499	0.4296315	0.4066016	0.0164666
Año 2	0.25	1.56	1.609990	1.387852	1.832128	-0.0499898	-0.5808651	-0.5551832	0.0213315
Año 2	0.50	1.70	1.698720	1.476542	1.920898	0.0012800	0.0148760	0.0139155	0.0000140
Año 2	1.00	1.75	1.876181	1.625872	2.126489	-0.1261805	-1.8019048	-2.1867069	1.1386417

# Comparación de líneas de regresión por muestras
# Opción 1
modelo1 <- aov(NHOJA ~ NFERTI * MUESTRA, data = hoja2)
modelo2 <- aov(NHOJA ~ NFERTI + MUESTRA, data = hoja2)
anova(modelo1, modelo2)

## Analysis of Variance Table
## 
## Model 1: NHOJA ~ NFERTI * MUESTRA
## Model 2: NHOJA ~ NFERTI + MUESTRA
##   Res.Df      RSS Df  Sum of Sq      F Pr(>F)
## 1      6 0.039407                            
## 2      7 0.049104 -1 -0.0096973 1.4765   0.27

MUESTRA1 <- subset(hoja2, MUESTRA == "Año 1")
MUESTRA2 <- hoja2[hoja2$MUESTRA   == "Año 2",]

regMUESTRA1 <- lm(NHOJA ~ NFERTI, data = MUESTRA1); summary(regMUESTRA1)

## Call:
## lm(formula = NHOJA ~ NFERTI, data = MUESTRA1)
## 
## Residuals:
##        1        2        3        4        5 
## -0.12038  0.03648  0.09891  0.02820 -0.04321 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.53038    0.06313  24.243 0.000154 ***
## NFERTI       0.44282    0.12254   3.614 0.036408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09707 on 3 degrees of freedom
## Multiple R-squared:  0.8132, Adjusted R-squared:  0.7509 
## F-statistic: 13.06 on 1 and 3 DF,  p-value: 0.03641

regMUESTRA2 <- lm(NHOJA ~ NFERTI, data = MUESTRA2); summary(regMUESTRA2)

## Call:
## lm(formula = NHOJA ~ NFERTI, data = MUESTRA2)
## 
## Residuals:
##        6        7        8        9       10 
## -0.06213  0.05582 -0.01889  0.05436 -0.02915 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.51213    0.03962  38.161 3.96e-05 ***
## NFERTI       0.26702    0.07692   3.471   0.0403 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06093 on 3 degrees of freedom
## Multiple R-squared:  0.8007, Adjusted R-squared:  0.7342 
## F-statistic: 12.05 on 1 and 3 DF,  p-value: 0.0403

# Opción 2
regMUESTRAS <- lm(NHOJA ~ MUESTRA / NFERTI - 1, data = hoja2)
summary(regMUESTRAS)

## Call:
## lm(formula = NHOJA ~ MUESTRA/NFERTI - 1, data = hoja2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.120384 -0.039694  0.004657  0.049886  0.098910 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## MUESTRAAño 1          1.5304     0.0527  29.038 1.11e-07 ***
## MUESTRAAño 2          1.5121     0.0527  28.692 1.19e-07 ***
## MUESTRAAño 1:NFERTI   0.4428     0.1023   4.328  0.00494 ** 
## MUESTRAAño 2:NFERTI   0.2670     0.1023   2.610  0.04012 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08104 on 6 degrees of freedom
## Multiple R-squared:  0.9986, Adjusted R-squared:  0.9976 
## F-statistic:  1048 on 4 and 6 DF,  p-value: 1.165e-08

Códigos en SAS

/***********************************************************

EXPERIMENTACIÓN EN AGRICULTURA

Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010

CAPÍTULO 15: REGRESIÓN LINEAL

************************************************************/

ODS HTML;

ODS GRAPHICS ON;

* Importar 'Hoja 1' de 'Tabla 15.1 & 15.2.xls';

PROC IMPORT

OUT = hoja1

DATAFILE = "C:\Users\Administrator\Desktop\Tabla 15.1 & 15.2.xls"

DBMS = XLS REPLACE;

SHEET = "Sheet1";

GETNAMES = YES;

DATAROW = 2;

RUN;

* Mostrar datos de 'Hoja de cálculo 1';

PROC PRINT DATA = hoja1;

RUN;

* Matrix de dispersión;

PROC INSIGHT DATA = hoja1;

SCATTER NHOJA NFERTI * NHOJA NFERTI;

RUN;

* Regressión simple (página 204);

PROC REG DATA = hoja1;

MODEL NHOJA = NFERTI / P R CLM CLI; *datos predecidos, residuos, 95% CI;

PLOT NHOJA * NFERTI;

PLOT NHOJA * NFERTI / PRED;

PLOT RESIDUAL. * PREDICTED.;

MODEL NHOJA = NFERTI / NOINT; * regresión sin intercepto (página 210);

RUN;

* Importar 'Hoja 2' de 'Tabla 15.1 & 15.2.xls';

PROC IMPORT

OUT = hoja2

DATAFILE = "C:\Users\Administrator\Desktop\Tabla 15.1 & 15.2.xls"

DBMS = XLS REPLACE;

SHEET = "Sheet2";

GETNAMES = YES;

DATAROW = 2;

RUN;

* Mostrar datos de 'Hoja de cálculo 2';

PROC PRINT DATA = hoja2;

RUN;

* Regressión simple (página 213);

PROC REG DATA = hoja2;

MODEL NHOJA = NFERTI / P R CLM CLI; *datos predecidos, residuos, 95% CI;

PLOT NHOJA * NFERTI;

PLOT NHOJA * NFERTI / PRED;

PLOT RESIDUAL. * PREDICTED.;

RUN;

* Comparación de líneas de regresión por muestras;

PROC REG DATA = hoja2;

BY MUESTRA;

MODEL NHOJA = NFERTI;

slope: test NHOJA = 1;

RUN;

* Comparación de pendientes. Ver tercera línea en el output de 'Type 3 Tests of Fixed Effects';

PROC MIXED DATA = hoja2;

CLASS MUESTRA;

MODEL NHOJA = NFERTI MUESTRA NFERTI*MUESTRA;

RUN;

QUIT;

ODS GRAPHICS OFF;

ODS HTML CLOSE;

-----------------------------------------------------------------------

En Excel, seleccionando ‘Data’ > ‘Data Analysis’ > ‘Regression’, se logra realizar los mismos análisis. Las regresiones se ejecutan con los datos ordenados de la misma manera. El único cambio necesario, para graficar dos líneas, es el de readecuar la variable MUESTRA (Años) como se observa en la siguiente imagen.

Archivos adjuntos:
Tabla 15.1 & 15.2.xls, capítulo 15.r, capítulo 15.sas

Estudiando Agricultura en Taiwán

Páginas

mayo 16, 2017

Regresión Lineal

No hay comentarios:

Publicar un comentario