Tarea No.1

La base leche.Rdata contiene informacion de 60 vacas de una finca ubicada en Zarcero, y se busca establecer las relaciones que existen entre ciertas variables de la vaca y la producción de leche. A continuación, se le da una descripción de las variables que se encuentran en la base en el orden que aparecen:

library(readxl)

## Warning: package 'readxl' was built under R version 4.3.3

base <- read_excel("leche.xlsx")

#Visualizacion de los datos

#que tenemos
str(base)

## tibble [60 × 9] (S3: tbl_df/tbl/data.frame)
##  $ ID        : num [1:60] 100 120 135 139 150 152 155 160 161 171 ...
##  $ Vaca      : chr [1:60] "T54" "U36" "V16" "V22" ...
##  $ produccion: num [1:60] 14.2 21.4 9.2 42.4 24.4 14.6 16.2 22.2 23.6 15.4 ...
##  $ lactancia : num [1:60] 419 143 412 247 118 429 126 65 138 300 ...
##  $ partos    : num [1:60] 9 9 7 7 7 7 6 8 8 6 ...
##  $ peso      : num [1:60] 395 510 370 425 450 385 470 425 455 375 ...
##  $ alimento  : num [1:60] 2.5 3.5 1.5 3 4 2.5 3 3.5 4 2.5 ...
##  $ edad      : num [1:60] 11.9 11.2 10.5 10.4 10.2 ...
##  $ especie   : chr [1:60] "raza2" "raza2" "raza1" "raza2" ...

producción - Cantidad de leche producida por una vaca en un día (kilogramos). Aunque la leche por lo general se mide en litros, los instrumentos no están graduados en litros sino que una vez depositada la leche en el molde, se procede a pesar directamente en la bascula.
lactancia - Número de días transcurridos desde el último parto de la vaca hasta la fecha de medición.
partos - Número total de partos que ha tenido la vaca.
peso - Peso de la vaca (kilogramos).
alimento - Cantidad de concentrado para producción lechera que se mezcla con un suplemento y miel, y que se le da a las vacas en cada ordeño (gramos).
edad - Edad de la vaca (años).
especie - Raza de la vaca; esta variable esta codificada como raza1 y raza2.

Lea el archivo de datos que está en modo xlsx. Ya que el archivo original viene en formato de excel, use la funcion read_excel de la librería readxl. Puede aprovechar para cambiar el formato de la variable especie de chr o caracter a factor, mediante la función factor(). Esto con el fin de definir la variable en un formato apropiado para sus respectivos análisis.

base$especie= factor(base$especie)
levels(base$especie)= c("raza1", "raza2")

Sálvelo en un nuevo archivo de tipo Rdata que se llame leche.Rdata usando la función save de la siguiente forma: save(base,file="leche.Rdata").

leche.Rdata = save(base, file = "leche.Rdata")

Defina la unidad de observación.

Y = base$produccion

Haga una nueva base que contenga sólo los datos de vacas que son raza2.

base1 = subset(base, especie == "raza2")
base1

## # A tibble: 34 × 9
##       ID Vaca  produccion lactancia partos  peso alimento  edad especie
##    <dbl> <chr>      <dbl>     <dbl>  <dbl> <dbl>    <dbl> <dbl> <fct>  
##  1   100 T54         14.2       419      9   395      2.5 12.0  raza2  
##  2   120 U36         21.4       143      9   510      3.5 11.2  raza2  
##  3   139 V22         42.4       247      7   425      3   10.4  raza2  
##  4   155 V53         16.2       126      6   470      3    9.97 raza2  
##  5   161 V64         23.6       138      8   455      4    9.82 raza2  
##  6   182 X57         21.4        92      7   420      3.5  8.88 raza2  
##  7   185 Y03         16.2       217      7   420      2.5  8.76 raza2  
##  8   217 Y30         26.8       127      5   415      4.5  7.83 raza2  
##  9   220 Z02         15.2       414      5   365      3    7.81 raza2  
## 10   259 Z36         20.4         4      6   510      3.5  6.96 raza2  
## # ℹ 24 more rows

¿Cuántas unidades contiene esta base?

nrow(base1)

## [1] 34

#####A PARTIR DE AHORA TRABAJE SOLO CON ESTA BASE

Familiarícese con las variables, sus unidades de medida, su rango, media y desviación estándar. Asegúrese de comprender el significado de cada variable y su construcción. Haga un análisis descriptivo de cada variable de interés (lactancia, partos, peso, alimento y edad) en su relación con la respuesta (produccion). Detecte valores extremos (no los elimine por ahora). Comente sobre si se dibuja algún tipo de relación clara entre la respuesta y cada predictor.

library(car)

## Warning: package 'car' was built under R version 4.3.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.3

varsindep = base1[,-c(1,2,9)]
sapply(varsindep,sd)

## produccion  lactancia     partos       peso   alimento       edad 
##  7.0645837 94.1976910  2.6171313 60.3838656  0.9005841  2.7922339

sapply(varsindep,mean)

## produccion  lactancia     partos       peso   alimento       edad 
##  17.288235 183.352941   3.382353 425.029412   2.852941   5.577647

sapply(varsindep, range)

##      produccion lactancia partos peso alimento  edad
## [1,]        6.8         4      1  310      1.5  2.43
## [2,]       42.4       419      9  565      4.5 11.95

scatterplotMatrix(varsindep)

Analisis descriptivo de las variables Producción: Esta variable está en un rango de 6.8 y 42.4. Representa la cantidad de leche producida en un día, su media es de 42.4 kologramos, la desviacion estandar de 17.28, una vaca tiene una alta produccion de leche en un día

Así con als demás

Elimine el valor extremo de la variable producción de leche y vuelva a revisar la relación entre la respuesta y los predictores.

base3= subset(varsindep, varsindep$produccion<30)

#####**A PARTIR DE AHORA USE LA BASE DONDE SE ELIMINARON LOS VALORES EXTREMOS

Obtenga los coeficientes de correlación entre la respuesta y cada predictor. Comente sobre la fuerza de la relación lineal que existe.

X= base3[,2:6]
Y= base3$produccion
R= cor(X,Y)
R

##                 [,1]
## lactancia -0.5036466
## partos     0.4075940
## peso       0.2093497
## alimento   0.9769086
## edad       0.3730765

Se puede observar que la fuerza que hay entre lactancia y produccion es fuerte y negativa, con partos,peso,alimetno y edad existe una relacion no muy fuerte y positiva y con el alimento, es una correlacion casi perfecta de 1, o sea, están muy correlacionas.

Estimación de los coeficientes del modelo.

Escriba un modelo lineal con los 4 predictores.

\[ \mu_{y\mid x_1,x_2,x_3,x_4,x_5}= \beta_0 + \beta_1*X_1 + \beta_2*X_2 + \beta_3*X_3+ \beta_4*X_4 +\beta_5*X_5 \\ \mu_{y\mid L,PA,PE,AL,E}= \beta_0 + \beta_1*L + \beta_2*PA + \beta_3*PE + \beta_4*AL +\beta_5*E \]

Cree un vector Y que contenga la variable respuesta.

Y= c(base3$produccion)
Y

##  [1] 14.2 21.4 16.2 23.6 21.4 16.2 26.8 15.2 20.4 25.0 21.6  9.8  9.6 26.4 23.8
## [16]  7.8 18.2 15.8 16.8 11.0 17.2 13.6 12.4 17.6 23.4 15.4  9.6  8.8  6.8 14.4
## [31] 18.2 10.6 16.2

Construya una matriz de estructura X que contenga los predictores de interés.

#MATRIZ DE ESTRUCTURA
matrizest = cbind(1, base3$lactancia,base3$partos,base3$peso,base3$alimento,base3$edad)
Y = base3$produccion
matrizest = as.matrix(matrizest)
str(matrizest)

##  num [1:33, 1:6] 1 1 1 1 1 1 1 1 1 1 ...

Obtenga los coeficientes estimados (betas) usando la matriz X y el vector Y.

\[ \hat{\beta}= (X^TX)^{-1}(X^TY) \]

#estimacion de coed por la matriz
beta= solve(t(matrizest) %*% matrizest) %*% t(matrizest) %*% Y
round(beta,3)

##        [,1]
## [1,]  2.861
## [2,] -0.003
## [3,]  0.344
## [4,] -0.005
## [5,]  5.750
## [6,] -0.196

Escriba la ecuación de regresión estimada.

\[ \mu_{y\mid x_1,x_2,x_3,x_4,x_5}= 2.861 -0.003*X_1 + 0.344*X_2 -0.005*X_3 + 5.750*X_4 -0.196 *X_5 \]

Interprete los coeficientes de la ecuación de regresión. Use unidades adecuadas de incremento en cada predictor. Previamente obtenga las desviaciones estándar de cada predictor para decidir el incremento. También tome en cuenta que la desviación estándar se ve afectada por valores extremos, entonces haga un boxplot de cada predictor para decidir si debe hacer un ajuste en su incremento.

sdlactancia = sd(base3$lactancia)
sdpartos = sd(base3$partos)
sdpeso = sd(base3$peso)
sdalimento = sd(base3$alimento)
sdedad = sd(base3$edad)
sdlactancia

## [1] 94.97402

sdpartos

## [1] 2.577217

sdpeso

## [1] 61.32011

sdalimento

## [1] 0.9141667

sdedad

## [1] 2.700846

par(mfrow = c(3, 2))
boxplot(base3$lactancia, ylab= "Lactancia")
boxplot(base3$partos, ylab= "partos")
boxplot(base3$peso, ylab= "peso")
boxplot(base3$alimento, ylab= "alimento")
boxplot(base3$edad, ylab= "edad")

lactancia - Número de días transcurridos desde el último parto de la vaca hasta la fecha de medición.
partos - Número total de partos que ha tenido la vaca.
peso - Peso de la vaca (kilogramos).
alimento - Cantidad de concentrado para producción lechera que se mezcla con un suplemento y miel, y que se le da a las vacas en cada ordeño (gramos).
edad - Edad de la vaca (años).

Obtenemos valores extremos solo en la lactancia, son 3 los podemos quitar y calcular una njeva desviacion estandar

sdlactancia1 = sd(subset(base3$lactancia, base3$lactancia < 300 ))
sdlactancia1

## [1] 67.8189

Nuevamente cargo las sd

sds <- round(cbind(sdlactancia1,sdpartos,sdpeso,sdalimento,sdedad),4)
sds

##      sdlactancia1 sdpartos  sdpeso sdalimento sdedad
## [1,]      67.8189   2.5772 61.3201     0.9142 2.7008

beta

##             [,1]
## [1,]  2.86078616
## [2,] -0.00268651
## [3,]  0.34404640
## [4,] -0.00537331
## [5,]  5.74997696
## [6,] -0.19643293

#Para las interpretaciones de las desviaciones

round(beta[2]*70,2)

## [1] -0.19

round(beta[3]*3,2)

## [1] 1.03

round(beta[4]*60,2)

## [1] -0.32

round(beta[5]*1,2)

## [1] 5.75

round(beta[6]*3,2)

## [1] -0.59

Entonces ahora si, despues de haber quitado los valores extremos, podemos decir que en realidad solo estaba afectando a la variable lactancia. Vamos a usar en lactancia 70 dias, en partos 3, en peso 60kg, en alimentos 1g y en peso 3kg

-Coeficiente de lactancia: Por cada aumento de 70 dias trascurridos despues del partos hasta la fecha de modificación, se espera que la produccion de leche de vacas de raza2 va a disminuir en promedio 0.19kg , manteniendo constante las demás variables

-Coeficiente de partos: Por cada aumento de 3 partos que haya tenido la vaca, se espera que la produccion de leche en vacas de raza dos aumente un 1.03kg, manteniendo constante las demás variables

-Coeficiente de peso: Por cada aumento de 60kg en las vacas de raza2, se espera que en promedio la produccion de leche disminuya 0.32kg, manteniendo constantes las demás variables

-Coeficiente de Alimento: Por cada aumento de 1g en el alimento de las vacas de raza dos, se espera que en promedio la produccion de leche aumente 5.74kg, manteniendo constante las demás variables

-Coeficiente de Edad: Por cada aumento de 3 años en las vacas de raza2, se espera que en promedio la produccion de leche disminuya 0.59kg, manteniendo constantes las demás variables

Usando la ecuación de regresión, obtenga el valor estimado para: lactancia = 200, partos = 5, peso = 400, alimento = 5 y edad = 4.

yh = t(beta) %*% c(1,200,5,400,5,4)
yh

##          [,1]
## [1,] 29.85855

Interprete este resultado (no ahorre palabras, esto debe estar bien claro).

La produccion de leche para las vacas de raza dos, con una lactancia de 200 días transcurridos despues del partor hasta la fecha de medicion, con 5 partos en total, con un peso de 400kg, un alimentos de 5g y una edad de 4 años, en promedio es de 30kg

Obtenga el valor estimado para: lactancia = 200, partos = 5, peso = 400, alimento = 5 y edad = 8.

yh1 = t(beta) %*% c(1,200,5,400,5,8)
yh1

##          [,1]
## [1,] 29.07281

¿Qué relación tiene este resultado con el coeficiente de edad?

Que cuando la edad es de 8 años, si mantenemos constantes las demás variables la produccion de leche va a disminuir en promedio 1.57kg

Obtenga el vector de valores estimados y el vector de residuales.

#valores estimados
ytechos = matrizest %*% beta
ytechos

##            [,1]
##  [1,] 14.736667
##  [2,] 20.753586
##  [3,] 17.352603
##  [4,] 23.868499
##  [5,] 21.145757
##  [6,] 15.083538
##  [7,] 26.346734
##  [8,] 17.223334
##  [9,] 20.931676
## [10,] 24.028588
## [11,] 19.712881
## [12,] 11.479987
## [13,]  8.226706
## [14,] 26.243981
## [15,] 23.458917
## [16,]  8.406938
## [17,] 19.741248
## [18,] 14.428182
## [19,] 16.390954
## [20,] 11.149210
## [21,] 17.452338
## [22,] 11.224316
## [23,] 11.628672
## [24,] 17.647617
## [25,] 22.783619
## [26,] 14.174053
## [27,]  8.903542
## [28,]  8.572563
## [29,]  8.210915
## [30,] 14.524869
## [31,] 20.124914
## [32,] 11.893891
## [33,] 17.548703

#residuales

residuales = Y - ytechos
residuales

##              [,1]
##  [1,] -0.53666738
##  [2,]  0.64641356
##  [3,] -1.15260300
##  [4,] -0.26849924
##  [5,]  0.25424336
##  [6,]  1.11646211
##  [7,]  0.45326592
##  [8,] -2.02333446
##  [9,] -0.53167641
## [10,]  0.97141164
## [11,]  1.88711924
## [12,] -1.67998696
## [13,]  1.37329365
## [14,]  0.15601891
## [15,]  0.34108299
## [16,] -0.60693829
## [17,] -1.54124847
## [18,]  1.37181811
## [19,]  0.40904588
## [20,] -0.14921004
## [21,] -0.25233805
## [22,]  2.37568434
## [23,]  0.77132837
## [24,] -0.04761734
## [25,]  0.61638147
## [26,]  1.22594742
## [27,]  0.69645790
## [28,]  0.22743689
## [29,] -1.41091549
## [30,] -0.12486894
## [31,] -1.92491406
## [32,] -1.29389065
## [33,] -1.34870299

Obtenga la suma de cuadrados residual (SCRes) usando el vector de residuales.

SCRes = t(residuales) %*% residuales
SCRes

##          [,1]
## [1,] 40.09785

Obtenga la estimación insesgada de la variancia condicional de la respuesta.

n= nrow(base3)
p=length(beta)
n

## [1] 33

## [1] 6

CMRes = SCRes / (n-p)
CMRes

##          [,1]
## [1,] 1.485105

¿Qué pasa con la SCRes si el coeficiente de edad aumenta un poco? Verifíquelo cambiando un poquito el coeficiente de edad.

beta1 = beta
beta1[6] = beta1[6]+0.002
beta1

##             [,1]
## [1,]  2.86078616
## [2,] -0.00268651
## [3,]  0.34404640
## [4,] -0.00537331
## [5,]  5.74997696
## [6,] -0.19443293

ytechos1 = matrizest %*% beta1
ytechos1

##            [,1]
##  [1,] 14.760567
##  [2,] 20.776026
##  [3,] 17.372543
##  [4,] 23.888139
##  [5,] 21.163517
##  [6,] 15.101058
##  [7,] 26.362394
##  [8,] 17.238954
##  [9,] 20.945596
## [10,] 24.042408
## [11,] 19.724901
## [12,] 11.491407
## [13,]  8.237326
## [14,] 26.254201
## [15,] 23.468577
## [16,]  8.416378
## [17,] 19.750268
## [18,] 14.436802
## [19,] 16.399534
## [20,] 11.156910
## [21,] 17.459658
## [22,] 11.231496
## [23,] 11.635592
## [24,] 17.654397
## [25,] 22.789839
## [26,] 14.180253
## [27,]  8.909702
## [28,]  8.578783
## [29,]  8.217015
## [30,] 14.530889
## [31,] 20.130874
## [32,] 11.898931
## [33,] 17.553563

residuales1 = Y - ytechos1
residuales1

##              [,1]
##  [1,] -0.56056738
##  [2,]  0.62397356
##  [3,] -1.17254300
##  [4,] -0.28813924
##  [5,]  0.23648336
##  [6,]  1.09894211
##  [7,]  0.43760592
##  [8,] -2.03895446
##  [9,] -0.54559641
## [10,]  0.95759164
## [11,]  1.87509924
## [12,] -1.69140696
## [13,]  1.36267365
## [14,]  0.14579891
## [15,]  0.33142299
## [16,] -0.61637829
## [17,] -1.55026847
## [18,]  1.36319811
## [19,]  0.40046588
## [20,] -0.15691004
## [21,] -0.25965805
## [22,]  2.36850434
## [23,]  0.76440837
## [24,] -0.05439734
## [25,]  0.61016147
## [26,]  1.21974742
## [27,]  0.69029790
## [28,]  0.22121689
## [29,] -1.41701549
## [30,] -0.13088894
## [31,] -1.93087406
## [32,] -1.29893065
## [33,] -1.35356299

SCRes1 = t(residuales1) %*% residuales1
SCRes1

##          [,1]
## [1,] 40.10268

SCRes

##          [,1]
## [1,] 40.09785

Como se usa el metodo de minimos cuadrados, ya se sabe que para la estimacion de los betas, se escontró el menos SCRes,entonces un aumento en la edad, solova a provocar que la SCRes aumente

Obtenga el coeficiente de determinación e interprételo.

SCTotal = var(base3$produccion) * (n-1)

#Calcilo de SCRegresion

#paso1
ybarra = mean(base3$produccion)
SCReg= sum((ytechos-ybarra)^2)
SCReg

## [1] 957.1676

R= SCReg/SCTotal
R

## [1] 0.9597922

El coeficiente de determinacion dio 0.96, lo que nos dicen que las variables explican un 96% de la variabilidad de nuestro modelo con respecto a la produccion de leche en vacas de raza 2

Obtenga la matriz de covariancias de los betas y los errores estándar de los betas.

#Matriz de varianza-covarianza
CMRes

##          [,1]
## [1,] 1.485105

matrizcov = c(CMRes)* solve(t(matrizest) %*% matrizest)
matrizcov

##              [,1]          [,2]          [,3]          [,4]          [,5]
## [1,]  4.641798238 -2.308944e-03  0.4257968760 -6.894041e-03 -2.547087e-01
## [2,] -0.002308944  1.103001e-05  0.0008510151  2.074753e-06  5.698895e-04
## [3,]  0.425796876  8.510151e-04  0.2705161550 -3.317544e-04  2.567340e-02
## [4,] -0.006894041  2.074753e-06 -0.0003317544  1.522261e-05 -7.084856e-06
## [5,] -0.254708718  5.698895e-04  0.0256733974 -7.084856e-06  9.636962e-02
## [6,] -0.352681801 -9.172735e-04 -0.2573064519  2.123567e-04 -3.759353e-02
##               [,6]
## [1,] -0.3526818013
## [2,] -0.0009172735
## [3,] -0.2573064519
## [4,]  0.0002123567
## [5,] -0.0375935290
## [6,]  0.2536936637

eebetas = sqrt(diag(matrizcov))
eebetas

## [1] 2.154483288 0.003321147 0.520111675 0.003901616 0.310434573 0.503680120

Construya los intervalos de confianza para los betas. Interprételos.

t= qt(0.975,  n-p)

ICBETA = cbind(beta-t*eebetas, beta+t*eebetas)
ICBETA

##             [,1]        [,2]
## [1,] -1.55984840 7.281420717
## [2,] -0.00950094 0.004127920
## [3,] -0.72313460 1.411227411
## [4,] -0.01337877 0.002632145
## [5,]  5.11301783 6.386936093
## [6,] -1.22989918 0.837033307

ICBETAALIMENTO = ICBETA[5,]*5
ICBETAALIMENTO

## [1] 25.56509 31.93468

Podemos interpretar la variables alimentos, se puede decir que con un 95% de confianza, el aumento en 5g de alimento para las vacas de raza 2, manteniendo constantes las demás variables, se espera que la produccion promedio de leche de las vacas de raza 2 aumente entre 26kg y 32kg

Ya que en la mayoria de los intervalos está contenido el cero (a excepción de la variable alimento), lo único que se puede decir es que al aumentar cualquiera de los predictores, no se puede demostrar que haya un incremento o disminución en el promedio de producción de leche. El coeficiente para el intercepto no se interpreta.**

Ajuste el modelo con 5 predictores utilizando la función lm y obtenga lo que se le pide en cada inciso. Compare los resultados con los obtenidos en las preguntas anteriores.

mod1= lm(produccion ~ lactancia+partos+peso+alimento+edad, data=base3)

Los coeficientes estimados (betas).

betas = mod1$coefficients
betas

## (Intercept)   lactancia      partos        peso    alimento        edad 
##  2.86078616 -0.00268651  0.34404640 -0.00537331  5.74997696 -0.19643293

El valor estimado para: lactancia = 200, partos = 5, peso = 400, alimento = 5 y edad = 4.

yhmod1= predict(mod1, data.frame(lactancia= 200,partos= 5,peso= 400,alimento= 5,edad= 4))

El vector de valores estimados y el vector de residuales.

ytechosmod1 = predict(mod1)
ytechosmod1

##         1         2         3         4         5         6         7         8 
## 14.736667 20.753586 17.352603 23.868499 21.145757 15.083538 26.346734 17.223334 
##         9        10        11        12        13        14        15        16 
## 20.931676 24.028588 19.712881 11.479987  8.226706 26.243981 23.458917  8.406938 
##        17        18        19        20        21        22        23        24 
## 19.741248 14.428182 16.390954 11.149210 17.452338 11.224316 11.628672 17.647617 
##        25        26        27        28        29        30        31        32 
## 22.783619 14.174053  8.903542  8.572563  8.210915 14.524869 20.124914 11.893891 
##        33 
## 17.548703

residualesmod1 = mod1$residuals
residualesmod1

##           1           2           3           4           5           6 
## -0.53666738  0.64641356 -1.15260300 -0.26849924  0.25424336  1.11646211 
##           7           8           9          10          11          12 
##  0.45326592 -2.02333446 -0.53167641  0.97141164  1.88711924 -1.67998696 
##          13          14          15          16          17          18 
##  1.37329365  0.15601891  0.34108299 -0.60693829 -1.54124847  1.37181811 
##          19          20          21          22          23          24 
##  0.40904588 -0.14921004 -0.25233805  2.37568434  0.77132837 -0.04761734 
##          25          26          27          28          29          30 
##  0.61638147  1.22594742  0.69645790  0.22743689 -1.41091549 -0.12486894 
##          31          32          33 
## -1.92491406 -1.29389065 -1.34870299

La suma de cuadrados residual y el cuadrado medio residual.

anovamod1= anova(mod1)
anovamod1

## Analysis of Variance Table
## 
## Response: produccion
##           Df Sum Sq Mean Sq  F value    Pr(>F)    
## lactancia  1 252.97  252.97 170.3356 3.548e-13 ***
## partos     1 164.18  164.18 110.5496 4.790e-11 ***
## peso       1   4.49    4.49   3.0254   0.09336 .  
## alimento   1 535.30  535.30 360.4489 < 2.2e-16 ***
## edad       1   0.23    0.23   0.1521   0.69960    
## Residuals 27  40.10    1.49                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SCResmod1 =anovamod1[6,2]
SCResmod1

## [1] 40.09785

La matriz de covariancia de los betas.

matrizvcovmod1= vcov(mod1)

Los errores estándar de los betas.

ee <- summary(mod1)$coef[,2]

Los grados de libertad de los residuales.

gl= summary(anovamod1)[6,1]
gl

## [1] "Max.   :27.000  "

El coeficiente de determinación.

summary(mod1)$r.sq

## [1] 0.9597922

Los intervalos de confianza de los betas.

confint(mod1)

##                   2.5 %      97.5 %
## (Intercept) -1.55984840 7.281420717
## lactancia   -0.00950094 0.004127920
## partos      -0.72313460 1.411227411
## peso        -0.01337877 0.002632145
## alimento     5.11301783 6.386936093
## edad        -1.22989918 0.837033307

Encuentre la variancia de la respuesta promedio con las siguientes características: lactancia = 200, partos = 5, peso = 400, alimento = 5 y edad = 4.

Yind = c(1,200,5,400,5,4)
varh= t(Yind) %*% matrizvcovmod1 %*% Yind
varh

##          [,1]
## [1,] 3.718463

Encuentre la variancia de los valores individuales de Y con las mismas características mencionadas en el punto anterior.

varindi = CMRes + varh
varindi

##          [,1]
## [1,] 5.203569

Obtenga un intervalo de confianza para los valores individuales de Y (hágalo paso a paso y con la función de R). Interprete el intervalo.

#manual
eevarindi= sqrt(varindi)
eevarindi

##          [,1]
## [1,] 2.281133

ICVindi = cbind(yhmod1-t*eevarindi,yhmod1+t*eevarindi)
ICVindi

##          [,1]     [,2]
## [1,] 25.17805 34.53904

#Automaticamente

ICVindiauto = predict(mod1, data.frame(lactancia= 200,partos= 5,peso= 400,alimento= 5,edad= 4),
                      interval = "prediction")
ICVindiauto

##        fit      lwr      upr
## 1 29.85855 25.17805 34.53904

Obtenga un intervalo de confianza para la respuesta promedio (hágalo paso a paso y con la función de R). Interprete el intervalo.

eemedia= sqrt(varh)
ICMEDIA = cbind(yhmod1-t*eemedia,yhmod1+t*eemedia)
ICMEDIA

##          [,1]     [,2]
## [1,] 25.90194 33.81516

ICVmedia = predict(mod1, data.frame(lactancia= 200,partos= 5,peso= 400,alimento= 5,edad= 4),
                      interval = "confidence")
ICVmedia

##        fit      lwr      upr
## 1 29.85855 25.90194 33.81516

Estime dos modelos: 1) con todos los predictores y 2) con sólo alimento como predictor.

mod2= lm(produccion~alimento,data=base3)

Compare los modelos usando el coeficiente de determinación.

summary(mod1)$r.sq

## [1] 0.9597922

summary(mod2)$r.sq

## [1] 0.9543503

Obtenga los coeficientes estandarizados para el modelo con los 5 predictores.

Hágalo usando los coeficientes sin estandarizar.

sdx =apply(base3[,-1],2,sd)
sdy = sd(base3$produccion)
sdx

##  lactancia     partos       peso   alimento       edad 
## 94.9740196  2.5772166 61.3201052  0.9141667  2.7008465

sdy

## [1] 5.582521

betasmod1 = mod1$coefficients
beta.s1 = betasmod1[-1]* sdx/sdy
beta.s1

##   lactancia      partos        peso    alimento        edad 
## -0.04570491  0.15883183 -0.05902207  0.94158841 -0.09503505

Hágalo estandarizando las variables.

base4= data.frame(scale(base3))
base4

##     produccion   lactancia     partos          peso   alimento       edad
## 1  -0.41688559  2.50148155  2.2222706 -0.4897301289 -0.3812049  2.4133848
## 2   0.87285420 -0.40457635  2.2222706  1.3856743507  0.7126875  2.1430992
## 3  -0.05862454 -0.58357267  1.0582241  0.7333597491  0.1657413  1.6802813
## 4   1.26694136 -0.45722233  1.8342551  0.4887417735  1.2596337  1.6247432
## 5   0.87285420 -0.94156531  1.4462396 -0.0820335029  0.7126875  1.2767041
## 6  -0.05862454  0.37458410  1.4462396 -0.0820335029 -0.3812049  1.2322736
## 7   1.84015904 -0.57304348  0.6702086 -0.1635728281  1.8065799  0.8879371
## 8  -0.23775506  2.44883557  0.6702086 -0.9789660801  0.1657413  0.8805320
## 9   0.69372367 -1.86813450  1.0582241  1.3856743507  0.7126875  0.5658159
## 10  1.51772409 -1.09950324  0.6702086 -0.8974267549  1.2596337  0.5473032
## 11  0.90868031 -0.21505084  0.2821931  2.2826069278  0.7126875  0.2140743
## 12 -1.20505990 -0.46775152  0.2821931  1.4672136759 -0.9281511  0.1029980
## 13 -1.24088601  1.02739421 -0.1058224  0.6518204239 -1.4750973 -0.0451037
## 14  1.76850683 -0.87839014 -0.1058224 -0.0820335029  1.8065799 -0.1191546
## 15  1.30276746 -1.12056163 -0.1058224 -0.0004941777  1.2596337 -0.2228258
## 16 -1.56332096  1.84867144 -0.4938379 -1.2235840557 -1.4750973 -0.2635537
## 17  0.29963652  0.23770456 -0.4938379  0.6518204239  0.7126875 -0.3413071
## 18 -0.13027675  0.15347100 -0.4938379 -0.4897301289 -0.3812049 -0.4153580
## 19  0.04885378  0.47987605 -0.4938379  2.0379889522  0.1657413 -0.4227631
## 20 -0.99010327  0.30087973 -0.4938379  0.8964383995 -0.9281511 -0.5856750
## 21  0.12050599 -0.35193038 -0.4938379 -0.1635728281  0.1657413 -0.6560233
## 22 -0.52436390  0.85892708 -0.8818534 -0.6528087793 -0.9281511 -0.6819411
## 23 -0.73932054  0.95368984 -0.8818534 -1.8758986573 -0.9281511 -0.7300741
## 24  0.19215820 -0.04658371 -0.8818534 -1.8758986573  0.1657413 -0.7559920
## 25  1.23111525 -0.03605452 -0.8818534  0.1462766076  1.2596337 -0.8596632
## 26 -0.20192896  0.10082502 -0.8818534 -0.0004941777 -0.3812049 -0.8633657
## 27 -1.24088601  0.09029583 -0.8818534 -1.4355863012 -1.4750973 -0.8707708
## 28 -1.38419043  0.03764985 -0.8818534 -0.4081908037 -1.4750973 -0.8596632
## 29 -1.74245148  0.86945628 -0.8818534  0.0810451475 -1.4750973 -0.8818784
## 30 -0.38105948 -0.46775152 -0.8818534 -0.5712694541 -0.3812049 -0.8966886
## 31  0.29963652 -0.38351796 -0.8818534 -0.1635728281  0.7126875 -0.9077962
## 32 -1.06175548 -1.04685727 -0.8818534 -0.5712694541 -0.9281511 -1.0781132
## 33 -0.05862454 -1.34167473 -0.8818534 -0.0004941777  0.1657413 -1.1114361

mod4 = lm(produccion~lactancia+partos+peso+alimento+edad, data=base4)
mod4

## 
## Call:
## lm(formula = produccion ~ lactancia + partos + peso + alimento + 
##     edad, data = base4)
## 
## Coefficients:
## (Intercept)    lactancia       partos         peso     alimento         edad  
##   1.414e-16   -4.570e-02    1.588e-01   -5.902e-02    9.416e-01   -9.504e-02

betaesta = mod4$coef

Use las matrices de correlaciones.

## # A tibble: 33 × 5
##    lactancia partos  peso alimento  edad
##        <dbl>  <dbl> <dbl>    <dbl> <dbl>
##  1       419      9   395      2.5 12.0 
##  2       143      9   510      3.5 11.2 
##  3       126      6   470      3    9.97
##  4       138      8   455      4    9.82
##  5        92      7   420      3.5  8.88
##  6       217      7   420      2.5  8.76
##  7       127      5   415      4.5  7.83
##  8       414      5   365      3    7.81
##  9         4      6   510      3.5  6.96
## 10        77      5   370      4    6.91
## # ℹ 23 more rows

##  [1] 14.2 21.4 16.2 23.6 21.4 16.2 26.8 15.2 20.4 25.0 21.6  9.8  9.6 26.4 23.8
## [16]  7.8 18.2 15.8 16.8 11.0 17.2 13.6 12.4 17.6 23.4 15.4  9.6  8.8  6.8 14.4
## [31] 18.2 10.6 16.2

RXY = cor(X,Y)
RXX = cor(X,X)
RXY

##                 [,1]
## lactancia -0.5036466
## partos     0.4075940
## peso       0.2093497
## alimento   0.9769086
## edad       0.3730765

RXX

##              lactancia       partos       peso   alimento      edad
## lactancia  1.000000000 -0.003679263 -0.2819211 -0.4928825 0.1042191
## partos    -0.003679263  1.000000000  0.3169245  0.3828466 0.9805259
## peso      -0.281921126  0.316924516  1.0000000  0.2445352 0.2641401
## alimento  -0.492882503  0.382846558  0.2445352  1.0000000 0.3533670
## edad       0.104219101  0.980525899  0.2641401  0.3533670 1.0000000

beta.s = solve(RXX) %*% RXY
beta.s

##                  [,1]
## lactancia -0.04570491
## partos     0.15883183
## peso      -0.05902207
## alimento   0.94158841
## edad      -0.09503505

beta.s1

##   lactancia      partos        peso    alimento        edad 
## -0.04570491  0.15883183 -0.05902207  0.94158841 -0.09503505

beta.s

##                  [,1]
## lactancia -0.04570491
## partos     0.15883183
## peso      -0.05902207
## alimento   0.94158841
## edad      -0.09503505

betaesta

##   (Intercept)     lactancia        partos          peso      alimento 
##  1.414396e-16 -4.570491e-02  1.588318e-01 -5.902207e-02  9.415884e-01 
##          edad 
## -9.503505e-02

Escriba la ecuación del modelo estimado con los coeficientes estandarizados.

\[ \hat{y}^s= -0.05*X_1 + 0.16*{X_2}^s -0.06*{X_3}^s+ 0.94*{X_4}^s -0.95*{X_5}^s \]

Interprete los coeficientes. Observe cuál variable tiene más impacto en la producción de leche.

MODELOS DE REGRESIÓN APLICADOS

Tarea No.1