<p><a href="http://arxiv.org/ftp/arxiv/papers/1007/1007.1012.pdf" rel="nofollow noreferrer">Rudmin (2010)</a> states that exact variance of pooled data set is the mean of the variances plus the variance of the means. <a href="https://stackoverflow.com/questions/9222056/existing-function-to-combine-standard-deviations-in-r?rq=1">flodel has already provided an answer and function</a> that gives similar values to Rudmin's statement. Using Rudmin's data set and <a href="https://stackoverflow.com/questions/9222056/existing-function-to-combine-standard-deviations-in-r?rq=1">flodel's function</a> based on <a href="http://en.wikipedia.org/wiki/Standard_deviation#Combining_standard_deviations" rel="nofollow noreferrer">Wikipedia</a>:</p>

<pre><code>df &lt;- data.frame(mean = c(30.66667, 31.14286, 40.33333), variance = c(8.555555, 13.26531, 1.555555), n = c(6,7,3))

grand.sd   &lt;- function(S, M, N) {sqrt(weighted.mean(S^2 + M^2, N) -
                                      weighted.mean(M, N)^2)}

grand.sd(sqrt(df$variance), df$mean, df$n)^2 

#[1] 22.83983 = Dp variance in Rudmin (2010). 
</code></pre>

<p>However this solution gives slightly different values compared to the <a href="http://www.talkstats.com/showthread.php/14523-An-average-of-standard-deviations" rel="nofollow noreferrer">function 5.38 from Headrick (2010)</a> (unless there is a mistake somewhere):</p>

<pre><code>dat &lt;- data.frame(variable = c(rep("x", 2), rep("y", 3)), replicate = c(1,2,1,2,3),
mean = c(3.4, 2.5, 6.5, 5.7, 5.1), sd = c(1.2, 0.7, 2.4, 4.0, 3.5),
n = c(3,3,5,4,6))

x &lt;- subset(dat, variable == "x")

((x$n[1]^2)*(x$sd[1]^2)+
(x$n[2]^2)*(x$sd[2]^2)-
(x$n[2])*(x$sd[1]^2) -
(x$n[2])*(x$sd[2]^2) -
(x$n[1])*(x$sd[1]^2) -
(x$n[1])*(x$sd[2]^2) +
(x$n[1])*(x$n[2])*(x$sd[1]^2) +
(x$n[1])*(x$n[2])*(x$sd[2]^2) +
(x$n[1])*(x$n[2])*(x$mean[1] - x$mean[2])^2)/
((x$n[1] + x$n[2] - 1)*(x$n[1] + x$n[2]))

#[1] 1.015

grand.sd(x$sd, x$mean, x$n)^2

#[1] 1.1675
</code></pre>

<p>To answer my own question,  the desired <code>data.frame</code> would be acquired followingly:</p>

<pre><code>library(plyr)
ddply(dat, c("variable"), function(dat) c(mean=with(dat,weighted.mean(mean, n)),  sd = with(dat, grand.sd(sd, mean, n))))   

  variable     mean       sd
1        x 2.950000 1.080509
2        y 5.726667 3.382793
</code></pre>
Standard deviation of combined data
<p>I have a data set with mean values, standard deviations and n. One of the variables has an equal sample size, while the sample size for the other one varies.</p>

<pre><code>dat &lt;- data.frame(variable = c(rep("x", 2), rep("y", 3)), replicate = c(1,2,1,2,3),
mean = c(3.4, 2.5, 6.5, 5.7, 5.1), sd = c(1.2, 0.7, 2.4, 4.0, 3.5),
n = c(3,3,5,4,6))
</code></pre>

<p>I need to combine <code>x</code> and <code>y</code> variables and am trying to find a code-sparing way to calculate combined standard deviation for instance using by <code>aggregate</code> function. <a href="http://www.talkstats.com/showthread.php/14523-An-average-of-standard-deviations" rel="nofollow noreferrer">The equation for combined standard deviation</a> is following: </p>

<p><img src="https://i.stack.imgur.com/bcfph.gif" alt="enter image description here"></p>

<p>And for unequal sample sizes (<a href="http://www.talkstats.com/showthread.php/14523-An-average-of-standard-deviations" rel="nofollow noreferrer">same source</a>):</p>

<p><img src="https://i.stack.imgur.com/qcoht.gif" alt="enter image description here"></p>

<p>My combined data frame should look like this:</p>

<pre><code>variable    mean    sd
x           2.95    sd_x
y           5.76    sd_y
</code></pre>

<p><strong>How to make a function in R that calculates the combined standard deviation?</strong> Or alternatively, if there is a package designed for this, it counts as an answer too =)</p>