大卫哈里斯提供了一个很好的答案，但由于该问题继续被编辑，也许这将有助于查看他的解决方案的细节。以下分析的重点是：
为此，让我们使用指定的公式创建一些实际数据，以便我们可以评估解决方案的准确性。这是通过以下方式完成的R：
set.seed(17)
n.names <- 1000
groupSize <- 3.5
n.cases <- 5 * n.names # Should exceed n.names
cv <- 0.10 # Must be 0 or greater
groupSize <- 3.5 # Must be greater than 0
proficiency <- round(rgamma(n.names, 20, scale=5)); hist(proficiency)
在这些初始步骤中，我们：
为随机数生成器设置一个种子，这样任何人都可以准确地重现结果。
指定有多少个工人n.names。
用 规定每组的预期工人人数groupSize。
指定有多少个案例（观察）可用n.cases。（稍后其中一些将被淘汰，因为它们随机发生，与我们合成劳动力中的任何一个工人都不对应。）
根据每组工作“熟练程度”的总和，安排工作量与预测的数量随机不同。的值cv是典型的比例变化；例如，这里给出的对应于典型的 10% 变化（在少数情况下可能超过 30%）。0.10
创建一支由具有不同工作能力的人员组成的劳动力队伍。这里给出的计算参数proficiency在最好和最差的工人之间创造了超过 4:1 的范围（根据我的经验，这对于技术和专业工作来说甚至可能有点窄，但对于日常制造工作来说可能很宽）。
有了这个合成劳动力，让我们模拟他们的工作。这相当于schedule为每个观察 ) 以反映将不可避免地发生的变化。（如果根本没有变化，我们会将这个问题提交给数学网站，其中受访者可以指出这个问题只是一组联立线性方程组，可以精确地解决熟练程度。）1
schedule <- matrix(rbinom(n.cases * n.names, 1, groupSize/n.names), nrow=n.cases)
schedule <- schedule[apply(schedule, 1, sum) > 0, ]
work <- round(schedule %*% proficiency * exp(rnorm(dim(schedule)[1], -cv^2/2, cv)))
hist(work)
我发现将所有工作组数据放入单个数据框中进行分析很方便，但要保持工作值分开：
data <- data.frame(schedule)
这是我们从真实数据开始的地方：data我们将使用（或）编码的工人分组和数组schedule中观察到的工作输出。work
不幸的是，如果某些工作人员总是配对，则R' 的lm过程将简单地失败并出现错误。 我们应该首先检查这种配对。一种方法是在时间表中找到完全相关的工人：
correlations <- cor(data)
outer(names(data), names(data), paste)[which(upper.tri(correlations) &
correlations >= 0.999999)]
输出将列出成对的工人对：这可以用来将这些工人组合成组，因为至少我们可以估计每个组的生产力，如果不是其中的个人的话。我们希望它只是吐出来character(0)。让我们假设它确实如此。
上述解释中隐含的一个微妙点是，所做工作的变化是相乘的，而不是相加的。这是现实的：在绝对规模上，一大群工人的产出变化将大于较小群体的变化。因此，我们将通过使用加权最小二乘而不是普通最小二乘来获得更好的估计。在这个特定模型中使用的最佳权重是工作量的倒数。（如果某些工作量为零，我会通过添加少量来避免除以零来伪造它。）
fit <- lm(work ~ . + 0, data=data, weights=1/(max(work)/10^3+work))
fit.sum <- summary(fit)
这应该只需要一两秒钟。
在继续之前，我们应该执行一些适合的诊断测试。虽然在这里讨论这些会让我们走得太远，但R产生有用诊断的一个命令是
plot(fit)
（这将需要几秒钟：这是一个大型数据集！）
尽管这几行代码完成了所有工作，并为每个工人提供了估计的熟练程度，但我们不想扫描所有 1000 行输出——至少不是马上。 让我们使用图形来显示结果。
fit.coef <- coef(fit.sum)
results <- cbind(fit.coef[, c("Estimate", "Std. Error")],
Actual=proficiency,
Difference=fit.coef[, "Estimate"] - proficiency,
Residual=(fit.coef[, "Estimate"] - proficiency)/fit.coef[, "Std. Error"])
hist(results[, "Residual"])
plot(results[, c("Actual", "Estimate")])
直方图（下图左下方）是估计和实际熟练程度之间的差异，表示为估计标准误差的倍数。对于一个好的过程，这些值几乎总是介于和之间，并且对称分布在周围。但是，在涉及 1000 名工人的情况下，我们完全希望看到其中一些标准化差异从甚至−220340. 这正是这里的情况：直方图与人们希望的一样漂亮。（当然这很好：毕竟这些是模拟数据。但是对称性证实了权重正确地完成了它们的工作。使用错误的权重往往会产生不对称的直方图。）
散点图（图右下方）直接将估计的熟练度与实际熟练度进行比较。当然这在现实中是不可能的，因为我们不知道实际的熟练程度：这就是计算机模拟的力量。观察：
如果工作中没有随机变化（设置cv=0并重新运行代码以查看这一点），散点图将是一条完美的对角线。所有的估计都是完全准确的。因此，这里看到的分散反映了这种变化。
有时，估计值与实际值相差甚远。例如，在 (110, 160) 附近有一个点，估计熟练程度比实际熟练程度高约 50%。这在任何大批量数据中几乎是不可避免的。如果估计将用于个人基础，例如用于评估工人，请记住这一点。总的来说，这些估计可能很好，但如果工作效率的变化是由于任何个人无法控制的原因造成的，那么对于少数工人来说，估计将是错误的：有些太高，有些太低。并且无法准确判断谁受到了影响。
以下是在此过程中生成的四个图。
最后，请注意，这种回归方法很容易适应控制可能与团队生产力相关的其他变量。这些可能包括小组规模、每项工作的持续时间、时间变量、每个小组经理的因素等等。只需将它们作为回归中的附加变量包括在内。