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.