estimatr
estimatr copied to clipboard
bad defaulting to match pair?
I have data like this:
> table(dat$TA1, dat$block)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
0 31 6 31 8 1 17 92 32 73 42 11 15 27 14 50 22 7
1 11 3 12 3 1 8 24 10 32 10 2 4 23 3 18 7 3
And output like this. Not appeared to reduce to matched pair, because of block 5 I suppose, and only 16 degrees of freedom, which seems to suggest just 2 units per block being used...
difference_in_means(ncb_wh ~ TA1, blocks = block, data=dat)
Design: Matched-pair
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
TA1 -0.01343043 0.05930286 -0.2264719 0.8237023 -0.1391469 0.112286 16
The issue I think is that when it see a block with a matched pair it is concluding that this is a matched pair design. We have discussed this before I know, but surely better behavior is to either drop such small blocks, or combine them with other small blocks.
I think we should probably error in the above case because we don't know how to calculate the standard error. Agree the output is for sure wrong.
I think the user should have to explicitly combine the small blocks or explicitly drop them. Could put that suggestion in the error maybe?
+1
Maybe consult with Peter or Cyrus
There are options, eg pool up and include warning. This would be like a local homoskedacticity assumption I suppose but could imagine being preferable to tanking. Maybe a force option to allow this.
Could easily imagine this happens by chance with some regularity in some schemes and choking seems undesirable
On Mon, 8 Oct 2018, 23:26 Graeme Blair, [email protected] wrote:
+1
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/DeclareDesign/estimatr/issues/260#issuecomment-427984360, or mute the thread https://github.com/notifications/unsubscribe-auth/AMJO_VlRyb2Gn3Fn4ey9lE2TcyS8TEsAks5ui8LpgaJpZM4XDdkQ .
I guess to me choking seems highly informative -- you may learn at the design stage you unintentionally have small blocks and can modify your blocking scheme. At the analysis stage, you may realize you want to change the blocking scheme and our collapsing wouldn't be exactly what you would choose. (In order to select for them, there are many options and I imagine it will be hard in general to select the optimal pooling up algorithm.
If we do allow this, this seems like a surprising behavior -- changing the blocks variable within the estimator! -- so think we should have a more verbose option like force_collapse_blocks to be very clear we are doing this.
I see that I'd be happy with a tank default and force collapse option
On Tue, 9 Oct 2018, 00:01 Graeme Blair, [email protected] wrote:
I guess to me choking seems highly informative -- you may learn at the design stage you unintentionally have small blocks and can modify your blocking scheme. At the analysis stage, you may realize you want to change the blocking scheme and our collapsing wouldn't be exactly what you would choose. (In order to select for them, there are many options and I imagine it will be hard in general to select the optimal pooling up algorithm.
If we do allow this, this seems like a surprising behavior -- changing the blocks variable within the estimator! -- so think we should have a more verbose option like force_collapse_blocks to be very clear we are doing this.
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/DeclareDesign/estimatr/issues/260#issuecomment-427993446, or mute the thread https://github.com/notifications/unsubscribe-auth/AMJO_cTleAMb2j6svIuBitQEwRclOBinks5ui8tNgaJpZM4XDdkQ .
I remember we discussed this previously. Whatever we discussed led us to implement the following warning that amounts to @macartan's original suggestion in this thread:
https://github.com/DeclareDesign/estimatr/blob/5850b9e1ac3dc1e582bc9cf3144f733b1aeda4e8/R/estimatr_difference_in_means.R#L305
See this and the next few lines. You should be seeing the warning, and indeed the test here properly generates this warning: https://github.com/DeclareDesign/estimatr/blob/5850b9e1ac3dc1e582bc9cf3144f733b1aeda4e8/tests/testthat/test-difference-in-means.R#L50
The following is from that test, see block "49"
> table(dat$Z, dat$bad_mp)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
> difference_in_means(Y ~ Z, blocks = bad_mp, data = dat)
Design: Matched-pair
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
Z 0.05100903 0.2016612 0.2529442 0.8013915 -0.354458 0.4564761 48
Warning message:
In difference_in_means(Y ~ Z, blocks = bad_mp, data = dat) :
Some `blocks` have two units/`clusters` while other blocks have more units/`clusters`. As standard variance estimates cannot be computed within blocks with two units, we use the matched pairs estimator of the variance.
Not wt machine now but if I remember I got the warning alright but the behavior was bad. I had all big blocks and only one block with a pair. It used matched pair which i think meant only two units per blog (which ones?) and other data seemed to be nit used
On Thu, 11 Oct 2018, 01:25 Luke Sonnet, [email protected] wrote:
I remember we had discussed this previously. Whatever we discussed led us to implement the following warning:
https://github.com/DeclareDesign/estimatr/blob/5850b9e1ac3dc1e582bc9cf3144f733b1aeda4e8/R/estimatr_difference_in_means.R#L305
See this and the next few lines. You should be seeing the warning, and indeed the test here properly generates this warning: https://github.com/DeclareDesign/estimatr/blob/5850b9e1ac3dc1e582bc9cf3144f733b1aeda4e8/tests/testthat/test-difference-in-means.R#L50 .
table(dat$Z, dat$bad_mp)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
difference_in_means(Y ~ Z, blocks = bad_mp, data = dat) Design: Matched-pair Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF Z 0.05100903 0.2016612 0.2529442 0.8013915 -0.354458 0.4564761 48 Warning message: In difference_in_means(Y ~ Z, blocks = bad_mp, data = dat) : Some
blocks
have two units/clusters
while other blocks have more units/clusters
. As standard variance estimates cannot be computed within blocks with two units, we use the matched pairs estimator of the variance.— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/DeclareDesign/estimatr/issues/260#issuecomment-428766103, or mute the thread https://github.com/notifications/unsubscribe-auth/AMJO_ZBCRgBgsKPad4V2y1EF8eyTfO9zks5ujoHegaJpZM4XDdkQ .
Apologies if we have gone over before but in mixed cases is it not possible to use matched pair approach for all matched pairs and block approach for larger blocks and then combine variances, given independence between blocks
On Thu, 11 Oct 2018, 06:42 Macartan Humphreys, [email protected] wrote:
Not wt machine now but if I remember I got the warning alright but the behavior was bad. I had all big blocks and only one block with a pair. It used matched pair which i think meant only two units per blog (which ones?) and other data seemed to be nit used
On Thu, 11 Oct 2018, 01:25 Luke Sonnet, [email protected] wrote:
I remember we had discussed this previously. Whatever we discussed led us to implement the following warning:
https://github.com/DeclareDesign/estimatr/blob/5850b9e1ac3dc1e582bc9cf3144f733b1aeda4e8/R/estimatr_difference_in_means.R#L305
See this and the next few lines. You should be seeing the warning, and indeed the test here properly generates this warning: https://github.com/DeclareDesign/estimatr/blob/5850b9e1ac3dc1e582bc9cf3144f733b1aeda4e8/tests/testthat/test-difference-in-means.R#L50 .
table(dat$Z, dat$bad_mp)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
difference_in_means(Y ~ Z, blocks = bad_mp, data = dat) Design: Matched-pair Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF Z 0.05100903 0.2016612 0.2529442 0.8013915 -0.354458 0.4564761 48 Warning message: In difference_in_means(Y ~ Z, blocks = bad_mp, data = dat) : Some
blocks
have two units/clusters
while other blocks have more units/clusters
. As standard variance estimates cannot be computed within blocks with two units, we use the matched pairs estimator of the variance.— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/DeclareDesign/estimatr/issues/260#issuecomment-428766103, or mute the thread https://github.com/notifications/unsubscribe-auth/AMJO_ZBCRgBgsKPad4V2y1EF8eyTfO9zks5ujoHegaJpZM4XDdkQ .
My intuition is that that should be correct. Do we need to wait for a proof?
On advice from Frederik, here's a piece on this problem: https://rss.onlinelibrary.wiley.com/doi/full/10.1111/rssb.12290
Sorry for the delay on this, and thanks for the paper @acoppock.
Macartan, I only mentioned that we had discussed this because I remember coming up with this solution, which I'm sure we agreed at the time was unsatisfactory. I also wanted to ensure you were getting the warning. However, your fear that it is only using two units per block is unfounded. It will apply the classic matched-pairs estimator, which is simply the sample variance of the block-level estimates
(1 / (J * (J-1))) \sum^J_j (\tauhat_j - \tauhat)^2
Unfortunately, this is not weighted by the size of your groups, and throws away the precision that could be gained by looking within those larger blocks.
As for your suggestion going forward, have no sense as to whether combining the MP estimator in each paired block and the classical estimator within other blocks would be appropriate. I've never seen that proposed anywhere, and the linked paper doesn't address it directly. That being said, my uncertainty comes from ignorance, and uncertainty about why you can combine two very different within block estimators.
To be clear, what you are proposing is the following, correct?
(1 / (J * (J - 1))) \sum^J_j (J * N_j / N)^2 \varhat(\tauhat_j)
where J is the number of blocks, N_j is the number of units in block j, N is the total number of units, and \varhat(\tauhat_j) is the estimated variance within the block which for MP is (\tauhat_j - \tauhat)^2 and for larger blocks is the classic variance estimator using the sample variance for treated and control units.
Is that right? I presume we use that correction (1 / (J * (J-1))) because otherwise there is no correction for using the estimate of tau directly in its variance estimation. That will be "overcorrecting" the variance from the larger blocks but I see no other solution.
As to the paper provided by Alex, this seems like a great general approach, although it seems to be largely speaking finely stratified experiments, where every block has one treatment condition with at most one unit, rather than the case where some blocks have treatment conditions with only one unit. Furthermore, because most of the gains from that estimator are from covariate adjustment, I'm not sure it applies in our case where we don't do covariate adjustment. However, the estimator S^2(Q_1) proposed for the case with groups of varying sizes has several desirable properties and could be tested against:
- Our current implementation, which is to just get the sample variance of the block estimates (the unweighted matched-pairs estimator)
- a slightly smarter estimator, which uses the sample variance of the block estimates but re-weights by the size of the blocks
- the (as far as I can tell) ad-hoc estimator suggested in this thread that I tried to clarify above.
Hi Luke -- sorry I haven't gotten to respond to this yet. My thought was to imagine this as two experiments. One among all those groups with >1 per stratum (group A) and the other with only 1 per stratum (group B): estimate \varhat(\tauhat_A) and \varhat(\tauhat_B) using the usual corrections for A and B and then combine ; now:
\tauhat_AB = w_A \tauhat_A + w_B \tauhat_B
so
var(\tauhat_AB) = w_A^2 var(\tauhat_A) + w_B^2 var(\tauhat_B)
I believe then that this last step is right; as this is the same as what we do in the usual blocked approach where we combine the varhats from different blocks:
varhat(\tauhat_AB) = w_A^2 varhat(\tauhat_A) + w_B^2 varhat(\tauhat_B)
On Thu, Oct 18, 2018 at 6:36 PM Luke Sonnet [email protected] wrote:
Sorry for the delay on this, and thanks for the paper @acoppock https://github.com/acoppock.
Macartan, I only mentioned that we had discussed this because I remember coming up with this solution, which I'm sure we agreed at the time was unsatisfactory. I also wanted to ensure you were getting the warning. However, your fear that it is only using two units per block is unfounded. It will apply the classic matched-pairs estimator, which is simply the sample variance of the block-level estimates (1 / (J * (J-1))) \sum^J_j (\tauhat_j - \tauhat)^2
Unfortunately, this is not weighted by the size of your groups, and throws away the precision that could be gained by looking within those larger blocks.
As for your suggestion going forward, have no sense as to whether combining the MP estimator in each paired block and the classical estimator within other blocks would be appropriate. I've never seen that proposed anywhere, and the linked paper doesn't address it directly. That being said, my uncertainty comes from ignorance, and uncertainty about why you can combine two very different within block estimators.
To be clear, what you are proposing is the following, correct?
(1 / (J * (J - 1))) \sum^J_j (J * N_j / N)^2 \varhat(\tauhat_j)
where J is the number of blocks, N_j is the number of units in block j, N is the total number of units, and \varhat(\tauhat_j) is the estimated variance within the block which for MP is (\tauhat_j - \tauhat)^2 and for larger blocks is the classic variance estimator using the sample variance for treated and control units.
Is that right? I presume we use that correction (1 / (J * (J-1))) because otherwise there is no correction for using the estimate of tau directly in its variance estimation. That will be "overcorrecting" the variance from the larger blocks but I see no other solution.
As to the paper provided by Alex, this seems like a great general approach, although it seems to be largely speaking finely stratified experiments, where every block has one treatment condition with at most one unit, rather than the case where some blocks have treatment conditions with only one unit. Furthermore, because most of the gains from that estimator are from covariate adjustment, I'm not sure it applies in our case where we don't do covariate adjustment. However, the estimator S^2(Q_1) proposed for the case with groups of varying sizes has several desirable properties and could be tested against:
- Our current implementation, which is to just get the sample variance of the block estimates (the unweighted matched-pairs estimator)
- a slightly smarter estimator, which uses the sample variance of the block estimates but re-weights by the size of the blocks
- the (as far as I can tell) ad-hoc estimator suggested in this thread that I tried to clarify above.
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/DeclareDesign/estimatr/issues/260#issuecomment-431187110, or mute the thread https://github.com/notifications/unsubscribe-auth/AMJO_WOx_t9MDN5fLPoMQa7akXrde8Daks5umQJygaJpZM4XDdkQ .
No worries Macartan. I just want to clarify:
Group B could have strata with 2 units or even 3+ units, as long as within each strata there is one unit that is alone it its treatment condition. Correct?
And then w_A is defined in the usual way n_A / N, where n_A is the total number of units in all group A strata and N is the total number of units in the sample. Correct?
Yes on both (by stratum I meant as defined by block and treatment status).
To be clear though I am just proposing this for discussion. It seems natural but I have not worked through it
On Tue, Oct 23, 2018 at 11:51 AM Luke Sonnet [email protected] wrote:
No worries Macartan. I just want to clarify:
Group B could have strata with 2 units or even 3+ units, as long as within each strata there is one unit that is alone it its treatment condition. Correct?
And then w_A is defined in the usual way n_A / N, where n_A is the total number of units in all group A strata and N is the total number of units in the sample. Correct?
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/DeclareDesign/estimatr/issues/260#issuecomment-432301517, or mute the thread https://github.com/notifications/unsubscribe-auth/AMJO_VVR0KI4HKx4uSe15PJ0znLNlYFKks5unzr-gaJpZM4XDdkQ .