# ANOVA

## 2901 days ago by ab277290

effect <- data.frame(geneexp = c(100, 99, 87, 98, 95, 94, 83, 92, 102, 80, 86, 83, 84, 82, 80, 83, 77, 78, 90, 76, 90, 76, 74, 85), agegroup = factor(rep(rep(1:6, rep(4,6)),1)), dose = factor(rep(c("D1", "D2", "D3", "D4"), c(1,1,1,1))) )
effect
 ``` geneexp agegroup dose 1 100 1 D1 2 99 1 D2 3 87 1 D3 4 98 1 D4 5 95 2 D1 6 94 2 D2 7 83 2 D3 8 92 2 D4 9 102 3 D1 10 80 3 D2 11 86 3 D3 12 83 3 D4 13 84 4 D1 14 82 4 D2 15 80 4 D3 16 83 4 D4 17 77 5 D1 18 78 5 D2 19 90 5 D3 20 76 5 D4 21 90 6 D1 22 76 6 D2 23 74 6 D3 24 85 6 D4```
with(effect, tapply(geneexp, agegroup, mean))
 ``` 1 2 3 4 5 6 96.00 91.00 87.75 82.25 80.25 81.25 ```
with(effect, tapply(geneexp, dose, mean))
 ``` D1 D2 D3 D4 91.33333 84.83333 83.33333 86.16667 ```
with(effect, tapply(geneexp, agegroup, var))
 ``` 1 2 3 4 5 6 36.666667 30.000000 96.250000 2.916667 42.916667 56.916667 ```
with(effect, tapply(geneexp, dose, var))
 ``` D1 D2 D3 D4 92.66667 88.16667 32.66667 59.76667 ```
boxplot(split(effect\$geneexp, effect\$agegroup), xlab="Age Group", ylab="Gene Expression") dev.off()
 ```null device 1 ```
boxplot(split(effect\$geneexp, effect\$dose), xlab="Age Group", ylab="Gene Expression") dev.off()
 ```null device 1 ```
fitexp <- lm(formula = geneexp ~ agegroup + dose, data=effect) fitexp
 ```Call: lm(formula = geneexp ~ agegroup + dose, data = effect) Coefficients: (Intercept) agegroup2 agegroup3 agegroup4 agegroup5 agegroup6 doseD2 100.917 -5.000 -8.250 -13.750 -15.750 -14.750 -6.500 doseD3 doseD4 -8.000 -5.167 ```
summary(fitexp)
 ```Call: lm(formula = geneexp ~ agegroup + dose, data = effect) Residuals: Min 1Q Median 3Q Max -8.1667 -4.0417 0.0833 2.6458 12.8333 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 100.917 3.806 26.513 5.14e-14 *** agegroup2 -5.000 4.395 -1.138 0.27312 agegroup3 -8.250 4.395 -1.877 0.08009 . agegroup4 -13.750 4.395 -3.129 0.00690 ** agegroup5 -15.750 4.395 -3.584 0.00272 ** agegroup6 -14.750 4.395 -3.356 0.00433 ** doseD2 -6.500 3.589 -1.811 0.09016 . doseD3 -8.000 3.589 -2.229 0.04150 * doseD4 -5.167 3.589 -1.440 0.17048 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.216 on 15 degrees of freedom Multiple R-squared: 0.6341, Adjusted R-squared: 0.439 F-statistic: 3.25 on 8 and 15 DF, p-value: 0.02353 ```
an <- anova(fitexp) an
 ```Analysis of Variance Table Response: geneexp Df Sum Sq Mean Sq F value Pr(>F) agegroup 5 786.83 157.367 4.0733 0.01551 * dose 3 217.50 72.500 1.8766 0.17691 Residuals 15 579.50 38.633 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ```
an[]
 ```Analysis of Variance Table Response: geneexp Df Sum Sq Mean Sq F value Pr(>F) agegroup 5 786.83 157.367 4.0733 0.01551 * dose 3 217.50 72.500 1.8766 0.17691 Residuals 15 579.50 38.633 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ```
fitexp2 <- lm(formula = geneexp ~ agegroup * dose, data=effect) summary(fitexp2)
 ```Call: lm(formula = geneexp ~ agegroup * dose, data = effect) Residuals: ALL 24 residuals are 0: no residual degrees of freedom! Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.000e+02 NA NA NA agegroup2 -5.000e+00 NA NA NA agegroup3 2.000e+00 NA NA NA agegroup4 -1.600e+01 NA NA NA agegroup5 -2.300e+01 NA NA NA agegroup6 -1.000e+01 NA NA NA doseD2 -1.000e+00 NA NA NA doseD3 -1.300e+01 NA NA NA doseD4 -2.000e+00 NA NA NA agegroup2:doseD2 7.864e-14 NA NA NA agegroup3:doseD2 -2.100e+01 NA NA NA agegroup4:doseD2 -1.000e+00 NA NA NA agegroup5:doseD2 2.000e+00 NA NA NA agegroup6:doseD2 -1.300e+01 NA NA NA agegroup2:doseD3 1.000e+00 NA NA NA agegroup3:doseD3 -3.000e+00 NA NA NA agegroup4:doseD3 9.000e+00 NA NA NA agegroup5:doseD3 2.600e+01 NA NA NA agegroup6:doseD3 -3.000e+00 NA NA NA agegroup2:doseD4 -1.000e+00 NA NA NA agegroup3:doseD4 -1.700e+01 NA NA NA agegroup4:doseD4 1.000e+00 NA NA NA agegroup5:doseD4 1.000e+00 NA NA NA agegroup6:doseD4 -3.000e+00 NA NA NA Residual standard error: NaN on 0 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: NaN F-statistic: NaN on 23 and 0 DF, p-value: NA ```
library(MASS) stepAIC(fitexp2, direction="f", scope = list(upper = ~agegroup*dose, lower = ~1))
 ```Start: AIC=-Inf geneexp ~ agegroup * dose Call: lm(formula = geneexp ~ agegroup * dose, data = effect) Coefficients: (Intercept) agegroup2 agegroup3 agegroup4 agegroup5 1.000e+02 -5.000e+00 2.000e+00 -1.600e+01 -2.300e+01 agegroup6 doseD2 doseD3 doseD4 agegroup2:doseD2 -1.000e+01 -1.000e+00 -1.300e+01 -2.000e+00 7.864e-14 agegroup3:doseD2 agegroup4:doseD2 agegroup5:doseD2 agegroup6:doseD2 agegroup2:doseD3 -2.100e+01 -1.000e+00 2.000e+00 -1.300e+01 1.000e+00 agegroup3:doseD3 agegroup4:doseD3 agegroup5:doseD3 agegroup6:doseD3 agegroup2:doseD4 -3.000e+00 9.000e+00 2.600e+01 -3.000e+00 -1.000e+00 agegroup3:doseD4 agegroup4:doseD4 agegroup5:doseD4 agegroup6:doseD4 -1.700e+01 1.000e+00 1.000e+00 -3.000e+00 ```