Estudiando Agricultura en Taiwán: Regresión con 2 variables independientes

Las fórmulas para los cálculos manuales de esta regresión se encuentran entre la página 215 y la 226.

library(psych); pairs.panels(NECTARINA)

En R, si desean obtener coeficientes de regresión sin demasiados códigos, summary(lm(PRODU ~ RAMOS + VIGOR, data = NECTARINA) bastará. Por el contrario, si quieren asegurarse que seleccionarán el mejor modelo, pueden revisar la siguiente colección de comandos.

Códigos para el software R

#***********************************************************
# EXPERIMENTACIÓN EN AGRICULTURA                              
# Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010       
# CAPÍTULO 16: REGRESIÓN MÚLTIPLE
# REGRESIÓN CON DOS VARIABLES INDEPENDIENTES
#************************************************************

#### Tabla 16.1. Producción, número de ramos fructíferos y vigor en once ramas principales de la nectarina 'Armking'
# Página 216
NECTARINA <- data.frame(
  RAMA  = c( 1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 10.0, 11.0),
  PRODU = c(18.5, 16.5, 22.0, 18.6, 15.1, 10.8, 32.4, 13.6, 17.4, 16.7, 14.9),
  RAMOS = c(34.0, 36.0, 43.0, 46.0, 35.0, 33.0, 52.0, 32.0, 29.0, 36.0, 39.0),
  VIGOR = c(26.0, 28.0, 25.0, 30.0, 25.0, 22.0, 32.0, 24.0, 24.0, 23.0, 26.0))
attach(NECTARINA); knitr::kable(NECTARINA)

RAMA	PRODU	RAMOS	VIGOR
1	18.5	34	26
2	16.5	36	28
3	22.0	43	25
4	18.6	46	30
5	15.1	35	25
6	10.8	33	22
7	32.4	52	32
8	13.6	32	24
9	17.4	29	24
10	16.7	36	23
11	14.9	39	26

#### Subconjunto, eliminando RAMA de datos iniciales
NECTARINA <- NECTARINA[, c("PRODU", "RAMOS", "VIGOR")]; knitr::kable(NECTARINA)

PRODU	RAMOS	VIGOR
18.5	34	26
16.5	36	28
22.0	43	25
18.6	46	30
15.1	35	25
10.8	33	22
32.4	52	32
13.6	32	24
17.4	29	24
16.7	36	23
14.9	39	26

#### Diagramas de dispersión matriciales creados por "psych", "PerformanceAnalytics", "GGally", "coorplot"

library(psych)

pairs.panels(NECTARINA)

library(PerformanceAnalytics)

chart.Correlation(NECTARINA)

library(GGally)

ggpairs(NECTARINA)

library(corrplot)

corrplot(cor(NECTARINA), type="upper", order="hclust")

#### Correlación Simple
cor(NECTARINA, method = "pearson") # Valor de Correlación Múltiple es tan grande como el mayor de los coeficientes simples

##           PRODU     RAMOS     VIGOR
## PRODU 1.0000000 0.7929390 0.7493343
## RAMOS 0.7929390 1.0000000 0.8102622
## VIGOR 0.7493343 0.8102622 1.0000000

#Correlación más probabilidades
library(Hmisc)

rcorr(as.matrix(NECTARINA), type="pearson")

##       PRODU RAMOS VIGOR
## PRODU  1.00  0.79  0.75
## RAMOS  0.79  1.00  0.81
## VIGOR  0.75  0.81  1.00
## 
## n= 11 
## 
## 
## P
##       PRODU  RAMOS  VIGOR 
## PRODU        0.0036 0.0079
## RAMOS 0.0036        0.0025
## VIGOR 0.0079 0.0025

#### Correlación Parcial
# Opción 1
library(ggm)

pcor(c("PRODU", "RAMOS", "VIGOR"), var(NECTARINA)) # Correlación Parcial entre PRODUCCION y RAMOS, controlando para VIGOR

## [1] 0.478709

pcor(c("PRODU", "VIGOR", "RAMOS"), var(NECTARINA)) # Correlación Parcial entre PRODUCCION y VIGOR, controlando para RAMOS

## [1] 0.299211

# Opción 2
library(ppcor)

pcor(NECTARINA)

## $estimate
##          PRODU    RAMOS    VIGOR
## PRODU 1.000000 0.478709 0.299211
## RAMOS 0.478709 1.000000 0.535563
## VIGOR 0.299211 0.535563 1.000000
## 
## $p.value
##           PRODU     RAMOS     VIGOR
## PRODU 0.0000000 0.1616020 0.4009931
## RAMOS 0.1616020 0.0000000 0.1106067
## VIGOR 0.4009931 0.1106067 0.0000000
## 
## $statistic
##           PRODU    RAMOS     VIGOR
## PRODU 0.0000000 1.542180 0.8869295
## RAMOS 1.5421800 0.000000 1.7937352
## VIGOR 0.8869295 1.793735 0.0000000
## 
## $n
## [1] 11
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"

# Opción 3
# Cargar función "pcor.test". 
# Debe estar dentro del directorio de trabajo, en este caso está en "C:/Users/Administrator/Desktop"
source("pcor.R")
pcor.test(PRODU, RAMOS, VIGOR)

##   estimate   p.value statistic  n gn  Method            Use
## 1 0.478709 0.1230299   1.54218 11  1 Pearson Var-Cov matrix

pcor.test(PRODU, VIGOR, RAMOS)

##   estimate   p.value statistic  n gn  Method            Use
## 1 0.299211 0.3751168 0.8869295 11  1 Pearson Var-Cov matrix

#### Regresión múltiple
fit <- lm(PRODU ~ RAMOS + VIGOR, data = NECTARINA)
anova(fit); summary(fit)

## Analysis of Variance Table
## 
## Response: PRODU
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## RAMOS      1 199.494 199.494 14.8812 0.004825 **
## VIGOR      1  10.546  10.546  0.7866 0.400993   
## Residuals  8 107.246  13.406                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Call:
## lm(formula = PRODU ~ RAMOS + VIGOR, data = NECTARINA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3565 -2.2359 -0.5819  2.2771  4.5864 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.1376    10.3861  -1.361    0.211
## RAMOS         0.4491     0.2912   1.542    0.162
## VIGOR         0.5811     0.6552   0.887    0.401
## 
## Residual standard error: 3.661 on 8 degrees of freedom
## Multiple R-squared:  0.662,  Adjusted R-squared:  0.5775 
## F-statistic: 7.834 on 2 and 8 DF,  p-value: 0.01305

# Diagnósticos (Opción 1)
library(broom)

tidy(fit); glance(fit)

##          term    estimate  std.error  statistic   p.value
## 1 (Intercept) -14.1375577 10.3860693 -1.3612039 0.2105495
## 2       RAMOS   0.4491260  0.2912280  1.5421800 0.1616020
## 3       VIGOR   0.5811434  0.6552306  0.8869295 0.4009931

##   r.squared adj.r.squared    sigma statistic    p.value df    logLik
## 1 0.6619891     0.5774863 3.661386  7.833937 0.01305338  3 -28.13309
##        AIC      BIC deviance df.residual
## 1 64.26617 65.85775  107.246           8

knitr::kable(augment(fit))

PRODU	RAMOS	VIGOR	.fitted	.se.fit	.resid	.hat	.sigma	.cooksd	.std.resid
18.5	34	26	16.24245	1.582818	2.2575474	0.1868835	3.798084	0.0358201	0.6837779
16.5	36	28	18.30299	2.113180	-1.8029913	0.3331056	3.824198	0.0605401	-0.6030032
22.0	43	25	19.70344	2.326746	2.2965570	0.4038378	3.749268	0.1490117	0.8123619
18.6	46	30	23.95654	1.934662	-5.3565377	0.2792023	3.103900	0.3833968	-1.7231854
15.1	35	25	16.11044	1.199032	-1.0104352	0.1072434	3.893260	0.0034159	-0.2920764
10.8	33	22	13.46875	1.990027	-2.6687532	0.2954113	3.725157	0.1053805	-0.8683506
32.4	52	32	27.81358	2.746127	4.5864198	0.5625360	2.907169	1.5374502	1.8939006
13.6	32	24	14.18191	1.477964	-0.5819139	0.1629434	3.906797	0.0019581	-0.1737144
17.4	29	24	12.83454	2.022662	4.5654640	0.3051796	3.321953	0.3276183	1.4959022
16.7	36	23	15.39727	1.884472	1.3027255	0.2649041	3.871827	0.0206868	0.4149879
14.9	39	26	18.48808	1.150590	-3.5880824	0.0987530	3.644192	0.0389202	-1.0322741

# Diagnósticos (Opción 2)
interv <- predict(fit, interval = "confidence")
pronos <- predict(fit, interval = "prediction")

DIAGNOSTICOS <- data.frame(PRODU, interv, pronos, fit$resid, stdres(fit), cooks.distance(fit), studres(fit), dfbetas(fit))
knitr::kable(DIAGNOSTICOS)

PRODU	fit	lwr	upr	fit.1	lwr.1	upr.1	fit.resid	stdres.fit.	cooks.distance.fit.	studres.fit.	X.Intercept.	RAMOS	VIGOR
18.5	16.24245	12.592468	19.89244	16.24245	7.044110	25.44080	2.2575474	0.6837779	0.0358201	0.6591678	-0.0435785	-0.2263558	0.1874924
16.5	18.30299	13.429990	23.17599	18.30299	8.554491	28.05149	-1.8029913	-0.6030032	0.0605401	-0.5773307	0.2088718	0.3114692	-0.3432323
22.0	19.70344	14.337958	25.06893	19.70344	9.699666	29.70722	2.2965570	0.8123619	0.1490117	0.7933202	0.2777194	0.5663529	-0.5163084
18.6	23.95654	19.495200	28.41787	23.95654	14.407158	33.50592	-5.3565377	-1.7231854	0.3833968	-2.0326835	0.8659652	-0.1552059	-0.4762846
15.1	16.11044	13.345463	18.87541	16.11044	7.226056	24.99481	-1.0104352	-0.2920764	0.0034159	-0.2746809	-0.0292751	0.0247419	-0.0038027
10.8	13.46875	8.879743	18.05776	13.46875	3.859063	23.07844	-2.6687532	-0.8683506	0.1053805	-0.8534851	-0.4836367	-0.1940198	0.4015206
32.4	27.81358	21.481000	34.14616	27.81358	17.259496	38.36766	4.5864198	1.8939006	1.5374502	2.3852413	-1.8490508	0.9090149	0.6136257
13.6	14.18191	10.773723	17.59011	14.18191	5.076811	23.28702	-0.5819139	-0.1737144	0.0019581	-0.1628023	-0.0313574	0.0318034	-0.0048880
17.4	12.83454	8.170270	17.49880	12.83454	3.188682	22.48039	4.5654640	1.4959022	0.3276183	1.6487513	0.2227642	-0.8255010	0.4367646
16.7	15.39727	11.051673	19.74288	15.39727	5.901414	24.89314	1.3027255	0.4149879	0.0206868	0.3924326	0.1831469	0.1301902	-0.1873336
14.9	18.48808	15.834817	21.14135	18.48808	9.637832	27.33833	-3.5880824	-1.0322741	0.0389202	-1.0371446	-0.0506718	-0.0961950	0.0718384

# Gráficos de diagnósticos
par(mfrow = c(1, 1))
plot(fit, 1) # Residuos vs. Estimados

plot(fit, 2) # Normal Q-Q

plot(fit, 3) # Residuos estandarizados vs Estimados

plot(fit, 4) # Cook's distance

plot(fit, 5) # Residuos vs. leverage

plot(fit, 6) # Cook's distance vs. leverage

#### Estimando el "mejor" modelo, usando los procedimientos "Forward", "Backward", "Stepwise"
# Opción 1
library(MASS)
fit  <- lm(PRODU ~ RAMOS + VIGOR, data = NECTARINA)
step <- stepAIC(fit, direction="both")

## Start:  AIC=31.05
## PRODU ~ RAMOS + VIGOR
## 
##         Df Sum of Sq    RSS    AIC
## - VIGOR  1    10.546 117.79 30.081
## <none>               107.25 31.050
## - RAMOS  1    31.883 139.13 31.913
## 
## Step:  AIC=30.08
## PRODU ~ RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## <none>               117.79 30.081
## + VIGOR  1    10.546 107.25 31.050
## - RAMOS  1   199.494 317.29 38.981

step$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## PRODU ~ RAMOS + VIGOR
## 
## Final Model:
## PRODU ~ RAMOS
## 
## 
##      Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                             8   107.2460 31.04952
## 2 - VIGOR  1 10.54555         9   117.7915 30.08123

# Opción 2
step(lm(PRODU ~ RAMOS + VIGOR, data = NECTARINA), direction = "forward")

## Start:  AIC=31.05
## PRODU ~ RAMOS + VIGOR

## 
## Call:
## lm(formula = PRODU ~ RAMOS + VIGOR, data = NECTARINA)
## 
## Coefficients:
## (Intercept)        RAMOS        VIGOR  
##    -14.1376       0.4491       0.5811

step(lm(PRODU ~ RAMOS + VIGOR, data = NECTARINA), direction = "backward")

## Start:  AIC=31.05
## PRODU ~ RAMOS + VIGOR
## 
##         Df Sum of Sq    RSS    AIC
## - VIGOR  1    10.546 117.79 30.081
## <none>               107.25 31.050
## - RAMOS  1    31.883 139.13 31.913
## 
## Step:  AIC=30.08
## PRODU ~ RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## <none>               117.79 30.081
## - RAMOS  1    199.49 317.29 38.981

## 
## Call:
## lm(formula = PRODU ~ RAMOS, data = NECTARINA)
## 
## Coefficients:
## (Intercept)        RAMOS  
##     -6.9766       0.6584

step(lm(PRODU ~ VIGOR + RAMOS, data = NECTARINA), direction = "forward")

## Start:  AIC=31.05
## PRODU ~ VIGOR + RAMOS

## 
## Call:
## lm(formula = PRODU ~ VIGOR + RAMOS, data = NECTARINA)
## 
## Coefficients:
## (Intercept)        VIGOR        RAMOS  
##    -14.1376       0.5811       0.4491

step(lm(PRODU ~ VIGOR + RAMOS, data = NECTARINA), direction = "backward")

## Start:  AIC=31.05
## PRODU ~ VIGOR + RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## - VIGOR  1    10.546 117.79 30.081
## <none>               107.25 31.050
## - RAMOS  1    31.883 139.13 31.913
## 
## Step:  AIC=30.08
## PRODU ~ RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## <none>               117.79 30.081
## - RAMOS  1    199.49 317.29 38.981

## 
## Call:
## lm(formula = PRODU ~ RAMOS, data = NECTARINA)
## 
## Coefficients:
## (Intercept)        RAMOS  
##     -6.9766       0.6584

step(lm(PRODU ~ VIGOR + RAMOS, data = NECTARINA), direction = "both")

## Start:  AIC=31.05
## PRODU ~ VIGOR + RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## - VIGOR  1    10.546 117.79 30.081
## <none>               107.25 31.050
## - RAMOS  1    31.883 139.13 31.913
## 
## Step:  AIC=30.08
## PRODU ~ RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## <none>               117.79 30.081
## + VIGOR  1    10.546 107.25 31.050
## - RAMOS  1   199.494 317.29 38.981

## 
## Call:
## lm(formula = PRODU ~ RAMOS, data = NECTARINA)
## 
## Coefficients:
## (Intercept)        RAMOS  
##     -6.9766       0.6584

# Opción 3
fwd.modelo = step(lm(PRODU ~ 1, data=NECTARINA), direction  = 'forward',  scope =~ RAMOS + VIGOR)

## Start:  AIC=38.98
## PRODU ~ 1
## 
##         Df Sum of Sq    RSS    AIC
## + RAMOS  1    199.49 117.79 30.081
## + VIGOR  1    178.16 139.13 31.913
## <none>               317.29 38.981
## 
## Step:  AIC=30.08
## PRODU ~ RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## <none>               117.79 30.081
## + VIGOR  1    10.546 107.25 31.049

summary(fwd.modelo)

## 
## Call:
## lm(formula = PRODU ~ RAMOS, data = NECTARINA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7105 -2.3848 -0.2264  1.8776  5.2825 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -6.9766     6.4553  -1.081   0.3079   
## RAMOS         0.6584     0.1686   3.904   0.0036 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.618 on 9 degrees of freedom
## Multiple R-squared:  0.6288, Adjusted R-squared:  0.5875 
## F-statistic: 15.24 on 1 and 9 DF,  p-value: 0.003596

fwd.modelo$anova

##      Step Df Deviance Resid. Df Resid. Dev      AIC
## 1         NA       NA        10   317.2855 38.98097
## 2 + RAMOS -1  199.494         9   117.7915 30.08123

bwd.modelo = step(lm(PRODU ~ 1, data=NECTARINA), direction  = 'backward', scope =~ RAMOS + VIGOR)

## Start:  AIC=38.98
## PRODU ~ 1

summary(bwd.modelo)

## 
## Call:
## lm(formula = PRODU ~ 1, data = NECTARINA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0636 -2.8636 -1.1636  0.6864 14.5364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   17.864      1.698   10.52 9.99e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.633 on 10 degrees of freedom

bwd.modelo$anova

##   Step Df Deviance Resid. Df Resid. Dev      AIC
## 1      NA       NA        10   317.2855 38.98097

both.modelo = step(lm(PRODU ~ 1, data=NECTARINA), direction = 'both',     scope =~ RAMOS + VIGOR)

## Start:  AIC=38.98
## PRODU ~ 1
## 
##         Df Sum of Sq    RSS    AIC
## + RAMOS  1    199.49 117.79 30.081
## + VIGOR  1    178.16 139.13 31.913
## <none>               317.29 38.981
## 
## Step:  AIC=30.08
## PRODU ~ RAMOS
## 
##         Df Sum of Sq    RSS    AIC
## <none>               117.79 30.081
## + VIGOR  1    10.546 107.25 31.050
## - RAMOS  1   199.494 317.29 38.981

summary(both.modelo)

## 
## Call:
## lm(formula = PRODU ~ RAMOS, data = NECTARINA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7105 -2.3848 -0.2264  1.8776  5.2825 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -6.9766     6.4553  -1.081   0.3079   
## RAMOS         0.6584     0.1686   3.904   0.0036 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.618 on 9 degrees of freedom
## Multiple R-squared:  0.6288, Adjusted R-squared:  0.5875 
## F-statistic: 15.24 on 1 and 9 DF,  p-value: 0.003596

both.modelo$anova

##      Step Df Deviance Resid. Df Resid. Dev      AIC
## 1         NA       NA        10   317.2855 38.98097
## 2 + RAMOS -1  199.494         9   117.7915 30.08123

# Opción 4
library(leaps)

saltos <- regsubsets(PRODU ~ RAMOS + VIGOR, data = NECTARINA, nbest=10)
summary(saltos)

## Subset selection object
## Call: regsubsets.formula(PRODU ~ RAMOS + VIGOR, data = NECTARINA, nbest = 10)
## 2 Variables  (and intercept)
##       Forced in Forced out
## RAMOS     FALSE      FALSE
## VIGOR     FALSE      FALSE
## 10 subsets of each size up to 2
## Selection Algorithm: exhaustive
##          RAMOS VIGOR
## 1  ( 1 ) "*"   " "  
## 1  ( 2 ) " "   "*"  
## 2  ( 1 ) "*"   "*"

plot(saltos, scale = "r2")

plot(saltos, scale = "Cp")

plot(saltos, scale = "adjr2")

plot(saltos, scale = "bic")

Códigos para el software SAS

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

EXPERIMENTACIÓN EN AGRICULTURA

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

CAPÍTULO 16: REGRESIÓN MÚLTIPLE

REGRESIÓN CON DOS VARIABLES INDEPENDIENTES

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

DATA NECTARINA;

INPUT RAMA PRODU RAMOS VIGOR;

LABEL

RAMA = "Rama No."

PRODU = "Producción (kg/rama)"

RAMOS = "No. ramos fructíferos/rama"

VIGOR = "Perímetro de la rama, cm";

CARDS;

1 18.5 34 26

2 16.5 36 28

3 22.0 43 25

4 18.6 46 30

5 15.1 35 25

6 10.8 33 22

7 32.4 52 32

8 13.6 32 24

9 17.4 29 24

10 16.7 36 23

11 14.9 39 26

;

ODS HTML;

* Mostrar datos a analizar;

PROC PRINT DATA = NECTARINA;

RUN;

* Matrix de dispersión;

PROC INSIGHT DATA = NECTARINA;

SCATTER PRODU RAMOS VIGOR * PRODU RAMOS VIGOR;

RUN;

* Correlación Simple;

PROC CORR PEARSON DATA = NECTARINA;

VAR PRODU RAMOS VIGOR; *Valor de Correlación Múltiple es tan grande como el mayor de los coeficientes simples;

RUN;

* Correlación Parcial PRODU-RAMOS.VIGOR;

PROC CORR DATA = NECTARINA;

VAR RAMOS;

WITH PRODU;

PARTIAL VIGOR;

RUN;

* Correlación Parcial PRODU-VIGOR.RAMOS;

PROC CORR DATA = NECTARINA;

VAR VIGOR;

WITH PRODU;

PARTIAL RAMOS;

RUN;

* Regresión múltiple;

PROC REG;

MODEL PRODU = RAMOS VIGOR/ R P CLM CLI influence;

RUN;

* Mallows' Cp;

PROC REG;

MODEL PRODU = RAMOS VIGOR/ SELECTION = CP;

RUN;

* Regresión múltiple con el procedimiento FORWARD;

PROC REG;

MODEL PRODU = RAMOS VIGOR/ SELECTION = FORWARD SLENTRY=0.05 DETAILS;

RUN;

* Regresión múltiple con el procedimiento STEPWISE;

PROC REG;

MODEL PRODU = RAMOS VIGOR/ SELECTION = STEPWISE;

RUN;

* Regresión múltiple con el procedimiento BACKWARD;

PROC REG;

MODEL PRODU = RAMOS VIGOR/ SELECTION = BACKWARD;

RUN;

QUIT;

ODS HTML CLOSE;

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

Archivos adjuntos: capítulo 16.r, pcor.r, capítulo 16.sas
Nota: El procedimiento “Stepwise” está en capilla ardiente. Las razones están enumeradas en esta publicación.

Estudiando Agricultura en Taiwán

Páginas

mayo 23, 2017

Regresión con 2 variables independientes

2 comentarios: