Ejemplos Modelos Lineales: análisis de la Varianza y la Covarianza
Angelo Santana
9 de noviembre de 2016
ejemplo2=read.csv2("datos/ejemplo2.csv")
ejemplo3=read.csv2("datos/ejemplo3.csv")
ejemplo4=read.csv2("datos/ejemplo4.csv")
ejemplo5=read.csv2("datos/ejemplo5.csv")
ejemplo6=read.csv2("datos/ejemplo6.csv")Primera exploración: datos y representaciones gráficas.
A continuación vemos el contenido de estos ficheros y representamos los datos.
Ejemplo 2
Una compañía algodonera, interesada en maximizar el rendimiento de la semilla de algodón, quiere evaluar si dicho rendimiento depende del tipo de fertilizante utilizado. Para ello dispone de cinco tipos de fertilizantes y para comparar su eficacia fumiga con cada uno un cierto número de parcelas de terreno de la misma calidad e igual superficie. Los datos resultantes se muestran a continuación.
ejemplo2## fertilizante rendimiento
## 1 F1 51
## 2 F1 49
## 3 F1 50
## 4 F1 49
## 5 F1 51
## 6 F1 50
## 7 F2 56
## 8 F2 60
## 9 F2 56
## 10 F2 56
## 11 F2 57
## 12 F3 48
## 13 F3 50
## 14 F3 53
## 15 F3 44
## 16 F3 45
## 17 F4 47
## 18 F4 48
## 19 F4 49
## 20 F4 44
## 21 F5 43
## 22 F5 43
## 23 F5 46
## 24 F5 47
## 25 F5 45
## 26 F5 46par(mfrow=c(1,2))
plot(ejemplo2)
plot(as.numeric(ejemplo2$fertilizante),ejemplo2$rendimiento,xlab="Fertilizante")
Ejemplo 3
Una industria química desea probar el efecto que tienen cuatro agentes químicos sobre la resistencia de un tipo particular de tela. Como puede haber variabilidad entre los rollos de tela, decide utilizar un diseño aleatorizado por bloques considerando los rollos de tela como bloques. Se seleccionan entonces 5 rollos de tela y a cada uno se le aplica los cuatro agentes químicos en orden aleatorio. Seguidamente se muestran los resultados de la resistencia a la tensión según tratamiento (agente químico) y rollo de tela (bloque).
ejemplo3## agenteQuim rolloTela resistencia
## 1 AQ1 RT1 73
## 2 AQ1 RT2 68
## 3 AQ1 RT3 74
## 4 AQ1 RT4 71
## 5 AQ1 RT5 67
## 6 AQ2 RT1 73
## 7 AQ2 RT2 67
## 8 AQ2 RT3 75
## 9 AQ2 RT4 72
## 10 AQ2 RT5 70
## 11 AQ3 RT1 75
## 12 AQ3 RT2 68
## 13 AQ3 RT3 78
## 14 AQ3 RT4 73
## 15 AQ3 RT5 68
## 16 AQ4 RT1 73
## 17 AQ4 RT2 71
## 18 AQ4 RT3 75
## 19 AQ4 RT4 75
## 20 AQ4 RT5 69par(mfrow=c(1,2))
boxplot(resistencia~agenteQuim,data=ejemplo3,ylab="Resistencia",xlab="Agente Químico")
boxplot(resistencia~rolloTela,data=ejemplo3,ylab="Resistencia",xlab="Rollo de Tela")
Ejemplo 4
Se investigan los efectos sobre la resistencia del papel que producen el porcentaje de la concentración de fibra de madera en la pulpa, la presión del tanque y el tiempo de cocción de la pulpa. Para ello se seleccionan tres niveles de la concentración de madera, tres de la presión y dos del tiempo de cocción. En este experimento factorial se realizaron dos réplicas por cada condición experimental. Los valores obtenidos se muestran a continuación.
ejemplo4## resistencia concFibra tiempo presion
## 1 196.6 2% 3 h. 400
## 2 197.7 2% 3 h. 500
## 3 199.8 2% 3 h. 650
## 4 198.4 2% 4 h. 400
## 5 199.6 2% 4 h. 500
## 6 200.6 2% 4 h. 650
## 7 196.0 2% 3 h. 400
## 8 196.0 2% 3 h. 500
## 9 199.4 2% 3 h. 650
## 10 198.6 2% 4 h. 400
## 11 200.4 2% 4 h. 500
## 12 200.9 2% 4 h. 650
## 13 198.5 4% 3 h. 400
## 14 196.0 4% 3 h. 500
## 15 198.4 4% 3 h. 650
## 16 197.5 4% 4 h. 400
## 17 198.7 4% 4 h. 500
## 18 199.6 4% 4 h. 650
## 19 197.2 4% 3 h. 400
## 20 196.9 4% 3 h. 500
## 21 197.6 4% 3 h. 650
## 22 198.1 4% 4 h. 400
## 23 198.0 4% 4 h. 500
## 24 199.0 4% 4 h. 650
## 25 197.5 8% 3 h. 400
## 26 195.6 8% 3 h. 500
## 27 197.4 8% 3 h. 650
## 28 197.6 8% 4 h. 400
## 29 197.0 8% 4 h. 500
## 30 198.5 8% 4 h. 650
## 31 196.6 8% 3 h. 400
## 32 196.2 8% 3 h. 500
## 33 198.1 8% 3 h. 650
## 34 198.4 8% 4 h. 400
## 35 197.8 8% 4 h. 500
## 36 199.8 8% 4 h. 650par(mfrow=c(1,3))
boxplot(resistencia~concFibra,data=ejemplo4,ylab="Resistencia",xlab="Concentración Fibra")
boxplot(resistencia~tiempo,data=ejemplo4,ylab="Resistencia",xlab="Tiempo")
boxplot(resistencia~presion,data=ejemplo4,ylab="Resistencia",xlab="Presion")
par(mfrow=c(1,3),las=2,cex.lab=1.5,cex.axis=1.2)
boxplot(resistencia~tiempo+presion,data=ejemplo4,ylab="Resistencia")
boxplot(resistencia~tiempo+concFibra,data=ejemplo4,ylab="Resistencia")
boxplot(resistencia~concFibra+presion,data=ejemplo4,ylab="Resistencia")
par(mfrow=c(2,3),las=2,cex.lab=1.5,cex.axis=1.2)
attach(ejemplo4)
interaction.plot(concFibra,tiempo,resistencia)
interaction.plot(tiempo,concFibra,resistencia)
interaction.plot(tiempo,presion,resistencia)
interaction.plot(presion,tiempo,resistencia)
interaction.plot(concFibra,presion,resistencia)
interaction.plot(presion,concFibra,resistencia)
detach(ejemplo4)Ejemplo 5
Se realizó un estudio experimental para evaluar el acabado superficial de partes metálicas producidas por cuatro máquinas. Cada máquina es operada por tres diferentes operarios, observándose dos piezas producidas por cada uno. Debido a que las máquinas están en localidades diferentes no es posible usar los mismos operarios en cada máquina. De esta forma los efectos de los operarios están anidados en la máquina que utilizan. Las respuestas observadas para cada una de las condiciones experimentales se muestran a continuación.
ejemplo5## acabado maquina operario
## 1 79 M1 op1
## 2 94 M1 op2
## 3 46 M1 op3
## 4 92 M2 op1
## 5 85 M2 op2
## 6 76 M2 op3
## 7 88 M3 op1
## 8 53 M3 op2
## 9 46 M3 op3
## 10 36 M4 op1
## 11 40 M4 op2
## 12 62 M4 op3
## 13 62 M1 op1
## 14 74 M1 op2
## 15 57 M1 op3
## 16 99 M2 op1
## 17 79 M2 op2
## 18 68 M2 op3
## 19 75 M3 op1
## 20 56 M3 op2
## 21 57 M3 op3
## 22 53 M4 op1
## 23 56 M4 op2
## 24 47 M4 op3par(mfrow=c(1,2))
boxplot(acabado~maquina,data=ejemplo5,ylab="acabado",xlab="maquina")
boxplot(acabado~operario,data=ejemplo5,ylab="acabado",xlab="operario")
par(mfrow=c(1,2))
attach(ejemplo5)
interaction.plot(maquina,operario,acabado)
interaction.plot(operario,maquina,acabado)
detach(ejemplo5)Ejemplo 6
Se estudia el voltaje generado en una dinamo en función de la velocidad de giro del rotor y del tipo de diseño del generador. Hay cuatro tipos de diseño distintos y cada uno ha sido observado entre 4 y seis veces. La velocidad de giro del motor se mide en rpm. Los datos se muestran a continuación.
ejemplo6## voltaje velocidad tipo
## 1 1.95 6336 T1
## 2 2.50 7099 T1
## 3 2.93 8026 T1
## 4 1.69 6230 T1
## 5 1.23 5369 T1
## 6 3.13 8343 T1
## 7 1.55 6522 T2
## 8 1.94 7310 T2
## 9 2.18 7974 T2
## 10 2.70 8501 T2
## 11 1.32 6646 T3
## 12 1.60 7384 T3
## 13 1.89 8000 T3
## 14 2.15 8545 T3
## 15 1.09 6755 T4
## 16 1.26 7362 T4
## 17 1.57 7934 T4
## 18 1.92 8554 T4require(car)
scatterplot(voltaje~velocidad|tipo,reg.line=lm,smooth=F,data=ejemplo6)
Ejemplo 2: Análisis de la varianza con un factor
Estimación:
Construimos las matrices X e Y:
Y= ejemplo2[,2]
X2=ifelse(ejemplo2$fertilizante=="F2",1,0)
X3=ifelse(ejemplo2$fertilizante=="F3",1,0)
X4=ifelse(ejemplo2$fertilizante=="F4",1,0)
X5=ifelse(ejemplo2$fertilizante=="F5",1,0)
X= cbind(1,X2,X3,X4,X5)
Y## [1] 51 49 50 49 51 50 56 60 56 56 57 48 50 53 44 45 47 48 49 44 43 43 46
## [24] 47 45 46X## X2 X3 X4 X5
## [1,] 1 0 0 0 0
## [2,] 1 0 0 0 0
## [3,] 1 0 0 0 0
## [4,] 1 0 0 0 0
## [5,] 1 0 0 0 0
## [6,] 1 0 0 0 0
## [7,] 1 1 0 0 0
## [8,] 1 1 0 0 0
## [9,] 1 1 0 0 0
## [10,] 1 1 0 0 0
## [11,] 1 1 0 0 0
## [12,] 1 0 1 0 0
## [13,] 1 0 1 0 0
## [14,] 1 0 1 0 0
## [15,] 1 0 1 0 0
## [16,] 1 0 1 0 0
## [17,] 1 0 0 1 0
## [18,] 1 0 0 1 0
## [19,] 1 0 0 1 0
## [20,] 1 0 0 1 0
## [21,] 1 0 0 0 1
## [22,] 1 0 0 0 1
## [23,] 1 0 0 0 1
## [24,] 1 0 0 0 1
## [25,] 1 0 0 0 1
## [26,] 1 0 0 0 1R nos puede generar directamente la matriz X del modelo:
XX=model.matrix(rendimiento~fertilizante,data=ejemplo2)
XX # Comprobamos que es la misma matriz X anterior.## (Intercept) fertilizanteF2 fertilizanteF3 fertilizanteF4 fertilizanteF5
## 1 1 0 0 0 0
## 2 1 0 0 0 0
## 3 1 0 0 0 0
## 4 1 0 0 0 0
## 5 1 0 0 0 0
## 6 1 0 0 0 0
## 7 1 1 0 0 0
## 8 1 1 0 0 0
## 9 1 1 0 0 0
## 10 1 1 0 0 0
## 11 1 1 0 0 0
## 12 1 0 1 0 0
## 13 1 0 1 0 0
## 14 1 0 1 0 0
## 15 1 0 1 0 0
## 16 1 0 1 0 0
## 17 1 0 0 1 0
## 18 1 0 0 1 0
## 19 1 0 0 1 0
## 20 1 0 0 1 0
## 21 1 0 0 0 1
## 22 1 0 0 0 1
## 23 1 0 0 0 1
## 24 1 0 0 0 1
## 25 1 0 0 0 1
## 26 1 0 0 0 1
## attr(,"assign")
## [1] 0 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$fertilizante
## [1] "contr.treatment"Calculamos el estimador de \(\beta\) utilizando directamente cálculo matricial:
solve(t(X)%*%X)%*%t(X)%*%Y## [,1]
## 50
## X2 7
## X3 -2
## X4 -3
## X5 -5Calculamos el estimador de \(\beta\) utilizando el programa “lm” de R:
lm(rendimiento~fertilizante,data=ejemplo2)##
## Call:
## lm(formula = rendimiento ~ fertilizante, data = ejemplo2)
##
## Coefficients:
## (Intercept) fertilizanteF2 fertilizanteF3 fertilizanteF4
## 50 7 -2 -3
## fertilizanteF5
## -5Presentamos un resumen de la estimación:
modelo2=lm(rendimiento~fertilizante,data=ejemplo2)
summary(modelo2)##
## Call:
## lm(formula = rendimiento ~ fertilizante, data = ejemplo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4 -1 0 1 5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.0000 0.8819 56.695 < 2e-16 ***
## fertilizanteF2 7.0000 1.3081 5.351 2.63e-05 ***
## fertilizanteF3 -2.0000 1.3081 -1.529 0.141203
## fertilizanteF4 -3.0000 1.3944 -2.151 0.043236 *
## fertilizanteF5 -5.0000 1.2472 -4.009 0.000636 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.16 on 21 degrees of freedom
## Multiple R-squared: 0.8178, Adjusted R-squared: 0.7831
## F-statistic: 23.57 on 4 and 21 DF, p-value: 1.649e-07El test estadístico \(F\) en este caso lleva a cabo el contraste: \[H_0: \alpha_2=\alpha_3=\dots=\alpha_p=0, \] La hipótesis nula en este caso significa que no hay diferencia entre usar un fertilizante u otro. En forma matricial:
\[H_{0}:\left(\begin{array}{ccccc} 0 & 1 & 0 & \dots & 0\\ 0 & 0 & 1 & \dots & 0\\ 0 & 0 & 0 & \dots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & & 1 \end{array}\right)\left(\begin{array}{c} \vartheta\\ \alpha_{2}\\ \alpha_{3}\\ \vdots\\ \alpha_{p} \end{array}\right)=\left(\begin{array}{c} 0\\ 0\\ 0\\ \vdots\\ 0 \end{array}\right)\]
Es más habitual mostrar este contraste en formato de tabla de análisis de la varianza:
modelo2.aov=aov(rendimiento~fertilizante,data=ejemplo2)
summary(modelo2.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## fertilizante 4 439.9 109.97 23.57 1.65e-07 ***
## Residuals 21 98.0 4.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1o bien:
anova(modelo2)## Analysis of Variance Table
##
## Response: rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## fertilizante 4 439.88 109.971 23.565 1.649e-07 ***
## Residuals 21 98.00 4.667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ejemplo 3: Análisis de la varianza con dos factores (sin interacción)
Estimación:
Construimos las matrices X e Y:
Y= ejemplo3[,3]
XAQ2=ifelse(ejemplo3$agenteQuim=="AQ2",1,0)
XAQ3=ifelse(ejemplo3$agenteQuim=="AQ3",1,0)
XAQ4=ifelse(ejemplo3$agenteQuim=="AQ4",1,0)
XRT2=ifelse(ejemplo3$rolloTela=="RT2",1,0)
XRT3=ifelse(ejemplo3$rolloTela=="RT3",1,0)
XRT4=ifelse(ejemplo3$rolloTela=="RT4",1,0)
XRT5=ifelse(ejemplo3$rolloTela=="RT5",1,0)
X= cbind(1,XAQ2,XAQ3,XAQ4,XRT2,XRT3,XRT4,XRT5)
Y## [1] 73 68 74 71 67 73 67 75 72 70 75 68 78 73 68 73 71 75 75 69X## XAQ2 XAQ3 XAQ4 XRT2 XRT3 XRT4 XRT5
## [1,] 1 0 0 0 0 0 0 0
## [2,] 1 0 0 0 1 0 0 0
## [3,] 1 0 0 0 0 1 0 0
## [4,] 1 0 0 0 0 0 1 0
## [5,] 1 0 0 0 0 0 0 1
## [6,] 1 1 0 0 0 0 0 0
## [7,] 1 1 0 0 1 0 0 0
## [8,] 1 1 0 0 0 1 0 0
## [9,] 1 1 0 0 0 0 1 0
## [10,] 1 1 0 0 0 0 0 1
## [11,] 1 0 1 0 0 0 0 0
## [12,] 1 0 1 0 1 0 0 0
## [13,] 1 0 1 0 0 1 0 0
## [14,] 1 0 1 0 0 0 1 0
## [15,] 1 0 1 0 0 0 0 1
## [16,] 1 0 0 1 0 0 0 0
## [17,] 1 0 0 1 1 0 0 0
## [18,] 1 0 0 1 0 1 0 0
## [19,] 1 0 0 1 0 0 1 0
## [20,] 1 0 0 1 0 0 0 1R nos puede generar también la matriz X del modelo:
XX=model.matrix(resistencia~agenteQuim+rolloTela,data=ejemplo3)
XX # Comprobamos que es la misma matriz X anterior.## (Intercept) agenteQuimAQ2 agenteQuimAQ3 agenteQuimAQ4 rolloTelaRT2
## 1 1 0 0 0 0
## 2 1 0 0 0 1
## 3 1 0 0 0 0
## 4 1 0 0 0 0
## 5 1 0 0 0 0
## 6 1 1 0 0 0
## 7 1 1 0 0 1
## 8 1 1 0 0 0
## 9 1 1 0 0 0
## 10 1 1 0 0 0
## 11 1 0 1 0 0
## 12 1 0 1 0 1
## 13 1 0 1 0 0
## 14 1 0 1 0 0
## 15 1 0 1 0 0
## 16 1 0 0 1 0
## 17 1 0 0 1 1
## 18 1 0 0 1 0
## 19 1 0 0 1 0
## 20 1 0 0 1 0
## rolloTelaRT3 rolloTelaRT4 rolloTelaRT5
## 1 0 0 0
## 2 0 0 0
## 3 1 0 0
## 4 0 1 0
## 5 0 0 1
## 6 0 0 0
## 7 0 0 0
## 8 1 0 0
## 9 0 1 0
## 10 0 0 1
## 11 0 0 0
## 12 0 0 0
## 13 1 0 0
## 14 0 1 0
## 15 0 0 1
## 16 0 0 0
## 17 0 0 0
## 18 1 0 0
## 19 0 1 0
## 20 0 0 1
## attr(,"assign")
## [1] 0 1 1 1 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$agenteQuim
## [1] "contr.treatment"
##
## attr(,"contrasts")$rolloTela
## [1] "contr.treatment"Calculamos el estimador de beta utilizando directamente el cálculo matricial:
solve(t(X)%*%X)%*%t(X)%*%Y## [,1]
## 72.35
## XAQ2 0.80
## XAQ3 1.80
## XAQ4 2.00
## XRT2 -5.00
## XRT3 2.00
## XRT4 -0.75
## XRT5 -5.00Calculamos el estimador de beta utilizando el programa “lm” de R:
lm(resistencia~agenteQuim+rolloTela,data=ejemplo3)##
## Call:
## lm(formula = resistencia ~ agenteQuim + rolloTela, data = ejemplo3)
##
## Coefficients:
## (Intercept) agenteQuimAQ2 agenteQuimAQ3 agenteQuimAQ4 rolloTelaRT2
## 72.35 0.80 1.80 2.00 -5.00
## rolloTelaRT3 rolloTelaRT4 rolloTelaRT5
## 2.00 -0.75 -5.00Resumen de la estimación y análisis de la varianza:
modelo3=lm(resistencia~agenteQuim+rolloTela,data=ejemplo3)
summary(modelo3)##
## Call:
## lm(formula = resistencia ~ agenteQuim + rolloTela, data = ejemplo3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3500 -0.7375 -0.3500 0.7000 1.8500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 72.3500 0.8524 84.873 < 2e-16 ***
## agenteQuimAQ2 0.8000 0.8524 0.938 0.366507
## agenteQuimAQ3 1.8000 0.8524 2.112 0.056374 .
## agenteQuimAQ4 2.0000 0.8524 2.346 0.036968 *
## rolloTelaRT2 -5.0000 0.9531 -5.246 0.000206 ***
## rolloTelaRT3 2.0000 0.9531 2.098 0.057699 .
## rolloTelaRT4 -0.7500 0.9531 -0.787 0.446584
## rolloTelaRT5 -5.0000 0.9531 -5.246 0.000206 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.348 on 12 degrees of freedom
## Multiple R-squared: 0.8863, Adjusted R-squared: 0.82
## F-statistic: 13.36 on 7 and 12 DF, p-value: 8.335e-05anova(modelo3)## Analysis of Variance Table
##
## Response: resistencia
## Df Sum Sq Mean Sq F value Pr(>F)
## agenteQuim 3 12.95 4.317 2.3761 0.1211
## rolloTela 4 157.00 39.250 21.6055 2.059e-05 ***
## Residuals 12 21.80 1.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ejemplo 4: Análisis de la varianza con tres factores (e interacciones)
Estimación:
Construimos las matrices X e Y:
Y= ejemplo4[,1]
# Modelo sin interacciones:
X=model.matrix(resistencia~concFibra+tiempo+presion,data=ejemplo4)
X## (Intercept) concFibra4% concFibra8% tiempo4 h. presion
## 1 1 0 0 0 400
## 2 1 0 0 0 500
## 3 1 0 0 0 650
## 4 1 0 0 1 400
## 5 1 0 0 1 500
## 6 1 0 0 1 650
## 7 1 0 0 0 400
## 8 1 0 0 0 500
## 9 1 0 0 0 650
## 10 1 0 0 1 400
## 11 1 0 0 1 500
## 12 1 0 0 1 650
## 13 1 1 0 0 400
## 14 1 1 0 0 500
## 15 1 1 0 0 650
## 16 1 1 0 1 400
## 17 1 1 0 1 500
## 18 1 1 0 1 650
## 19 1 1 0 0 400
## 20 1 1 0 0 500
## 21 1 1 0 0 650
## 22 1 1 0 1 400
## 23 1 1 0 1 500
## 24 1 1 0 1 650
## 25 1 0 1 0 400
## 26 1 0 1 0 500
## 27 1 0 1 0 650
## 28 1 0 1 1 400
## 29 1 0 1 1 500
## 30 1 0 1 1 650
## 31 1 0 1 0 400
## 32 1 0 1 0 500
## 33 1 0 1 0 650
## 34 1 0 1 1 400
## 35 1 0 1 1 500
## 36 1 0 1 1 650
## attr(,"assign")
## [1] 0 1 1 2 3
## attr(,"contrasts")
## attr(,"contrasts")$concFibra
## [1] "contr.treatment"
##
## attr(,"contrasts")$tiempo
## [1] "contr.treatment"# OJO!!: la presión es un factor y se ha leido como variable numérica. Debemos modificar la sintaxis:
X=model.matrix(resistencia~concFibra+tiempo+factor(presion),data=ejemplo4)
X## (Intercept) concFibra4% concFibra8% tiempo4 h. factor(presion)500
## 1 1 0 0 0 0
## 2 1 0 0 0 1
## 3 1 0 0 0 0
## 4 1 0 0 1 0
## 5 1 0 0 1 1
## 6 1 0 0 1 0
## 7 1 0 0 0 0
## 8 1 0 0 0 1
## 9 1 0 0 0 0
## 10 1 0 0 1 0
## 11 1 0 0 1 1
## 12 1 0 0 1 0
## 13 1 1 0 0 0
## 14 1 1 0 0 1
## 15 1 1 0 0 0
## 16 1 1 0 1 0
## 17 1 1 0 1 1
## 18 1 1 0 1 0
## 19 1 1 0 0 0
## 20 1 1 0 0 1
## 21 1 1 0 0 0
## 22 1 1 0 1 0
## 23 1 1 0 1 1
## 24 1 1 0 1 0
## 25 1 0 1 0 0
## 26 1 0 1 0 1
## 27 1 0 1 0 0
## 28 1 0 1 1 0
## 29 1 0 1 1 1
## 30 1 0 1 1 0
## 31 1 0 1 0 0
## 32 1 0 1 0 1
## 33 1 0 1 0 0
## 34 1 0 1 1 0
## 35 1 0 1 1 1
## 36 1 0 1 1 0
## factor(presion)650
## 1 0
## 2 0
## 3 1
## 4 0
## 5 0
## 6 1
## 7 0
## 8 0
## 9 1
## 10 0
## 11 0
## 12 1
## 13 0
## 14 0
## 15 1
## 16 0
## 17 0
## 18 1
## 19 0
## 20 0
## 21 1
## 22 0
## 23 0
## 24 1
## 25 0
## 26 0
## 27 1
## 28 0
## 29 0
## 30 1
## 31 0
## 32 0
## 33 1
## 34 0
## 35 0
## 36 1
## attr(,"assign")
## [1] 0 1 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$concFibra
## [1] "contr.treatment"
##
## attr(,"contrasts")$tiempo
## [1] "contr.treatment"
##
## attr(,"contrasts")$`factor(presion)`
## [1] "contr.treatment"Calculamos el estimador de beta utilizando directamente el cálculo matricial:
solve(t(X)%*%X)%*%t(X)%*%Y## [,1]
## (Intercept) 197.44444444
## concFibra4% -0.70833333
## concFibra8% -1.12500000
## tiempo4 h. 1.50000000
## factor(presion)500 -0.09166667
## factor(presion)650 1.50833333Calculamos el estimador de beta utilizando el programa “lm” de R:
lm(resistencia~concFibra+tiempo+factor(presion),data=ejemplo4)##
## Call:
## lm(formula = resistencia ~ concFibra + tiempo + factor(presion),
## data = ejemplo4)
##
## Coefficients:
## (Intercept) concFibra4% concFibra8%
## 197.44444 -0.70833 -1.12500
## tiempo4 h. factor(presion)500 factor(presion)650
## 1.50000 -0.09167 1.50833Modelo con interacciones:
X=model.matrix(resistencia~concFibra+tiempo+factor(presion)+concFibra*tiempo
+concFibra*factor(presion)+tiempo*factor(presion),data=ejemplo4)
X## (Intercept) concFibra4% concFibra8% tiempo4 h. factor(presion)500
## 1 1 0 0 0 0
## 2 1 0 0 0 1
## 3 1 0 0 0 0
## 4 1 0 0 1 0
## 5 1 0 0 1 1
## 6 1 0 0 1 0
## 7 1 0 0 0 0
## 8 1 0 0 0 1
## 9 1 0 0 0 0
## 10 1 0 0 1 0
## 11 1 0 0 1 1
## 12 1 0 0 1 0
## 13 1 1 0 0 0
## 14 1 1 0 0 1
## 15 1 1 0 0 0
## 16 1 1 0 1 0
## 17 1 1 0 1 1
## 18 1 1 0 1 0
## 19 1 1 0 0 0
## 20 1 1 0 0 1
## 21 1 1 0 0 0
## 22 1 1 0 1 0
## 23 1 1 0 1 1
## 24 1 1 0 1 0
## 25 1 0 1 0 0
## 26 1 0 1 0 1
## 27 1 0 1 0 0
## 28 1 0 1 1 0
## 29 1 0 1 1 1
## 30 1 0 1 1 0
## 31 1 0 1 0 0
## 32 1 0 1 0 1
## 33 1 0 1 0 0
## 34 1 0 1 1 0
## 35 1 0 1 1 1
## 36 1 0 1 1 0
## factor(presion)650 concFibra4%:tiempo4 h. concFibra8%:tiempo4 h.
## 1 0 0 0
## 2 0 0 0
## 3 1 0 0
## 4 0 0 0
## 5 0 0 0
## 6 1 0 0
## 7 0 0 0
## 8 0 0 0
## 9 1 0 0
## 10 0 0 0
## 11 0 0 0
## 12 1 0 0
## 13 0 0 0
## 14 0 0 0
## 15 1 0 0
## 16 0 1 0
## 17 0 1 0
## 18 1 1 0
## 19 0 0 0
## 20 0 0 0
## 21 1 0 0
## 22 0 1 0
## 23 0 1 0
## 24 1 1 0
## 25 0 0 0
## 26 0 0 0
## 27 1 0 0
## 28 0 0 1
## 29 0 0 1
## 30 1 0 1
## 31 0 0 0
## 32 0 0 0
## 33 1 0 0
## 34 0 0 1
## 35 0 0 1
## 36 1 0 1
## concFibra4%:factor(presion)500 concFibra8%:factor(presion)500
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0
## 7 0 0
## 8 0 0
## 9 0 0
## 10 0 0
## 11 0 0
## 12 0 0
## 13 0 0
## 14 1 0
## 15 0 0
## 16 0 0
## 17 1 0
## 18 0 0
## 19 0 0
## 20 1 0
## 21 0 0
## 22 0 0
## 23 1 0
## 24 0 0
## 25 0 0
## 26 0 1
## 27 0 0
## 28 0 0
## 29 0 1
## 30 0 0
## 31 0 0
## 32 0 1
## 33 0 0
## 34 0 0
## 35 0 1
## 36 0 0
## concFibra4%:factor(presion)650 concFibra8%:factor(presion)650
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0
## 7 0 0
## 8 0 0
## 9 0 0
## 10 0 0
## 11 0 0
## 12 0 0
## 13 0 0
## 14 0 0
## 15 1 0
## 16 0 0
## 17 0 0
## 18 1 0
## 19 0 0
## 20 0 0
## 21 1 0
## 22 0 0
## 23 0 0
## 24 1 0
## 25 0 0
## 26 0 0
## 27 0 1
## 28 0 0
## 29 0 0
## 30 0 1
## 31 0 0
## 32 0 0
## 33 0 1
## 34 0 0
## 35 0 0
## 36 0 1
## tiempo4 h.:factor(presion)500 tiempo4 h.:factor(presion)650
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 1 0
## 6 0 1
## 7 0 0
## 8 0 0
## 9 0 0
## 10 0 0
## 11 1 0
## 12 0 1
## 13 0 0
## 14 0 0
## 15 0 0
## 16 0 0
## 17 1 0
## 18 0 1
## 19 0 0
## 20 0 0
## 21 0 0
## 22 0 0
## 23 1 0
## 24 0 1
## 25 0 0
## 26 0 0
## 27 0 0
## 28 0 0
## 29 1 0
## 30 0 1
## 31 0 0
## 32 0 0
## 33 0 0
## 34 0 0
## 35 1 0
## 36 0 1
## attr(,"assign")
## [1] 0 1 1 2 3 3 4 4 5 5 5 5 6 6
## attr(,"contrasts")
## attr(,"contrasts")$concFibra
## [1] "contr.treatment"
##
## attr(,"contrasts")$tiempo
## [1] "contr.treatment"
##
## attr(,"contrasts")$`factor(presion)`
## [1] "contr.treatment"Calculamos el estimador de beta utilizando directamente el cálculo matricial:
solve(t(X)%*%X)%*%t(X)%*%Y## [,1]
## (Intercept) 196.5500000
## concFibra4% 0.9833333
## concFibra8% 0.5666667
## tiempo4 h. 1.7000000
## factor(presion)500 0.4500000
## factor(presion)650 2.6500000
## concFibra4%:tiempo4 h. -1.1166667
## concFibra8%:tiempo4 h. -0.8833333
## concFibra4%:factor(presion)500 -1.4500000
## concFibra8%:factor(presion)500 -1.9000000
## concFibra4%:factor(presion)650 -1.9500000
## concFibra8%:factor(presion)650 -1.8500000
## tiempo4 h.:factor(presion)500 1.1500000
## tiempo4 h.:factor(presion)650 0.2500000Calculamos el estimador de beta utilizando el programa “lm” de R:
lm(resistencia~concFibra+tiempo+factor(presion)+concFibra*tiempo
+concFibra*factor(presion)+tiempo*factor(presion),data=ejemplo4)##
## Call:
## lm(formula = resistencia ~ concFibra + tiempo + factor(presion) +
## concFibra * tiempo + concFibra * factor(presion) + tiempo *
## factor(presion), data = ejemplo4)
##
## Coefficients:
## (Intercept) concFibra4%
## 196.5500 0.9833
## concFibra8% tiempo4 h.
## 0.5667 1.7000
## factor(presion)500 factor(presion)650
## 0.4500 2.6500
## concFibra4%:tiempo4 h. concFibra8%:tiempo4 h.
## -1.1167 -0.8833
## concFibra4%:factor(presion)500 concFibra8%:factor(presion)500
## -1.4500 -1.9000
## concFibra4%:factor(presion)650 concFibra8%:factor(presion)650
## -1.9500 -1.8500
## tiempo4 h.:factor(presion)500 tiempo4 h.:factor(presion)650
## 1.1500 0.2500Comparación de ambos modelos:
modelo4.sin=lm(resistencia~concFibra+tiempo+factor(presion),data=ejemplo4)
modelo4.con=lm(resistencia~concFibra+tiempo+factor(presion)+concFibra*tiempo
+concFibra*factor(presion)+tiempo*factor(presion),data=ejemplo4)
summary(modelo4.sin)##
## Call:
## lm(formula = resistencia ~ concFibra + tiempo + factor(presion),
## data = ejemplo4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.44444 -0.63194 0.02222 0.45139 1.76389
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 197.44444 0.32422 608.987 < 2e-16 ***
## concFibra4% -0.70833 0.32422 -2.185 0.0369 *
## concFibra8% -1.12500 0.32422 -3.470 0.0016 **
## tiempo4 h. 1.50000 0.26472 5.666 3.55e-06 ***
## factor(presion)500 -0.09167 0.32422 -0.283 0.7793
## factor(presion)650 1.50833 0.32422 4.652 6.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7942 on 30 degrees of freedom
## Multiple R-squared: 0.7147, Adjusted R-squared: 0.6671
## F-statistic: 15.03 on 5 and 30 DF, p-value: 2.091e-07summary(modelo4.con)##
## Call:
## lm(formula = resistencia ~ concFibra + tiempo + factor(presion) +
## concFibra * tiempo + concFibra * factor(presion) + tiempo *
## factor(presion), data = ejemplo4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00000 -0.49167 0.01667 0.39583 0.96667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.5500 0.3888 505.480 < 2e-16 ***
## concFibra4% 0.9833 0.5091 1.931 0.06641 .
## concFibra8% 0.5667 0.5091 1.113 0.27770
## tiempo4 h. 1.7000 0.4648 3.658 0.00138 **
## factor(presion)500 0.4500 0.5091 0.884 0.38631
## factor(presion)650 2.6500 0.5091 5.205 3.21e-05 ***
## concFibra4%:tiempo4 h. -1.1167 0.5091 -2.193 0.03914 *
## concFibra8%:tiempo4 h. -0.8833 0.5091 -1.735 0.09672 .
## concFibra4%:factor(presion)500 -1.4500 0.6235 -2.325 0.02966 *
## concFibra8%:factor(presion)500 -1.9000 0.6235 -3.047 0.00591 **
## concFibra4%:factor(presion)650 -1.9500 0.6235 -3.127 0.00490 **
## concFibra8%:factor(presion)650 -1.8500 0.6235 -2.967 0.00712 **
## tiempo4 h.:factor(presion)500 1.1500 0.5091 2.259 0.03414 *
## tiempo4 h.:factor(presion)650 0.2500 0.5091 0.491 0.62825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6235 on 22 degrees of freedom
## Multiple R-squared: 0.871, Adjusted R-squared: 0.7948
## F-statistic: 11.43 on 13 and 22 DF, p-value: 6.146e-07anova(modelo4.con)## Analysis of Variance Table
##
## Response: resistencia
## Df Sum Sq Mean Sq F value Pr(>F)
## concFibra 2 7.7639 3.8819 9.9847 0.0008211 ***
## tiempo 1 20.2500 20.2500 52.0850 3.123e-07 ***
## factor(presion) 2 19.3739 9.6869 24.9158 2.224e-06 ***
## concFibra:tiempo 2 2.0817 1.0408 2.6771 0.0910703 .
## concFibra:factor(presion) 4 6.0911 1.5228 3.9167 0.0150436 *
## tiempo:factor(presion) 2 2.1950 1.0975 2.8229 0.0810470 .
## Residuals 22 8.5533 0.3888
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(modelo4.sin,modelo4.con)## Analysis of Variance Table
##
## Model 1: resistencia ~ concFibra + tiempo + factor(presion)
## Model 2: resistencia ~ concFibra + tiempo + factor(presion) + concFibra *
## tiempo + concFibra * factor(presion) + tiempo * factor(presion)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 18.9211
## 2 22 8.5533 8 10.368 3.3334 0.0119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo sólo con interacción entre concentración de Fibra y presión:
modelo4b.con=lm(resistencia~concFibra+tiempo+factor(presion)
+concFibra*factor(presion),data=ejemplo4)
anova(modelo4.con,modelo4b.con)## Analysis of Variance Table
##
## Model 1: resistencia ~ concFibra + tiempo + factor(presion) + concFibra *
## tiempo + concFibra * factor(presion) + tiempo * factor(presion)
## Model 2: resistencia ~ concFibra + tiempo + factor(presion) + concFibra *
## factor(presion)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 22 8.5533
## 2 26 12.8300 -4 -4.2767 2.75 0.05395 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(modelo4.sin,modelo4b.con)## Analysis of Variance Table
##
## Model 1: resistencia ~ concFibra + tiempo + factor(presion)
## Model 2: resistencia ~ concFibra + tiempo + factor(presion) + concFibra *
## factor(presion)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 18.921
## 2 26 12.830 4 6.0911 3.0859 0.03322 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(modelo4b.con)## Analysis of Variance Table
##
## Response: resistencia
## Df Sum Sq Mean Sq F value Pr(>F)
## concFibra 2 7.7639 3.8819 7.8668 0.00213 **
## tiempo 1 20.2500 20.2500 41.0366 8.711e-07 ***
## factor(presion) 2 19.3739 9.6869 19.6306 6.370e-06 ***
## concFibra:factor(presion) 4 6.0911 1.5228 3.0859 0.03322 *
## Residuals 26 12.8300 0.4935
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modelo4b.con)##
## Call:
## lm(formula = resistencia ~ concFibra + tiempo + factor(presion) +
## concFibra * factor(presion), data = ejemplo4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6750 -0.3438 0.0000 0.4000 1.4250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.6500 0.3702 531.152 < 2e-16 ***
## concFibra4% 0.4250 0.4967 0.856 0.4000
## concFibra8% 0.1250 0.4967 0.252 0.8033
## tiempo4 h. 1.5000 0.2342 6.406 8.71e-07 ***
## factor(presion)500 1.0250 0.4967 2.064 0.0492 *
## factor(presion)650 2.7750 0.4967 5.587 7.20e-06 ***
## concFibra4%:factor(presion)500 -1.4500 0.7025 -2.064 0.0491 *
## concFibra8%:factor(presion)500 -1.9000 0.7025 -2.705 0.0119 *
## concFibra4%:factor(presion)650 -1.9500 0.7025 -2.776 0.0101 *
## concFibra8%:factor(presion)650 -1.8500 0.7025 -2.634 0.0140 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7025 on 26 degrees of freedom
## Multiple R-squared: 0.8065, Adjusted R-squared: 0.7395
## F-statistic: 12.04 on 9 and 26 DF, p-value: 3.153e-07Ejemplo 6: Análisis de la covarianza
Estimación:
Construimos las matrices X e Y:
Y= ejemplo6[,1]
XT2=ifelse(ejemplo6$tipo=="T2",1,0)
XT3=ifelse(ejemplo6$tipo=="T3",1,0)
XT4=ifelse(ejemplo6$tipo=="T4",1,0)
X=cbind(1,ejemplo6$velocidad,XT2,XT3,XT4)Modelo sin interacciones:
XX=model.matrix(voltaje~velocidad+tipo,data=ejemplo6)
XX## (Intercept) velocidad tipoT2 tipoT3 tipoT4
## 1 1 6336 0 0 0
## 2 1 7099 0 0 0
## 3 1 8026 0 0 0
## 4 1 6230 0 0 0
## 5 1 5369 0 0 0
## 6 1 8343 0 0 0
## 7 1 6522 1 0 0
## 8 1 7310 1 0 0
## 9 1 7974 1 0 0
## 10 1 8501 1 0 0
## 11 1 6646 0 1 0
## 12 1 7384 0 1 0
## 13 1 8000 0 1 0
## 14 1 8545 0 1 0
## 15 1 6755 0 0 1
## 16 1 7362 0 0 1
## 17 1 7934 0 0 1
## 18 1 8554 0 0 1
## attr(,"assign")
## [1] 0 1 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$tipo
## [1] "contr.treatment"Calculamos el estimador de beta utilizando directamente el cálculo matricial:
solve(t(X)%*%X)%*%t(X)%*%Y## [,1]
## -1.6954187213
## 0.0005700677
## XT2 -0.5313416048
## XT3 -0.9220361396
## XT4 -1.2063116472Calculamos el estimador de beta utilizando el programa “lm” de R:
lm(voltaje~velocidad+tipo,data=ejemplo6)##
## Call:
## lm(formula = voltaje ~ velocidad + tipo, data = ejemplo6)
##
## Coefficients:
## (Intercept) velocidad tipoT2 tipoT3 tipoT4
## -1.6954187 0.0005701 -0.5313416 -0.9220361 -1.2063116Mostramos la inferencia asociada al modelo:
summary(lm(voltaje~velocidad+tipo,data=ejemplo6))##
## Call:
## lm(formula = voltaje ~ velocidad + tipo, data = ejemplo6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16610 -0.05424 0.00382 0.06670 0.14879
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.695e+00 2.269e-01 -7.471 4.69e-06 ***
## velocidad 5.701e-04 3.219e-05 17.709 1.75e-10 ***
## tipoT2 -5.313e-01 7.657e-02 -6.940 1.02e-05 ***
## tipoT3 -9.220e-01 7.721e-02 -11.942 2.21e-08 ***
## tipoT4 -1.206e+00 7.728e-02 -15.609 8.44e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1137 on 13 degrees of freedom
## Multiple R-squared: 0.9716, Adjusted R-squared: 0.9629
## F-statistic: 111.3 on 4 and 13 DF, p-value: 6.42e-10
Contrastes de hipótesis lineales generales
Supongamos que en el ejemplo 2 queremos comparar entre sí los rendimientos de los distintos fertilizantes. Utilizaremos para ello la función glht (general linear hypothesis testing) del paquete multcomp:
require(multcomp)
modelo2.aov=aov(rendimiento~ fertilizante,data=ejemplo2)
compara=glht(modelo2,linfct=mcp(fertilizante=c("F1-F2=0",
"F1-F3=0",
"F1-F4=0",
"F1-F5=0",
"F2-F3=0",
"F2-F4=0",
"F2-F5=0")))
summary(compara)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: lm(formula = rendimiento ~ fertilizante, data = ejemplo2)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## F1 - F2 == 0 -7.000 1.308 -5.351 < 0.001 ***
## F1 - F3 == 0 2.000 1.308 1.529 0.49488
## F1 - F4 == 0 3.000 1.394 2.151 0.19548
## F1 - F5 == 0 5.000 1.247 4.009 0.00374 **
## F2 - F3 == 0 9.000 1.366 6.587 < 0.001 ***
## F2 - F4 == 0 10.000 1.449 6.901 < 0.001 ***
## F2 - F5 == 0 12.000 1.308 9.174 < 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(compara)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: lm(formula = rendimiento ~ fertilizante, data = ejemplo2)
##
## Quantile = 2.8606
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## F1 - F2 == 0 -7.0000 -10.7419 -3.2581
## F1 - F3 == 0 2.0000 -1.7419 5.7419
## F1 - F4 == 0 3.0000 -0.9889 6.9889
## F1 - F5 == 0 5.0000 1.4322 8.5678
## F2 - F3 == 0 9.0000 5.0917 12.9083
## F2 - F4 == 0 10.0000 5.8546 14.1454
## F2 - F5 == 0 12.0000 8.2581 15.7419