ANOVA

2275 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