When adding means of different sample sizes with known standard deviations, how can I compute the standard deviation of the combined mean? I have three means and standard deviations that represent subscales of a factor. Each subscale has a different number of items. I am summing the three subscale means together to form a factor, and am wondering if there is any way to obtain a standard deviation for the resulting factor score. Below is my example, with means listed in the first column and SDs in the second:
Sub1 (18 items)     49.80    8.16
Sub2 (13 items)     42.10    5.50
Sub3 (14 items)     36.27    7.91
Fact (45 items)     128.17   XXX
 A: If the "factor" score of interest ($Z$) can be view as adding together the three subscales of this particular interest and each subscale is assumed to be normally distributed and independent of the others, then the we simply see $Z$ as being the linear combination of three independent Gaussian variables $A$, $B$ and $C$.
The general case is such that if our variable of interest $Z$ is the mixture of $K$ independent Gaussian variables, ie: $Z = \sum_k^{K} w_k N(\mu_{k}, \sigma_k^2)$, then $Z$ is itself normal-distributed as $Z\sim N(\sum_k^K w_k \mu_{k}, \sum_k^K (w_k \sigma_k)^2)$. The Wikipedia article on the sum of normally distributed random variables covers nicely if you want to read more details.
Now, using the above for the case of three Gaussian random variables $A$, $B$ and $C$, distributed as $N(\mu_A,\sigma_A^2)$, $N(\mu_B,\sigma_B^2)$ and $N(\mu_C,\sigma_C^2)$ respectively and assuming that $w_A = w_B = w_C = 1$, then we can immediately substitute $\mu_{A} = 49.80$, $\sigma_A = 8.16$,
$\mu_{B} = 42.1$, $\sigma_b = 5.5$, $\mu_{C} = 36.27$ and $\sigma_C = 7.91$, to get:
$\mu_Z = w_A \mu_A + w_B \mu_B + w_C \mu_C = 128.17$ (as you have derived yourself) and
$\sigma_Z^2 = w_A^2 \sigma_A^2 + w_B^2 \sigma_B^2 + w_C^2 \sigma_C^2 = 159.4037$ (or $\sigma_Z = 12.62552$)
Here is a quick R snippet showing this in action:
N = 10^6
mu = c(49.8, 42.1, 36.27)
sigma = c(8.16, 5.5, 7.91) 
set.seed(710)
A =  rnorm(N, mu[1], sd = sigma[1])
B =  rnorm(N, mu[2], sd = sigma[2])
C =  rnorm(N, mu[3], sd = sigma[3])
Z = A+B+C
round(c(mean(Z), sd(Z)),3)
# 128.170  12.625 as expected

