Modelos gamlss aplicados a la producción de la lombriz roja californiana

Datos

Abajo las instrucciones para cargar los datos del experimento.

datos <- read.table('datoslomb', header = TRUE)
head(datos)

##   Cama Comb Replica Humedad Temp Dieta Peso  pH Tiempo
## 1    1   11      II      80   15     1 5221 6.3     14
## 2    2   26      IV      80   25     2 5342 6.6     14
## 3    3   11     III      80   15     1 5221 6.3     14
## 4    4   18      IV      89   25     1 4903 6.7     14
## 5    5   24      II      89   20     2 5283 6.4     14
## 6    6    3     III      89   15     0 5110 6.6     14

Vamos a convertir a factor la variable Dieta.

datos$Dieta <- factor(datos$Dieta)

Vamos a quedarnos con las filas distintas de la base datos y la almacenamos en datis.

library(dplyr)
datis <- distinct(datos, Humedad, Temp, Dieta, Tiempo, Peso)

Distribuciones que mejor explican el Peso

library(gamlss)
mod <- fitDist(y=datis$Peso, type="realplus", k=log(nrow(datis)))

## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[1,1] = 0
## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[3,3] = 0
## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[4,4] = 0
## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[4,4] = 0

mod$fits

##           GG       IGAMMA           IG        LOGNO       LOGNO2 
##     721.0011     733.9338     734.3123     734.3525     734.3525 
##           GA          GB2         BCCG        BCCGo       exGAUS 
##     734.8076     735.6651     736.4313     736.4313     739.7765 
##         BCTo          BCT          WEI         WEI3          EXP 
##     740.3025     740.3025     741.4169     741.4169     926.3476 
##           GP      PARETO2     PARETO2o         WEI2 
##     930.2188     930.2188     930.2189 5200099.1739

mod$fails

## NULL

Explorando algunas relaciones entre Peso y las otras variables

par(mfrow=c(2, 2))
with(datis, plot(Peso ~ Tiempo, las=1, pch=20))
with(datis, plot(Peso ~ Humedad, las=1, pch=20))
with(datis, plot(Peso ~ Temp, las=1, pch=20))
with(datis, plot(Peso ~ Dieta, las=1, pch=20))

Aplicación de modelos gamlss

El código mostrado a continuación modifica algunos parámteros internos para hacer más exhaustiva la estimación de los parámetros.

library(gamlss)
control1 <- gamlss.control(n.cyc=15000, c.crit=0.00001, trace=FALSE)
control2 <- glim.control(cyc=15000, bf.cyc=3000, trace=FALSE)

Para la construcción de los modelos se utilizará un proceso de selección de variables hacia adelante y el horizonte de modelación para cada parámetro será el horizonte de abajo.

horizonte <- formula(~ Temp + Humedad + Tiempo + Dieta +
                       I(Temp^2) + I(Humedad^2) + I(Tiempo^2))

Modelo con respuesta normal modelando \(\mu\) solamente (de referencia)

El parámetro \(\mu\) es la media y \(\sigma\) es la desviación de la distribución.

control1 <- gamlss.control(n.cyc=15000, c.crit=0.001, trace=FALSE)
control2 <- glim.control(cyc=15000, bf.cyc=3000, trace=FALSE)
mod0 <- gamlss(Peso ~ 1,
               data=datis, control=control1, i.control=control2,
               family=NO)

mod0 <- stepGAICAll.A(object=mod0, trace=0,
                      scope=horizonte, 
                      sigma.scope=formula(~1))

## --------------------------------------------------- 
## Start:  AIC= 732.16 
##  Peso ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 699.56 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 699.56 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
## ---------------------------------------------------

summary(mod0)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = Peso ~ Tiempo + Dieta + I(Tiempo^2) +  
##     I(Temp^2), family = NO, data = datis, control = control1,  
##     i.control = control2, trace = FALSE, sigma.formula = ~1) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 2284.5252   781.6918   2.923  0.00563 **
## Tiempo       267.6698    79.2660   3.377  0.00162 **
## Dieta1      -291.2293   102.7449  -2.834  0.00709 **
## Dieta2       -75.4635   122.7581  -0.615  0.54213   
## I(Tiempo^2)   -5.2409     1.9293  -2.716  0.00962 **
## I(Temp^2)      0.6468     0.2785   2.322  0.02528 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.7223     0.1021   56.07   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  41 
##                       at cycle:  2 
##  
## Global Deviance:     685.5562 
##             AIC:     699.5562 
##             SBC:     712.6546 
## ******************************************************************

A continuación se construye el gráfico de gusano el cual es un herramienta para explorar el comportamiento de los residuales. Si el modelo está bien ajustado se espera que los puntos (residuales) queden dentro de las hipérbolas.

wp(mod0)

De la figura anterior se observa que los residuales caen en la parte central (dentro de las hipérbolas).

A continuación calculamos el coeficiente de correlación lineal entre los \(y_i\) y \(\hat{E}(y | x)\) para el modelo en cuestión.

y0 <- fitted(object=mod0, what='mu')
cor(datis$Peso, y0)

## [1] 0.7670557

El valor de \(R^2_{pseudo}\) es:

Rsq(mod0)

## [1] 0.5883744

Para ver si los valores \(\hat{E}(y | x)\) acompaña el valor de \(y_i\) podemos hacer el siguiente diagrama de dispersión.

plot(x=datis$Peso, y=y0, las=1)
abline(a=0, b=1, lty='dashed')

Modelo con respuesta normal modelando \(\mu\) y \(\sigma\)

El parámetro \(\mu\) es la media y \(\sigma\) es la desviación de la distribución.

mod1 <- gamlss(Peso ~ 1,
               sigma.fo= ~ 1, 
               data=datis, control=control1, i.control=control2,
               family=NO)

mod1 <- stepGAICAll.A(object=mod1, trace=0,
                      scope=horizonte, 
                      sigma.scope=horizonte)

## --------------------------------------------------- 
## Start:  AIC= 732.16 
##  Peso ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 699.56 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 596.71 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
## ---------------------------------------------------

summary(mod1)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = Peso ~ Tiempo + Dieta + I(Tiempo^2) +  
##     I(Temp^2), sigma.formula = ~I(Temp^2) + Temp +  
##     Dieta + I(Humedad^2) + Humedad, family = NO, data = datis,  
##     control = control1, i.control = control2, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) -2.674e+02  2.434e+01   -10.987 6.83e-13 ***
## Tiempo       4.724e+02  5.674e-02  8325.956  < 2e-16 ***
## Dieta1      -1.786e+02  4.047e+00   -44.133  < 2e-16 ***
## Dieta2       1.055e+02  3.488e+01     3.025  0.00464 ** 
## I(Tiempo^2) -9.388e+00  1.391e-03 -6751.324  < 2e-16 ***
## I(Temp^2)    9.203e-01  3.900e-02    23.601  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.431e+01  1.816e+00 -40.925  < 2e-16 ***
## I(Temp^2)    -7.331e-02  3.747e-03 -19.564  < 2e-16 ***
## Temp          2.293e+00  1.510e-01  15.184  < 2e-16 ***
## Dieta1       -1.651e+00  2.285e-01  -7.227 1.95e-08 ***
## Dieta2        2.734e+00  3.196e-01   8.556 4.25e-10 ***
## I(Humedad^2) -1.235e-02  1.823e-04 -67.760  < 2e-16 ***
## Humedad       1.795e+00  2.704e-02  66.375  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  13
##       Residual Deg. of Freedom:  35 
##                       at cycle:  16 
##  
## Global Deviance:     570.709 
##             AIC:     596.709 
##             SBC:     621.0346 
## ******************************************************************

mod1 <- update(mod1, what='mu',
               formula=Peso ~ Tiempo + I(Tiempo^2) + Dieta + Temp + I(Temp^2))

A continuación se construye el gráfico de gusano el cual es un herramienta para explorar el comportamiento de los residuales. Si el modelo está bien ajustado se espera que los puntos (residuales) queden dentro de las hipérbolas.

wp(mod1)

De la figura anterior se observa que los residuales caen en la parte central (dentro de las hipérbolas).

A continuación calculamos el coeficiente de correlación lineal entre los \(y_i\) y \(\hat{E}(y | x)\) para el modelo en cuestión.

y1 <- fitted(object=mod1, what='mu')
cor(datis$Peso, y1)

## [1] 0.609807

El valor de \(R^2_{pseudo}\) es:

Rsq(mod1)

## [1] 0.8963201

Para ver si los valores \(\hat{E}(y | x)\) acompaña el valor de \(y_i\) podemos hacer el siguiente diagrama de dispersión.

plot(x=datis$Peso, y=y1, las=1)
abline(a=0, b=1, lty='dashed')

Modelo con respuesta IGAMMA (Inverse Gamma)

El parámetro \(\mu\) es la moda de la distribución.

mod2 <- gamlss(Peso ~ 1,
               sigma.fo= ~ 1, 
               nu.fo= ~ 1,
               data=datis, control=control1, i.control=control2, 
               family=IGAMMA)

mod2 <- stepGAICAll.A(object=mod2,
                      scope = horizonte, 
                      sigma.scope = horizonte,
                      nu.scope = horizonte)

## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 730.19 
##  Peso ~ 1 
## 
##                Df    AIC
## + Tiempo        1 711.80
## + I(Tiempo^2)   1 714.40
## + Dieta         2 726.26
## + Temp          1 728.40
## + I(Temp^2)     1 728.46
## <none>            730.19
## + Humedad       1 730.57
## + I(Humedad^2)  1 730.88
## 
## Step:  AIC= 711.8 
##  Peso ~ Tiempo 
## 
##                Df    AIC
## + Dieta         2 704.67
## + Temp          1 706.80
## + I(Temp^2)     1 706.84
## + I(Tiempo^2)   1 708.08
## + Humedad       1 711.78
## <none>            711.80
## + I(Humedad^2)  1 712.19
## 
## Step:  AIC= 704.67 
##  Peso ~ Tiempo + Dieta 
## 
##                Df    AIC
## + I(Tiempo^2)   1 700.68
## + I(Temp^2)     1 702.56
## + Temp          1 702.69
## <none>            704.67
## + Humedad       1 706.10
## + I(Humedad^2)  1 706.32
## 
## Step:  AIC= 700.68 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) 
## 
##                Df    AIC
## + I(Temp^2)     1 697.71
## + Temp          1 697.82
## <none>            700.68
## + Humedad       1 702.12
## + I(Humedad^2)  1 702.34
## 
## Step:  AIC= 697.71 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
##                Df    AIC
## <none>            697.71
## + Humedad       1 698.68
## + I(Humedad^2)  1 699.01
## + Temp          1 699.64
## --------------------------------------------------- 
## Distribution parameter:  sigma 
## Start:  AIC= 697.71 
##  ~1 
## 
##                Df    AIC
## + I(Temp^2)     1 681.16
## + Temp          1 685.44
## + Tiempo        1 694.95
## + I(Tiempo^2)   1 695.97
## + Dieta         2 696.93
## <none>            697.71
## + I(Humedad^2)  1 698.62
## + Humedad       1 698.77
## 
## Step:  AIC= 681.16 
##  ~I(Temp^2) 
## 
##                Df    AIC
## + Temp          1 667.83
## <none>            681.16
## + Tiempo        1 681.88
## + I(Tiempo^2)   1 682.48
## + Dieta         2 683.12
## + Humedad       1 683.14
## + I(Humedad^2)  1 683.15
## 
## Step:  AIC= 667.83 
##  ~I(Temp^2) + Temp 
## 
##                Df    AIC
## + Dieta         2 652.93
## <none>            667.83
## + Tiempo        1 669.35
## + I(Tiempo^2)   1 669.61
## + I(Humedad^2)  1 669.77
## + Humedad       1 669.79
## 
## Step:  AIC= 652.93 
##  ~I(Temp^2) + Temp + Dieta 
## 
##                Df    AIC
## + I(Humedad^2)  1 613.88
## + Humedad       1 615.76
## <none>            652.93
## + I(Tiempo^2)   1 654.79
## + Tiempo        1 654.85
## 
## Step:  AIC= 613.88 
##  ~I(Temp^2) + Temp + Dieta + I(Humedad^2) 
## 
##               Df    AIC
## <none>           613.88
## + Tiempo       1 615.33
## + I(Tiempo^2)  1 615.63
## + Humedad      1 680.30
## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 613.88 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
##               Df    AIC
## <none>           613.88
## - I(Temp^2)    1 614.42
## - Dieta        2 692.94
## - I(Tiempo^2)  1 707.76
## - Tiempo       1 712.05
## ---------------------------------------------------

summary(mod2)

## Warning in summary.gamlss(mod2): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("IGAMMA", "Inverse Gamma") 
## 
## Call:  gamlss(formula = Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2),  
##     sigma.formula = ~I(Temp^2) + Temp + Dieta + I(Humedad^2),  
##     nu.formula = ~1, family = IGAMMA, data = datis, control = control1,  
##     i.control = control2, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  7.596e+00  2.646e-02   287.103   <2e-16 ***
## Tiempo       8.818e-02  1.202e-05  7337.215   <2e-16 ***
## Dieta1      -3.262e-02  1.567e-03   -20.823   <2e-16 ***
## Dieta2       1.798e-02  7.408e-03     2.427   0.0196 *  
## I(Tiempo^2) -1.757e-03  2.944e-07 -5967.233   <2e-16 ***
## I(Temp^2)    7.116e-05  4.239e-05     1.679   0.1006    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.923e+00  3.638e+00  -1.628    0.111    
## I(Temp^2)    -5.949e-02  9.117e-03  -6.525 6.96e-08 ***
## Temp          1.740e+00  3.687e-01   4.718 2.64e-05 ***
## Dieta1       -1.830e+00  2.392e-01  -7.650 1.73e-09 ***
## Dieta2        3.350e+00  2.901e-01  11.548 1.29e-14 ***
## I(Humedad^2) -1.213e-03  8.537e-05 -14.212  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  12
##       Residual Deg. of Freedom:  36 
##                       at cycle:  26 
##  
## Global Deviance:     589.8823 
##             AIC:     613.8823 
##             SBC:     636.3367 
## ******************************************************************

De la tabla de resumen anterior se observa que \(Temp^2\) no es significativo en la modelación de \(\mu\), por esta razón se va a ajustar nuevamente el modelo si él.

mod2 <- update(mod2, what='mu',
               formula=Peso ~ Tiempo + Dieta + I(Tiempo^2))

A continuación se construye el gráfico de gusano el cual es un herramienta para explorar el comportamiento de los residuales. Si el modelo está bien ajustado se espera que los puntos (residuales) queden dentro de las hipérbolas.

wp(mod2)

De la figura anterior se observa que los residuales caen en la parte central (dentro de las hipérbolas).

A continuación calculamos el coeficiente de correlación lineal entre los \(y_i\) y \(\hat{E}(y | x)\) para el modelo en cuestión.

y2 <- fitted(object=mod2, what='mu')
cor(datis$Peso, y2)

## [1] 0.7078815

El valor de \(R^2_{pseudo}\) es:

Rsq(mod2)

## [1] 0.9383854

Para ver si los valores \(\hat{E}(y | x)\) acompaña el valor de \(y_i\) podemos hacer el siguiente diagrama de dispersión.

plot(x=datis$Peso, y=y2, las=1)
abline(a=0, b=1, lty='dashed')

Modelo con respuesta IG (Inverse Gaussian)

El parámetro \(\mu\) es la media de la distribución.

mod3 <- gamlss(Peso ~ 1,
               sigma.fo= ~ 1, 
               data=datis, control=control1, i.control=control2, 
               family=IG)

mod3 <- stepGAICAll.A(object=mod3,
                      scope = horizonte, 
                      sigma.scope = horizonte,
                      nu.scope = horizonte)

## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 730.57 
##  Peso ~ 1 
## 
##                Df    AIC
## + Tiempo        1 710.29
## + I(Tiempo^2)   1 713.21
## + Dieta         2 726.93
## + Temp          1 728.25
## + I(Temp^2)     1 728.32
## <none>            730.57
## + Humedad       1 731.19
## + I(Humedad^2)  1 731.50
## 
## Step:  AIC= 710.29 
##  Peso ~ Tiempo 
## 
##                Df    AIC
## + Dieta         2 703.78
## + Temp          1 705.95
## + I(Temp^2)     1 706.04
## + I(Tiempo^2)   1 706.13
## <none>            710.29
## + Humedad       1 710.72
## + I(Humedad^2)  1 711.10
## 
## Step:  AIC= 703.78 
##  Peso ~ Tiempo + Dieta 
## 
##                Df    AIC
## + I(Tiempo^2)   1 699.32
## + I(Temp^2)     1 702.44
## + Temp          1 702.56
## <none>            703.78
## + Humedad       1 705.00
## + I(Humedad^2)  1 705.27
## 
## Step:  AIC= 699.32 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) 
## 
##                Df    AIC
## + I(Temp^2)     1 697.63
## + Temp          1 697.72
## <none>            699.32
## + Humedad       1 700.55
## + I(Humedad^2)  1 700.82
## 
## Step:  AIC= 697.63 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
##                Df    AIC
## <none>            697.63
## + Humedad       1 698.57
## + I(Humedad^2)  1 698.92
## + Temp          1 699.56
## --------------------------------------------------- 
## Distribution parameter:  sigma 
## Start:  AIC= 697.63 
##  ~1 
## 
##                Df    AIC
## + I(Temp^2)     1 679.99
## + Temp          1 684.28
## + Tiempo        1 695.12
## + I(Tiempo^2)   1 696.04
## + Dieta         2 696.56
## <none>            697.63
## + I(Humedad^2)  1 698.30
## + Humedad       1 698.38
## 
## Step:  AIC= 679.99 
##  ~I(Temp^2) 
## 
##                Df    AIC
## + Temp          1 667.16
## + Dieta         2 674.30
## + Tiempo        1 677.43
## + I(Tiempo^2)   1 679.18
## <none>            679.99
## + I(Humedad^2)  1 681.99
## + Humedad       1 681.99
## 
## Step:  AIC= 667.16 
##  ~I(Temp^2) + Temp 
## 
##                Df    AIC
## + Dieta         2 649.52
## + Tiempo        1 666.91
## <none>            667.16
## + I(Tiempo^2)   1 668.45
## + I(Humedad^2)  1 669.08
## + Humedad       1 669.10
## 
## Step:  AIC= 649.52 
##  ~I(Temp^2) + Temp + Dieta 
## 
##                Df    AIC
## + I(Humedad^2)  1 610.16
## + Humedad       1 611.24
## + Tiempo        1 646.20
## + I(Tiempo^2)   1 648.03
## <none>            649.52
## 
## Step:  AIC= 610.16 
##  ~I(Temp^2) + Temp + Dieta + I(Humedad^2) 
## 
##               Df    AIC
## <none>           610.16
## + Tiempo       1 611.75
## + I(Tiempo^2)  1 611.87
## + Humedad      1 650.89
## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 610.16 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
##               Df    AIC
## - I(Temp^2)    1 608.36
## <none>           610.16
## - Dieta        2 692.03
## - I(Tiempo^2)  1 706.94
## - Tiempo       1 709.64
## 
## Step:  AIC= 608.36 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) 
## 
##               Df    AIC
## <none>           608.36
## - Dieta        2 692.71
## - I(Tiempo^2)  1 710.10
## - Tiempo       1 715.26
## ---------------------------------------------------

summary(mod3)

## Warning in summary.gamlss(mod3): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("IG", "Inverse Gaussian") 
## 
## Call:  
## gamlss(formula = Peso ~ Tiempo + Dieta + I(Tiempo^2), sigma.formula = ~I(Temp^2) +  
##     Temp + Dieta + I(Humedad^2), family = IG, data = datis, control = control1,  
##     i.control = control2, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  7.640e+00  1.512e-03  5053.882  < 2e-16 ***
## Tiempo       8.818e-02  1.241e-05  7106.686  < 2e-16 ***
## Dieta1      -3.284e-02  1.507e-03   -21.791  < 2e-16 ***
## Dieta2       2.120e-02  7.780e-03     2.725  0.00927 ** 
## I(Tiempo^2) -1.757e-03  3.053e-07 -5754.272  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.085e+01  3.732e+00  -2.907  0.00581 ** 
## I(Temp^2)    -6.083e-02  9.300e-03  -6.541 6.63e-08 ***
## Temp          1.791e+00  3.763e-01   4.760 2.31e-05 ***
## Dieta1       -1.857e+00  2.411e-01  -7.700 1.47e-09 ***
## Dieta2        3.374e+00  2.965e-01  11.379 2.07e-14 ***
## I(Humedad^2) -1.188e-03  8.742e-05 -13.589  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  11
##       Residual Deg. of Freedom:  37 
##                       at cycle:  18 
##  
## Global Deviance:     586.3558 
##             AIC:     608.3558 
##             SBC:     628.939 
## ******************************************************************

A continuación se construye el gráfico de gusano el cual es un herramienta para explorar el comportamiento de los residuales. Si el modelo está bien ajustado se espera que los puntos (residuales) queden dentro de las hipérbolas.

wp(mod3)

De la figura anterior se observa que los residuales caen en la parte central (dentro de las hipérbolas).

A continuación calculamos el coeficiente de correlación lineal entre los \(y_i\) y \(\hat{E}(y | x)\) para el modelo en cuestión.

y3 <- fitted(object=mod3, what='mu')
cor(datis$Peso, y3)

## [1] 0.708554

El valor de \(R^2_{pseudo}\) es:

Rsq(mod3)

## [1] 0.9461271

Para ver si los valores \(\hat{E}(y | x)\) acompaña el valor de \(y_i\) podemos hacer el siguiente diagrama de dispersión.

plot(x=datis$Peso, y=y3, las=1)
abline(a=0, b=1, lty='dashed')

Modelo con respuesta GA (Gamma)

El parámetro \(\mu\) es la media de la distribución.

mod4 <- gamlss(Peso ~ 1,
               sigma.fo= ~ 1, 
               data=datis, control=control1, i.control=control2, 
               family=GA)

mod4 <- stepGAICAll.A(object=mod4,
                      scope = horizonte, 
                      sigma.scope = horizonte,
                      nu.scope = horizonte)

## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 731.07 
##  Peso ~ 1 
## 
##                Df    AIC
## + Tiempo        1 711.31
## + I(Tiempo^2)   1 714.18
## + Dieta         2 727.49
## + Temp          1 728.77
## + I(Temp^2)     1 728.85
## <none>            731.07
## + Humedad       1 731.74
## + I(Humedad^2)  1 732.04
## 
## Step:  AIC= 711.31 
##  Peso ~ Tiempo 
## 
##                Df    AIC
## + Dieta         2 704.80
## + Temp          1 706.48
## + I(Temp^2)     1 706.56
## + I(Tiempo^2)   1 707.04
## <none>            711.31
## + Humedad       1 711.77
## + I(Humedad^2)  1 712.14
## 
## Step:  AIC= 704.8 
##  Peso ~ Tiempo + Dieta 
## 
##                Df    AIC
## + I(Tiempo^2)   1 700.23
## + I(Temp^2)     1 702.91
## + Temp          1 703.03
## <none>            704.80
## + Humedad       1 706.20
## + I(Humedad^2)  1 706.42
## 
## Step:  AIC= 700.23 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) 
## 
##                Df    AIC
## + I(Temp^2)     1 697.78
## + Temp          1 697.87
## <none>            700.23
## + Humedad       1 701.64
## + I(Humedad^2)  1 701.87
## 
## Step:  AIC= 697.78 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
##                Df    AIC
## <none>            697.78
## + Humedad       1 698.78
## + I(Humedad^2)  1 699.12
## + Temp          1 699.72
## --------------------------------------------------- 
## Distribution parameter:  sigma 
## Start:  AIC= 697.78 
##  ~1 
## 
##                Df    AIC
## + I(Temp^2)     1 680.55
## + Temp          1 684.81
## + Tiempo        1 695.02
## + I(Tiempo^2)   1 695.99
## + Dieta         2 697.01
## <none>            697.78
## + I(Humedad^2)  1 698.55
## + Humedad       1 698.65
## 
## Step:  AIC= 680.55 
##  ~I(Temp^2) 
## 
##                Df    AIC
## + Temp          1 667.29
## + Dieta         2 676.54
## + Tiempo        1 678.68
## + I(Tiempo^2)   1 680.23
## <none>            680.55
## + Humedad       1 682.53
## + I(Humedad^2)  1 682.53
## 
## Step:  AIC= 667.29 
##  ~I(Temp^2) + Temp 
## 
##                Df    AIC
## + Dieta         2 650.38
## <none>            667.29
## + Tiempo        1 668.30
## + I(Tiempo^2)   1 669.12
## + I(Humedad^2)  1 669.22
## + Humedad       1 669.24
## 
## Step:  AIC= 650.38 
##  ~I(Temp^2) + Temp + Dieta 
## 
##                Df    AIC
## + I(Humedad^2)  1 610.70
## + Humedad       1 611.92
## + Tiempo        1 648.91
## + I(Tiempo^2)   1 650.16
## <none>            650.38
## 
## Step:  AIC= 610.7 
##  ~I(Temp^2) + Temp + Dieta + I(Humedad^2) 
## 
##               Df    AIC
## <none>           610.70
## + I(Tiempo^2)  1 612.64
## + Tiempo       1 612.68
## + Humedad      1 652.47
## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 610.7 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) + I(Temp^2) 
## 
##               Df    AIC
## - I(Temp^2)    1 609.26
## <none>           610.70
## - Dieta        2 691.28
## - I(Tiempo^2)  1 707.86
## - Tiempo       1 711.05
## 
## Step:  AIC= 609.26 
##  Peso ~ Tiempo + Dieta + I(Tiempo^2) 
## 
##               Df    AIC
## <none>           609.26
## - Dieta        2 640.48
## - I(Tiempo^2)  1 711.88
## - Tiempo       1 716.62
## ---------------------------------------------------

summary(mod4)

## Warning in summary.gamlss(mod4): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  
## gamlss(formula = Peso ~ Tiempo + Dieta + I(Tiempo^2), sigma.formula = ~I(Temp^2) +  
##     Temp + Dieta + I(Humedad^2), family = GA, data = datis, control = control1,  
##     i.control = control2, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  7.640e+00  1.549e-03  4932.799   <2e-16 ***
## Tiempo       8.818e-02  1.160e-05  7598.815   <2e-16 ***
## Dieta1      -3.258e-02  1.545e-03   -21.088   <2e-16 ***
## Dieta2       1.921e-02  7.570e-03     2.538   0.0149 *  
## I(Tiempo^2) -1.757e-03  2.843e-07 -6179.988   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.585e+00  3.713e+00  -1.773   0.0835 .  
## I(Temp^2)    -6.118e-02  9.263e-03  -6.605 5.35e-08 ***
## Temp          1.804e+00  3.747e-01   4.814 1.94e-05 ***
## Dieta1       -1.906e+00  2.408e-01  -7.917 7.31e-10 ***
## Dieta2        3.358e+00  2.952e-01  11.375 2.09e-14 ***
## I(Humedad^2) -1.195e-03  8.704e-05 -13.729  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  11
##       Residual Deg. of Freedom:  37 
##                       at cycle:  17 
##  
## Global Deviance:     587.2553 
##             AIC:     609.2553 
##             SBC:     629.8385 
## ******************************************************************

A continuación se construye el gráfico de gusano el cual es un herramienta para explorar el comportamiento de los residuales. Si el modelo está bien ajustado se espera que los puntos (residuales) queden dentro de las hipérbolas.

wp(mod4)

De la figura anterior se observa que los residuales caen en la parte central (dentro de las hipérbolas).

A continuación calculamos el coeficiente de correlación lineal entre los \(y_i\) y \(\hat{E}(y | x)\) para el modelo en cuestión.

y4 <- fitted(object=mod4, what='mu')
cor(datis$Peso, y4)

## [1] 0.7081088

El valor de \(R^2_{pseudo}\) es:

Rsq(mod4)

## [1] 0.9456714

Para ver si los valores \(\hat{E}(y | x)\) acompaña el valor de \(y_i\) podemos hacer el siguiente diagrama de dispersión.

plot(x=datis$Peso, y=y4, las=1)
abline(a=0, b=1, lty='dashed')

Modelo con respuesta GG (Generalized Gamma)

Para este modelo se tuvo que cambiar el párametro de control control1.

control1 <- gamlss.control(n.cyc=15000, c.crit=0.01, trace=FALSE)
control2 <- glim.control(cyc=5000, bf.cyc=3000, trace=FALSE)

mod5 <- gamlss(Peso ~ 1,
               sigma.fo= ~ 1, 
               nu.fo= ~ 1,
               data=datis, control=control1, i.control=control2,
               family=GG)

mod5 <- stepGAICAll.A(object=mod5, trace=0,
                      scope = horizonte, 
                      sigma.scope = horizonte,
                      nu.scope = horizonte)

## --------------------------------------------------- 
## Start:  AIC= 715.69 
##  Peso ~ 1 
## 
## Model with term  Dieta has failed 
## Model with term  Dieta has failed 
## --------------------------------------------------- 
## Start:  AIC= 698.42 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 611.22 
##  ~1 
## 
## Model with term  Temp has failed 
## Model with term  Humedad has failed 
## Model with term  I(Temp^2) has failed 
## Model with term  I(Humedad^2) has failed 
## --------------------------------------------------- 
## Start:  AIC= 595.21 
##  ~I(Temp^2) + Temp + Dieta + I(Humedad^2) 
## 
## Model with term  I(Humedad^2) has failed 
## --------------------------------------------------- 
## Start:  AIC= 595.21 
##  Peso ~ Tiempo + I(Tiempo^2) + Dieta 
## 
## Model with term  I(Tiempo^2) has failed 
## ---------------------------------------------------

summary(mod5)

## Warning in summary.gamlss(mod5): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("GG", "generalised Gamma Lopatatsidis-Green") 
## 
## Call:  
## gamlss(formula = Peso ~ Tiempo + I(Tiempo^2) + Dieta, sigma.formula = ~I(Temp^2) +  
##     Temp + Dieta + I(Humedad^2), nu.formula = ~Tiempo + I(Tiempo^2) +  
##     Dieta, family = GG, data = datis, control = control1, i.control = control2,  
##     trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  7.639e+00  9.147e-04  8351.824   <2e-16 ***
## Tiempo       8.818e-02  1.346e-05  6553.090   <2e-16 ***
## I(Tiempo^2) -1.757e-03  3.297e-07 -5329.520   <2e-16 ***
## Dieta1      -3.157e-02  9.056e-04   -34.863   <2e-16 ***
## Dieta2       3.417e-03  1.845e-03     1.853   0.0708 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.117e+01  2.835e+00  -3.940 0.000302 ***
## I(Temp^2)    -5.595e-02  7.059e-03  -7.927 7.09e-10 ***
## Temp          1.752e+00  2.840e-01   6.169 2.27e-07 ***
## Dieta1       -2.512e+00  1.985e-01 -12.654 6.40e-16 ***
## Dieta2        1.721e+00  2.403e-01   7.159 8.61e-09 ***
## I(Humedad^2) -7.681e-04  7.451e-05 -10.308 4.50e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7008.715   1269.699  -5.520 1.82e-06 ***
## Tiempo        593.838    106.784   5.561 1.58e-06 ***
## I(Tiempo^2)   -12.211      2.203  -5.544 1.68e-06 ***
## Dieta1        598.584    187.876   3.186  0.00269 ** 
## Dieta2        -19.762     36.702  -0.538  0.59304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  16
##       Residual Deg. of Freedom:  32 
##                       at cycle:  60 
##  
## Global Deviance:     563.2089 
##             AIC:     595.2089 
##             SBC:     625.1481 
## ******************************************************************

A continuación se construye el gráfico de gusano el cual es un herramienta para explorar el comportamiento de los residuales. Si el modelo está bien ajustado se espera que los puntos (residuales) queden dentro de las hipérbolas.

wp(mod5)

De la figura anterior se observa que los residuales caen en la parte central (dentro de las hipérbolas).

A continuación calculamos el coeficiente de correlación lineal entre los \(y_i\) y \(\hat{E}(y | x)\) para el modelo en cuestión.

y5 <- fitted(object=mod5, what='mu')
cor(datis$Peso, y5)

## [1] 0.7035563

El valor de \(R^2_{pseudo}\) es:

Rsq(mod5)

## [1] 0.9524221

Para ver si los valores \(\hat{E}(y | x)\) acompaña el valor de \(y_i\) podemos hacer el siguiente diagrama de dispersión.

plot(x=datis$Peso, y=y5, las=1)
abline(a=0, b=1, lty='dashed')

Comparación los modelos

En esta sección se muestran los resultados comparando los modelos.

Usando \(\rho(y, \hat{E}(y | x))\) y \(R^2_{pseudo}\)

correl <- c(cor(datis$Peso, y0), cor(datis$Peso, y1), cor(datis$Peso, y2),
            cor(datis$Peso, y3), cor(datis$Peso, y4), cor(datis$Peso, y5))

r2_pseudo <- c(Rsq(mod0), Rsq(mod1), Rsq(mod2), 
               Rsq(mod3), Rsq(mod4), Rsq(mod5))

tabla <- cbind(correlacion=correl, r2_pseudo=r2_pseudo)
tabla <- round(tabla, 2)
rownames(tabla) <- paste0('mod', 0:5)

library(kableExtra)

## Warning: package 'kableExtra' was built under R version 3.6.3

kable(tabla) %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)

	correlacion	r2_pseudo
mod0	0.77	0.59
mod1	0.61	0.90
mod2	0.71	0.94
mod3	0.71	0.95
mod4	0.71	0.95
mod5	0.70	0.95

De la tabla anterior se observa que mod1 tiene el segundo mejor \(\rho\) y el mejor \(R^2_{pseudo}\).

Usando el criterio del GAIC

Otra forma de comparar modelos no anidados es por medio del GAIC.

aic <- AIC(mod0, mod1, mod2, mod3, mod4, mod5, k=log(nrow(datis)))

kable(aic) %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

	df	AIC
mod5	16	625.1481
mod3	11	628.9390
mod4	11	629.8385
mod2	11	635.0056
mod1	14	673.5703
mod0	7	712.6546

De la tabla anterior se observa que el mejor modelo es mod1 ya que presenta el menor GAIC.

Residuales

Los residuales de los modelos considerados se muestran a continuación. Si el modelo fue bien ajustado, los residuales deben tener una distribución \(N(0, 1)\).

library(car)
par(mfrow=c(3, 2))
qqPlot(resid(mod0), dist="norm", mean=0, sd=1, pch=20, id=FALSE)
qqPlot(resid(mod1), dist="norm", mean=0, sd=1, pch=20, id=FALSE)
qqPlot(resid(mod2), dist="norm", mean=0, sd=1, pch=20, id=FALSE)
qqPlot(resid(mod3), dist="norm", mean=0, sd=1, pch=20, id=FALSE)
qqPlot(resid(mod4), dist="norm", mean=0, sd=1, pch=20, id=FALSE)
qqPlot(resid(mod5), dist="norm", mean=0, sd=1, pch=20, id=FALSE)

Los residuales del modelo mod1 no violan gravemente el supuesto de normalidad.

El mejor modelo

El mejor modelo identificado hasta ahora fue el modelo mod1 que tiene respuesta normal y para el cual se modeló \(\mu\) y \(\sigma\).

A continuación la tabla de resumen del modelo.

summary(mod1)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = Peso ~ Tiempo + I(Tiempo^2) + Dieta +  
##     Temp + I(Temp^2), sigma.formula = ~I(Temp^2) +  
##     Temp + Dieta + I(Humedad^2) + Humedad, family = NO,  
##     data = datis, control = control1, i.control = control2,  
##     trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) 9098.00719   98.36123    92.496  < 2e-16 ***
## Tiempo       357.89672    0.12770  2802.583  < 2e-16 ***
## I(Tiempo^2)   -7.05358    0.00311 -2268.099  < 2e-16 ***
## Dieta1      -183.14956   58.33909    -3.139 0.003492 ** 
## Dieta2       139.11145   34.22307     4.065 0.000269 ***
## Temp        -771.08472    1.89764  -406.340  < 2e-16 ***
## I(Temp^2)     18.80639    0.19103    98.446  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)   1.329e+02  1.341e+00   99.093   <2e-16 ***
## I(Temp^2)    -1.221e-01  2.888e-03  -42.276   <2e-16 ***
## Temp          4.625e+00  1.174e-01   39.408   <2e-16 ***
## Dieta1       -2.344e-01  2.982e-01   -0.786    0.437    
## Dieta2       -4.002e+00  2.142e-01  -18.681   <2e-16 ***
## I(Humedad^2)  2.729e-02  1.406e-04  194.132   <2e-16 ***
## Humedad      -4.307e+00  1.786e-02 -241.149   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  48 
## Degrees of Freedom for the fit:  14
##       Residual Deg. of Freedom:  34 
##                       at cycle:  10 
##  
## Global Deviance:     619.3734 
##             AIC:     647.3734 
##             SBC:     673.5703 
## ******************************************************************

La media estimada para el peso está dada por:

\[ \widehat{E}(y) = 9098.01 + 357.9 \times Time + -7.05 \times Time^2 + -183.15 \times Dieta1 + 139.11 \times Dieta2 + -771.08 \times Temp + 18.81 \times Temp^2 \]

La desviación estimada para el peso está dada por:

\[ \log(\widehat{Des}(y)) = 132.92 + 4.62 \times Temp + -0.12 \times Temp^2 + -0.23 \times Dieta1 + -4 \times Dieta2 + -4.31 \times Hum + 0.03 \times Hum^2 \]

Maximizando

A continuación se muestra tres diagramas de contornos para \(\widehat{E}(y | tiempo, temp)\) para cada una de las dietas. Los puntos en rojo representan la combinación de tiempo y temperatura que maximizan el peso.

A continuación las tres superficies.

Los valores de tiempo y temperatura que maximizan el peso son los siguientes:

	Tiempo	Temperatura	Peso máximo
Dieta 0	25.36987	15	6303.072
Dieta 1	25.36987	15	6119.923
Dieta 2	25.36987	15	6442.184