Modelos Reducidos

Modelos Lineales

De: Luis Maron & Manuel Condori

Data:

Se tiene la información recopilada de medias biométricas: peso vivo(kg), largo de cuerpo (cm) y amplitud toraxica (cm) de las vicuñas machos En condiciones de pastoreo natural del departamento de Puno - Perú. 
dat<-("
PEVI  LCUERP    AMTORX
    31.5    74    76
    37      73    79
    46      75    89
    46      74    86
    42.5    77    77
    21      61    66
    46.5    75    90
    45.5    80    88
    44      80    86
    51.5    84    94
    38.5    77    81
    34.5    69    77
    46      75    87
    47.5    78    97
    28.5    71    72")     
data1 = read.table(textConnection(dat),header=TRUE)
library(sjPlot)
tab_df(data1)
PEVI LCUERP AMTORX
31.5 74 76
37 73 79
46 75 89
46 74 86
42.5 77 77
21 61 66
46.5 75 90
45.5 80 88
44 80 86
51.5 84 94
38.5 77 81
34.5 69 77
46 75 87
47.5 78 97
28.5 71 72

Modelo

attach(data1)
model1<-lm(PEVI~.,data=data1)
summary(model1)
## 
## Call:
## lm(formula = PEVI ~ ., data = data1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2180 -1.3433 -0.7693  1.1717  5.3393 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -51.7890    10.6982  -4.841 0.000404 ***
## LCUERP        0.4498     0.2245   2.004 0.068185 .  
## AMTORX        0.7054     0.1414   4.988 0.000316 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.836 on 12 degrees of freedom
## Multiple R-squared:  0.9031, Adjusted R-squared:  0.8869 
## F-statistic: 55.89 on 2 and 12 DF,  p-value: 8.303e-07

tab_model(model1)
  PEVI
Predictors Estimates CI p
(Intercept) -51.79 -75.10 – -28.48 <0.001
LCUERP 0.45 -0.04 – 0.94 0.068
AMTORX 0.71 0.40 – 1.01 <0.001
Observations 15
R2 / adjusted R2 0.903 / 0.887

Modelo:  Y(PEVI)= − 51.7890 + 0.4498 * LCUERP + 0.7054 * AMTORX

Tabla ANOVA

ANOVA<-aov(lm(PEVI~.,data1))
summary(ANOVA)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## LCUERP       1  698.9   698.9   86.90 7.61e-07 ***
## AMTORX       1  200.1   200.1   24.88 0.000316 ***
## Residuals   12   96.5     8.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

X <- model.matrix(model1); X
##    (Intercept) LCUERP AMTORX
## 1            1     74     76
## 2            1     73     79
## 3            1     75     89
## 4            1     74     86
## 5            1     77     77
## 6            1     61     66
## 7            1     75     90
## 8            1     80     88
## 9            1     80     86
## 10           1     84     94
## 11           1     77     81
## 12           1     69     77
## 13           1     75     87
## 14           1     78     97
## 15           1     71     72
## attr(,"assign")
## [1] 0 1 2

XT <- t(X); XT #X'
##              1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## (Intercept)  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
## LCUERP      74 73 75 74 77 61 75 80 80 84 77 69 75 78 71
## AMTORX      76 79 89 86 77 66 90 88 86 94 81 77 87 97 72
## attr(,"assign")
## [1] 0 1 2

XTX <- XT%*%X; XTX
##             (Intercept) LCUERP AMTORX
## (Intercept)          15   1123   1245
## LCUERP             1123  84477  93704
## AMTORX             1245  93704 104347

XTY <- XT%*%PEVI; XTY #X'Y
##                [,1]
## (Intercept)   606.5
## LCUERP      45936.5
## AMTORX      51276.0

BETAS<-solve(XTX)%*%XTY; BETAS
##                    [,1]
## (Intercept) -51.7889714
## LCUERP        0.4498109
## AMTORX        0.7053791

pp <- mean(PEVI); pp
## [1] 40.43333

SSM <- length(PEVI)*pp^2; SSM
## [1] 24522.82

SSR <- t(BETAS)%*%XTY; SSR
##          [,1]
## [1,] 25421.74

SSRm <- t(BETAS)%*%XTY- SSM ; SSRm 
##          [,1]
## [1,] 898.9273

SSTm <- PEVI%*%PEVI- SSM; SSTm
##          [,1]
## [1,] 995.4333

SSE <- SSTm- SSRm; SSE
##          [,1]
## [1,] 96.50604

MSE <- SSE/ANOVA$df.residual; MSE
##         [,1]
## [1,] 8.04217


$${\\}$$

Modelo reducido(Ejemplos)

I) Ha :  b0 = 0

modelo <- lm(PEVI~0+LCUERP+AMTORX)
X <- model.matrix(modelo); X
##    LCUERP AMTORX
## 1      74     76
## 2      73     79
## 3      75     89
## 4      74     86
## 5      77     77
## 6      61     66
## 7      75     90
## 8      80     88
## 9      80     86
## 10     84     94
## 11     77     81
## 12     69     77
## 13     75     87
## 14     78     97
## 15     71     72
## attr(,"assign")
## [1] 1 2

XT <- t(X); XT #X'
##         1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## LCUERP 74 73 75 74 77 61 75 80 80 84 77 69 75 78 71
## AMTORX 76 79 89 86 77 66 90 88 86 94 81 77 87 97 72
## attr(,"assign")
## [1] 1 2

XTX <- XT%*%X; XTX
##        LCUERP AMTORX
## LCUERP  84477  93704
## AMTORX  93704 104347

XTY <- XT%*%PEVI; XTY #X'Y
##           [,1]
## LCUERP 45936.5
## AMTORX 51276.0

BETAS_p<-solve(XTX)%*%XTY; BETAS_p
##              [,1]
## LCUERP -0.3314880
## AMTORX  0.7890764

pp <- mean(PEVI); pp
## [1] 40.43333

SSM <- length(PEVI)*pp^2; SSM
## [1] 24522.82

SSR_ <- t(BETAS_p)%*%XTY; SSR_
##          [,1]
## [1,] 25233.28

SSRm <- t(BETAS_p)%*%XTY- SSM ; SSRm 
##          [,1]
## [1,] 710.4646

SSTm <- PEVI%*%PEVI- SSM; SSTm
##          [,1]
## [1,] 995.4333

SSE <- SSTm- SSRm; SSE
##          [,1]
## [1,] 284.9688

MSE <- SSE/ANOVA$df.residual; MSE
##         [,1]
## [1,] 23.7474

Q <- SSR - SSR_; Q
##          [,1]
## [1,] 188.4627

tab_model(modelo)
  PEVI
Predictors Estimates CI p
LCUERP -0.33 -0.89 – 0.22 0.221
AMTORX 0.79 0.29 – 1.29 0.005
Observations 15
R2 / adjusted R2 0.989 / 0.987 
summary(aov(modelo))
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## LCUERP     1  24979   24979 1139.52 4.77e-14 ***
## AMTORX     1    254     254   11.59   0.0047 ** 
## Residuals 13    285      22                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo)
## 
## Call:
## lm(formula = PEVI ~ 0 + LCUERP + AMTORX)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.8583  -3.2853   0.3447   2.6194   7.2657 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)   
## LCUERP  -0.3315     0.2576  -1.287   0.2205   
## AMTORX   0.7891     0.2317   3.405   0.0047 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.682 on 13 degrees of freedom
## Multiple R-squared:  0.9888, Adjusted R-squared:  0.9871 
## F-statistic: 575.6 on 2 and 13 DF,  p-value: 2.05e-13

Modelo:  Y(PEVI)= − 0.3315 * LCUERP + 0.7891 * AMTORX

Decicion:

Tomamos el p-value para contranstar: Dado que p(0.00...205)<α = 0.05, rechazamos la H0. Por lo que  b0 ≠ 0.
$${\\}$$

II) Ha :  b0 = b1 = 0

modelo <- lm(PEVI~0+0+AMTORX)
X <- model.matrix(modelo); X
##    AMTORX
## 1      76
## 2      79
## 3      89
## 4      86
## 5      77
## 6      66
## 7      90
## 8      88
## 9      86
## 10     94
## 11     81
## 12     77
## 13     87
## 14     97
## 15     72
## attr(,"assign")
## [1] 1

XT <- t(X); XT #X'
##         1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## AMTORX 76 79 89 86 77 66 90 88 86 94 81 77 87 97 72
## attr(,"assign")
## [1] 1

XTX <- XT%*%X; XTX
##        AMTORX
## AMTORX 104347

XTY <- XT%*%PEVI; XTY #X'Y
##         [,1]
## AMTORX 51276

BETAS_p<-solve(XTX)%*%XTY; BETAS_p
##             [,1]
## AMTORX 0.4913989

pp <- mean(PEVI); pp
## [1] 40.43333

SSM <- length(PEVI)*pp^2; SSM
## [1] 24522.82

SSR_ <- t(BETAS_p)%*%XTY; SSR_
##          [,1]
## [1,] 25196.97

SSRm <- t(BETAS_p)%*%XTY- SSM ; SSRm 
##          [,1]
## [1,] 674.1528

SSTm <- PEVI%*%PEVI- SSM; SSTm
##          [,1]
## [1,] 995.4333

SSE <- SSTm- SSRm; SSE
##          [,1]
## [1,] 321.2805

MSE <- SSE/ANOVA$df.residual; MSE
##          [,1]
## [1,] 26.77338

Q <- SSR - SSR_; Q
##          [,1]
## [1,] 224.7745

f <- Q/MSE;f
##          [,1]
## [1,] 8.395447

tab_model(modelo)
  PEVI
Predictors Estimates CI p
AMTORX 0.49 0.46 – 0.52 <0.001
Observations 15
R2 / adjusted R2 0.987 / 0.987 
summary(aov(modelo))
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## AMTORX     1  25197   25197    1098 1.06e-14 ***
## Residuals 14    321      23                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo)
## 
## Call:
## lm(formula = PEVI ~ 0 + 0 + AMTORX)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.432  -2.579   1.740   2.761   5.309 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## AMTORX  0.49140    0.01483   33.14 1.06e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.79 on 14 degrees of freedom
## Multiple R-squared:  0.9874, Adjusted R-squared:  0.9865 
## F-statistic:  1098 on 1 and 14 DF,  p-value: 1.056e-14

Modelo:  Y(PEVI)=0.49140 * AMTORX

Decicion:

Tomamos el p-value para contranstar: Dado que p(0.00...106)<α = 0.05, rechazamos la H0. Por lo que  b0 ≠ b1 ≠ 0; esto quiere decir que al menos una de las Xs explican la variabilidad de Y.
$${\\}$$

III) Ha :  b0 = b2 = 0

modelo <- lm(PEVI~0+LCUERP+0)
X <- model.matrix(modelo); X
##    LCUERP
## 1      74
## 2      73
## 3      75
## 4      74
## 5      77
## 6      61
## 7      75
## 8      80
## 9      80
## 10     84
## 11     77
## 12     69
## 13     75
## 14     78
## 15     71
## attr(,"assign")
## [1] 1

XT <- t(X); XT #X'
##         1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## LCUERP 74 73 75 74 77 61 75 80 80 84 77 69 75 78 71
## attr(,"assign")
## [1] 1

XTX <- XT%*%X; XTX
##        LCUERP
## LCUERP  84477

XTY <- XT%*%PEVI; XTY #X'Y
##           [,1]
## LCUERP 45936.5

BETAS_p<-solve(XTX)%*%XTY; BETAS_p
##             [,1]
## LCUERP 0.5437752

pp <- mean(PEVI); pp
## [1] 40.43333

SSM <- length(PEVI)*pp^2; SSM
## [1] 24522.82

SSR_ <- t(BETAS_p)%*%XTY; SSR_
##          [,1]
## [1,] 24979.13

SSRm <- t(BETAS_p)%*%XTY- SSM ; SSRm 
##          [,1]
## [1,] 456.3141

SSTm <- PEVI%*%PEVI- SSM; SSTm
##          [,1]
## [1,] 995.4333

SSE <- SSTm- SSRm; SSE
##          [,1]
## [1,] 539.1192

MSE <- SSE/ANOVA$df.residual; MSE
##         [,1]
## [1,] 44.9266

Q <- SSR - SSR_; Q
##          [,1]
## [1,] 442.6132

f <- Q/MSE;f
##          [,1]
## [1,] 9.851917

tab_model(modelo)
  PEVI
Predictors Estimates CI p
LCUERP 0.54 0.50 – 0.59 <0.001
Observations 15
R2 / adjusted R2 0.979 / 0.977 
summary(aov(modelo))
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## LCUERP     1  24979   24979   648.7 3.97e-13 ***
## Residuals 14    539      39                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo)
## 
## Call:
## lm(formula = PEVI ~ 0 + LCUERP + 0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1703  -3.1956   0.6293   5.2169   5.8229 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## LCUERP  0.54378    0.02135   25.47 3.97e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.206 on 14 degrees of freedom
## Multiple R-squared:  0.9789, Adjusted R-squared:  0.9774 
## F-statistic: 648.7 on 1 and 14 DF,  p-value: 3.972e-13

Modelo:  Y(PEVI)=0.54378 * LCUERP

Decicion:

Tomamos el p-value para contranstar: Dado que p(0.00...397)<α = 0.05, rechazamos la H0. Por lo que  b0 ≠ b1 ≠ 0; esto quiere decir que al menos una de las Xs explican la variabilidad de Y.
$${\\}$$

Fin ☺