Calculate within and between variances and confidence intervals in R

You have four groups of three observations:

> run1 = c(9.85, 9.95, 10.00)
> run2 = c(9.90, 8.80, 9.50)
> run3 = c(11.20, 11.10, 9.80)
> run4 = c(9.70, 10.10, 10.00)
> runs = c(run1, run2, run3, run4)
> runs
 [1]  9.85  9.95 10.00  9.90  8.80  9.50 11.20 11.10  9.80  9.70 10.10 10.00

Make some labels:

> n = rep(3, 4)
> group = rep(1:4, n)
> group
 [1] 1 1 1 2 2 2 3 3 3 4 4 4

Calculate within-run stats:

> withinRunStats = function(x) c(sum = sum(x), mean = mean(x), var = var(x), n = length(x))
> tapply(runs, group, withinRunStats)
$`1`
         sum         mean          var            n 
29.800000000  9.933333333  0.005833333  3.000000000 

$`2`
  sum  mean   var     n 
28.20  9.40  0.31  3.00 

$`3`
  sum  mean   var     n 
32.10 10.70  0.61  3.00 

$`4`
        sum        mean         var           n 
29.80000000  9.93333333  0.04333333  3.00000000 

You can do some ANOVA here:

> data = data.frame(y = runs, group = factor(group))
> data
       y group
1   9.85     1
2   9.95     1
3  10.00     1
4   9.90     2
5   8.80     2
6   9.50     2
7  11.20     3
8  11.10     3
9   9.80     3
10  9.70     4
11 10.10     4
12 10.00     4

> fit = lm(runs ~ group, data)
> fit

Call:
lm(formula = runs ~ group, data = data)

Coefficients:
(Intercept)       group2       group3       group4  
  9.933e+00   -5.333e-01    7.667e-01   -2.448e-15 

> anova(fit)
Analysis of Variance Table

Response: runs
          Df  Sum Sq Mean Sq F value  Pr(>F)  
group      3 2.57583 0.85861  3.5437 0.06769 .
Residuals  8 1.93833 0.24229                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

> degreesOfFreedom = anova(fit)[, "Df"]
> names(degreesOfFreedom) = c("treatment", "error")
> degreesOfFreedom
treatment     error 
        3         8

Error or within-group variance:

> anova(fit)["Residuals", "Mean Sq"]
[1] 0.2422917

Treatment or between-group variance:

> anova(fit)["group", "Mean Sq"]
[1] 0.8586111

This should give you enough confidence to do confidence intervals.


If you want to apply a function (such as var) across a factor such as Run or Rep, you can use tapply:

> with(variance, tapply(Value, Run, var))
          1           2           3           4 
0.005833333 0.310000000 0.610000000 0.043333333 
> with(variance, tapply(Value, Rep, var))
          1          2          3 
0.48562500 0.88729167 0.05583333 

I'm going to take a crack at this when I have more time, but meanwhile, here's the dput() for Kiar's data structure:

structure(list(Run = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4), Rep = c(1, 
2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3), Value = c(9.85, 9.95, 10, 9.9, 
8.8, 9.5, 11.2, 11.1, 9.8, 9.7, 10.1, 10)), .Names = c("Run", 
"Rep", "Value"), row.names = c(NA, -12L), class = "data.frame")

... in case you'd like to take a quick shot at it.

Tags:

R

Statistics