-
Notifications
You must be signed in to change notification settings - Fork 0
/
3_fertility_robustness2.html
865 lines (815 loc) · 67.9 KB
/
3_fertility_robustness2.html
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<title>3_fertility_robustness2.utf8.md</title>
<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/bootstrap.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/default.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type="text/css">
h1 {
font-size: 34px;
}
h1.title {
font-size: 38px;
}
h2 {
font-size: 30px;
}
h3 {
font-size: 24px;
}
h4 {
font-size: 18px;
}
h5 {
font-size: 16px;
}
h6 {
font-size: 12px;
}
.table th:not([align]) {
text-align: left;
}
</style>
</head>
<body>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
code {
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
img {
max-width:100%;
height: auto;
}
.tabbed-pane {
padding-top: 12px;
}
.html-widget {
margin-bottom: 20px;
}
button.code-folding-btn:focus {
outline: none;
}
</style>
<style type="text/css">
/* padding for bootstrap navbar */
body {
padding-top: 51px;
padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar) */
.section h1 {
padding-top: 56px;
margin-top: -56px;
}
.section h2 {
padding-top: 56px;
margin-top: -56px;
}
.section h3 {
padding-top: 56px;
margin-top: -56px;
}
.section h4 {
padding-top: 56px;
margin-top: -56px;
}
.section h5 {
padding-top: 56px;
margin-top: -56px;
}
.section h6 {
padding-top: 56px;
margin-top: -56px;
}
</style>
<script>
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark it active
menuAnchor.parent().addClass('active');
// if it's got a parent navbar menu mark it active as well
menuAnchor.closest('li.dropdown').addClass('active');
});
</script>
<div class="container-fluid main-container">
<!-- tabsets -->
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
</script>
<!-- code folding -->
<div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">Fertility diary</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="1_power_analysis.html">Power analysis</a>
</li>
<li>
<a href="1_researcher_df_analysis.html">Researcher degree of freedom simulation</a>
</li>
<li>
<a href="1_wrangle_data.html">Data wrangling</a>
</li>
<li>
<a href="2_descriptives.html">Descriptives</a>
</li>
<li>
<a href="3_fertility_as_prereg.html">Preregistered analyses</a>
</li>
<li>
<a href="3_fertility_robustness.html">Robustness tests</a>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<div class="fluid-row" id="header">
</div>
<div id="robustness-checks-2-of-ovulatory-changes" class="section level1">
<h1>Robustness checks 2 <small>of ovulatory changes</small></h1>
<p>These were conducted after further communication with Steve Gangestad, after publication.</p>
<p><span style="background:red;width:20px;height:20px;display:inline-block;"></span> Cycling women (not on hormonal birth control)</p>
<p><span style="background:black;width:20px;height:20px;display:inline-block;"></span> Women on hormonal birth control</p>
<pre class="r"><code>library(knitr)
opts_chunk$set(fig.width = 8, fig.height = 8, cache = T, warning = T, message = F, cache = F)</code></pre>
<pre class="r"><code>source("0_helpers.R")
load("full_data.rdata")
diary = diary %>%
mutate(
included = included_all,
fertile = if_else(is.na(prc_stirn_b_squished), prc_stirn_b_backward_inferred, prc_stirn_b_squished),
contraceptive_methods = factor(contraceptive_method, levels =
c("barrier_or_abstinence", "fertility_awareness", "none", "hormonal")),
relationship_status_clean = factor(relationship_status_clean),
cohabitation = factor(cohabitation),
certainty_menstruation = as.numeric(as.character(certainty_menstruation)),
partner_st_vs_lt = partner_attractiveness_shortterm - partner_attractiveness_longterm
) %>% group_by(person) %>%
mutate(
fertile_mean = mean(fertile, na.rm = T)
)
opts_chunk$set(warning = F)
library(Cairo)
library(effects)
opts_chunk$set(dev = "CairoPNG")
options(width=8000)
diary$age_group = cut(diary$age,c(18,20,25,30,35,70), include.lowest = T)</code></pre>
<p>Steve Gangestad claimed (via email) that when we wrote the following in the paper, we made a mistake.</p>
<p><strong>Edit</strong>: Steve Gangestad has asked me to clarify that we tested several terms when testing the interaction. As he explains, this is the better explanation for the p value difference than differences between Wald and LR ratio tests - and he’s right, I should have this clearer, but rushed a response. He also says, it underpowers the test to include a test for the interaction term ST:LT:fertile and ST:LT:fertile:HC. I am not sure if he agrees with our inclusion of the ST:fertile:HC interaction as part of the test (i.e., moderation should be absent for HC users); he may disagree. I added tests without the ST:LT:fertile interaction here. I’ve lightly edited this text to clarify the points, the <a href="https://github.com/rubenarslan/ovulatory_shifts/commits/master/3_fertility_robustness2.html">changeset</a> can be viewed here. We can only say that we thought we were testing the predicted model, but even testing the model he now said was predicted, there is too much uncertainty to conclude evidence for moderation. He also writes that 2000 naturally cycling women would be needed to power the test for a moderator at level 2 that attenuates a level 1 predictor at 80%. I am working on a power analysis for the planned collaborative study (as he knows), and given that we claimed no absence of effects, I think we do not actually disagree, but this also means (and he seems to agree) that previous work reporting such moderations was doomed to give us misleading answers. So, we’re left with a theory that makes predictions that reasonable people can disagree a lot about how to test specifically.</p>
<p>In results:</p>
<blockquote>
<p>Further, none of the predicted moderation patterns turned significant when adding more women, and using slightly different items for the partner attractiveness moderator variables did not change the pattern.</p>
</blockquote>
<p>In discussion:</p>
<blockquote>
<p>There was insufficient evidence for moderation of male mate retention behavior, extra-pair, or in-pair desire by the partner’s attractiveness (no matter whether it was assessed as relative to self, sexual, or physical), as predicted by the good genes ovulatory shift hypothesis. Although some patterns descriptively pointed in the predicted direction, none of the predicted patterns were significant, and some were opposite to our predictions. Because only 144 naturally cycling women remained for our preregistered analyses, statistical power may have been insufficient to detect plausible moderation effect sizes. However, we found no evidence for moderation effects in the more inclusive sample of our robustness tests either. Although our sample sizes are bigger than many published studies that reported a moderation effect (Haselton & Gangestad, 2006; Pillsworth & Haselton, 2006b), we would ideally prefer to exceed their power by a wider margin owing to winner’s curse, that is, effect sizes being overestimated through selection and publication bias.</p>
</blockquote>
<blockquote>
<p>Although further tests should be conducted, the good genes ovulatory shift hypothesis could be wrong, given that we could not replicate previously reported moderators.</p>
</blockquote>
<blockquote>
<p>In contrast to previous reports, we found no evidence that sexual desire shifted more strongly among women who deemed their partner less sex- ually attractive.</p>
</blockquote>
<p>In one model, we allow two moderators of the fertility effect on extra-pair desire, namely partner’s long-term and short-term attractiveness, <a href="3_fertility_robustness.html#partners-short-term-vs-long-term-attractiveness">summarised here</a>. In that model summary, for one of the moderator specifications (ST vs. LT attractiveness), Steve Gangestad saw a t value of -2.7 for the hypothesis-relevant interaction <code>fertile:partner_attractiveness_shortterm</code>. He argues that this invalidates our statement in the results section, that “none of the predicted moderation patterns turned significant when adding more women” and hence our discussion statement that we found “no evidence for moderation effects in the more inclusive sample of our robustness tests”.</p>
<p>The confusion here is understandable, but we think we were consistent in how we presented and interpreted the evidence. We did not report a Wald t-test with Satterthwaite’s degrees of freedom method for the moderator analyses (as computed by the <code>lmerTest</code> package) to test the question of moderation. Instead, we formulated a baseline model without the focal interaction terms and compared the moderation model to that using the <code>anova</code> function, which performs a likelihood ratio test.</p>
<p>This likelihood ratio test yielded a reported p value of .026, which was not significant according to the threshold of .01 we set (ignoring for the moment the issues with interpreting p values in our non-preregistered robustness checks). We reported the model comparison that we thought was the appropriate test of the prediction and we reported this consistently in the preregistered analyses and the robustness checks. It’s clearly the only p value we reported for moderator analyses in our supplement. Still, I noticed one real (but ultimately unimportant) mistake in our analysis, and took this discussion with Steve as a prompt to examine this moderation in more detail.</p>
<p>The baseline model we specified in the publication did not include the interaction between partners’ long-term attractiveness and fertile. If moderation by long-term attractiveness is not the focal test, as I understand Steve Gangestad to argue (and I think I agree), it should be omitted. Additionally, if the theory does not predict an interaction ST:LT:fertile (which wasn’t clear to me given that some studies used difference score ST-LT which are inappropriate, but seem to map to including ST:LT:fertile).</p>
<div id="model-comparisons-extra-pair-desire-and-behaviour" class="section level2">
<h2>Model comparisons extra-pair desire and behaviour</h2>
<pre class="r"><code>epd1_m0_old = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)</code></pre>
<p>We can specify this new baseline model here.</p>
<pre class="r"><code>epd1_m0 = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)</code></pre>
<p>Now, the model we want to test. We specified an interaction between ST and LT attractiveness, a hormonal contraception dummy and the fertility predictor. So, two terms are added to the baseline model, the interaction between ST:LT:HC:fertile and ST:HC:fertile.</p>
<pre class="r"><code>epd1_m = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)
epd1_m2 = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + (partner_attractiveness_shortterm + partner_attractiveness_longterm) * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)</code></pre>
<p>In the original robustness checks, we did not report p values for coefficients, because we weren’t seeking to make inferences about the covariates, but asking whether the moderator helps us explain the data better, for which conducted a likelihood ratio test. For this, we reported one focal p value against a baseline comparison. However, I’ll do so here, to compare the two different p values we can get.</p>
<p>Are coefficients different from zero?</p>
<pre class="r"><code>summary(epd1_m)</code></pre>
<pre><code>## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + (1 | person)
## Data: diary
##
## REML criterion at convergence: 48514
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.285 -0.557 -0.148 0.403 7.987
##
## Random effects:
## Groups Name Variance Std.Dev.
## person (Intercept) 0.289 0.537
## Residual 0.320 0.566
## Number of obs: 26680, groups: person, 1054
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.81689 0.04690 1312.71869 38.74 < 2e-16 ***
## partner_attractiveness_longterm -0.13238 0.03233 1114.48399 -4.09 0.00004544 ***
## includedhorm_contra -0.09586 0.03937 1253.39361 -2.43 0.01503 *
## partner_attractiveness_shortterm -0.04705 0.02979 1102.58699 -1.58 0.11456
## fertile 0.16947 0.03654 25909.66540 4.64 0.00000353 ***
## menstruationpre -0.08939 0.01729 25913.25563 -5.17 0.00000024 ***
## menstruationyes -0.06950 0.01632 26009.82271 -4.26 0.00002061 ***
## fertile_mean -0.08072 0.20816 1438.55478 -0.39 0.69822
## partner_attractiveness_longterm:includedhorm_contra -0.02470 0.04123 1110.42509 -0.60 0.54926
## includedhorm_contra:partner_attractiveness_shortterm 0.02739 0.03886 1091.40471 0.70 0.48105
## partner_attractiveness_longterm:partner_attractiveness_shortterm 0.03371 0.02380 1111.71749 1.42 0.15695
## partner_attractiveness_shortterm:fertile -0.08784 0.03203 26012.67475 -2.74 0.00610 **
## partner_attractiveness_longterm:fertile 0.06901 0.03534 26107.65069 1.95 0.05085 .
## includedhorm_contra:fertile -0.17245 0.04607 26011.54183 -3.74 0.00018 ***
## includedhorm_contra:menstruationpre 0.06851 0.02220 25909.80477 3.09 0.00203 **
## includedhorm_contra:menstruationyes 0.08475 0.02138 25989.74398 3.96 0.00007397 ***
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm 0.00761 0.03716 1109.36373 0.20 0.83778
## partner_attractiveness_longterm:partner_attractiveness_shortterm:fertile -0.01932 0.02617 26036.70540 -0.74 0.46028
## includedhorm_contra:partner_attractiveness_shortterm:fertile 0.05629 0.04106 25972.25932 1.37 0.17040
## partner_attractiveness_longterm:includedhorm_contra:fertile -0.00301 0.04481 26065.14013 -0.07 0.94638
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm:fertile 0.02688 0.04051 25982.80289 0.66 0.50705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>And what does the likelihood ratio test say?</p>
<p>This is the old way of testing it, with an arguably lenient baseline. With our threshold of .01, not significant, but not far from the threshold.</p>
<pre class="r"><code>anova(epd1_m0_old, epd1_m)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0_old: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m0_old: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m0_old: included * (menstruation + fertile) + fertile_mean + (1 |
## epd1_m0_old: person)
## epd1_m: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m: partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m: included * fertile + included * (menstruation + fertile) +
## epd1_m: fertile_mean + (1 | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0_old 16 48448 48579 -24208 48416
## epd1_m 23 48446 48634 -24200 48400 16 7 0.026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>But this is with the more stringent baseline, where we take into account variance explained by LT:HC:fertile interaction already.</p>
<pre class="r"><code>anova(epd1_m0, epd1_m)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m0: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m0: partner_attractiveness_longterm * included * fertile + included *
## epd1_m0: (menstruation + fertile) + fertile_mean + (1 | person)
## epd1_m: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m: partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m: included * fertile + included * (menstruation + fertile) +
## epd1_m: fertile_mean + (1 | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0 18 48445 48592 -24204 48409
## epd1_m 23 48446 48634 -24200 48400 9.28 5 0.098 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>What about the re-specification, where we omit the terms <code>LT:ST:fertile:HC</code> which Steve Gangestad argues aren’t predicted by theory.</p>
<pre class="r"><code>anova(epd1_m0, epd1_m2)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m0: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m0: partner_attractiveness_longterm * included * fertile + included *
## epd1_m0: (menstruation + fertile) + fertile_mean + (1 | person)
## epd1_m2: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m2: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m2: (partner_attractiveness_shortterm + partner_attractiveness_longterm) *
## epd1_m2: included * fertile + included * (menstruation + fertile) +
## epd1_m2: fertile_mean + (1 | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0 18 48445 48592 -24204 48409
## epd1_m2 20 48440 48604 -24200 48400 8.57 2 0.014 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre class="r"><code>summary(epd1_m2)</code></pre>
<pre><code>## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + (partner_attractiveness_shortterm + partner_attractiveness_longterm) * included * fertile + included * (menstruation + fertile) + fertile_mean + (1 | person)
## Data: diary
##
## REML criterion at convergence: 48499
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.283 -0.558 -0.148 0.403 7.989
##
## Random effects:
## Groups Name Variance Std.Dev.
## person (Intercept) 0.288 0.537
## Residual 0.320 0.566
## Number of obs: 26680, groups: person, 1054
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.8166 0.0465 1312.2463 39.10 < 2e-16 ***
## partner_attractiveness_longterm -0.1316 0.0315 1112.4575 -4.18 0.00003181 ***
## includedhorm_contra -0.0939 0.0378 1276.4053 -2.48 0.01310 *
## partner_attractiveness_shortterm -0.0467 0.0296 1103.1164 -1.58 0.11546
## fertile 0.1620 0.0351 25917.5850 4.62 0.00000394 ***
## menstruationpre -0.0895 0.0173 25915.8282 -5.18 0.00000022 ***
## menstruationyes -0.0698 0.0163 26012.3159 -4.28 0.00001909 ***
## fertile_mean -0.0834 0.2080 1440.4054 -0.40 0.68851
## partner_attractiveness_longterm:includedhorm_contra -0.0257 0.0404 1114.4358 -0.64 0.52537
## includedhorm_contra:partner_attractiveness_shortterm 0.0272 0.0388 1092.4121 0.70 0.48383
## partner_attractiveness_longterm:partner_attractiveness_shortterm 0.0354 0.0180 1043.0915 1.97 0.04886 *
## partner_attractiveness_shortterm:fertile -0.0845 0.0317 26015.9778 -2.67 0.00769 **
## partner_attractiveness_longterm:fertile 0.0774 0.0335 26103.5791 2.31 0.02080 *
## includedhorm_contra:fertile -0.1636 0.0446 26025.5272 -3.67 0.00024 ***
## includedhorm_contra:menstruationpre 0.0687 0.0222 25912.3930 3.09 0.00199 **
## includedhorm_contra:menstruationyes 0.0850 0.0214 25992.2371 3.98 0.00006998 ***
## includedhorm_contra:partner_attractiveness_shortterm:fertile 0.0537 0.0407 25970.3975 1.32 0.18692
## partner_attractiveness_longterm:includedhorm_contra:fertile -0.0116 0.0433 26061.4705 -0.27 0.78885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>Here, we see that we have different conclusions about the ST:fertile interaction depending on how we look at it. Why might this be? One potential reason why an effect might be different from zero but not explain more variance, is that our model currenty ignores that the extra-pair variable has a lowest possible value and is not normally distributed.</p>
<pre class="r"><code>qplot(diary$extra_pair)</code></pre>
<p><img src="3_fertility_robustness2_files/figure-html/unnamed-chunk-7-1.png" width="288" /></p>
<p>Let’s plot this to see the shape of the effect. The predictors are z-scored, but especially LT attractiveness is not normally distributed and ranges from -4 to 0.8. For ST attractiveness, the range is -3.3 to 1.6 so we take slightly odd cutpoints that map most of the range of the data.</p>
<pre class="r"><code>effs = allEffects(epd1_m)
effs = data.frame(effs$`partner_attractiveness_longterm:included:partner_attractiveness_shortterm:fertile`) %>%
filter(round(partner_attractiveness_longterm,1) %in% c(-3, -0.4,0.8),round(partner_attractiveness_shortterm,1) %in% c(-2,0.4,2))
ggplot(effs, aes(fertile, fit, ymin = lower, ymax = upper, color = included)) +
facet_grid(partner_attractiveness_shortterm ~ partner_attractiveness_longterm) +
geom_smooth(stat='identity') +
scale_color_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
scale_fill_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
ggtitle("Moderation", "top-to-bottom: short-term,\nleft-to-right: long-term attractiveness of the partner")+
ylab("Extra-pair desire and behaviour")</code></pre>
<p><img src="3_fertility_robustness2_files/figure-html/unnamed-chunk-8-1.png" width="768" /></p>
<p>Now, in our robustness checks, we omitted random slopes for menstruation and fertile effects.</p>
<p>I believe this leads us to overestimate the certainty in results, but didn’t worry about it, because I was short on computational time (given that we estimated thousands of mixed models) and thought the models already showed sufficiently clearly that estimates were uncertain and evidence did not favour moderators. However, Steve Gangestad disagreed with this somewhat, thinking the Wald t-test based p value of .006 above is evidence that we disregarded or dismissed. We thought we tested the right theoretical prediction using the model comparison above, and think this is reasonable, but for the sake of the argument, what happens to the NHST p value when we allow for random slopes?</p>
<pre class="r"><code>epd1_m0_r = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)
epd1_m_r = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)
epd1_m_r2 = lmer(extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + (partner_attractiveness_shortterm + partner_attractiveness_longterm) * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)</code></pre>
<ol style="list-style-type: lower-alpha">
<li>random slopes make the model fit better</li>
</ol>
<pre class="r"><code>anova(epd1_m, epd1_m_r)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m: partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m: included * fertile + included * (menstruation + fertile) +
## epd1_m: fertile_mean + (1 | person)
## epd1_m_r: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m_r: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m_r: partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m_r: included * fertile + included * (menstruation + fertile) +
## epd1_m_r: fertile_mean + (1 + menstruation + fertile | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m 23 48446 48634 -24200 48400
## epd1_m_r 32 48124 48386 -24030 48060 339 9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<ol start="2" style="list-style-type: lower-alpha">
<li>there is no evidence according to the LR test that we explain more variance with the ST:HC:fertile interaction after including random slopes</li>
</ol>
<pre class="r"><code>anova(epd1_m0_r, epd1_m_r)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0_r: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m0_r: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m0_r: partner_attractiveness_longterm * included * fertile + included *
## epd1_m0_r: (menstruation + fertile) + fertile_mean + (1 + menstruation +
## epd1_m0_r: fertile | person)
## epd1_m_r: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m_r: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m_r: partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m_r: included * fertile + included * (menstruation + fertile) +
## epd1_m_r: fertile_mean + (1 + menstruation + fertile | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0_r 27 48119 48341 -24033 48065
## epd1_m_r 32 48124 48386 -24030 48060 5.16 5 0.4</code></pre>
<ol start="3" style="list-style-type: lower-alpha">
<li>the p value for the ST:fertile interaction is no longer significant according to our .01 threshold</li>
</ol>
<pre class="r"><code>summary(epd1_m_r)</code></pre>
<pre><code>## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + (1 + menstruation + fertile | person)
## Data: diary
##
## REML criterion at convergence: 48169
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.751 -0.550 -0.137 0.390 7.902
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## person (Intercept) 0.3102 0.557
## menstruationpre 0.0439 0.210 -0.33
## menstruationyes 0.0673 0.259 -0.26 0.77
## fertile 0.2907 0.539 -0.17 0.44 0.57
## Residual 0.3028 0.550
## Number of obs: 26680, groups: person, 1054
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.82834 0.04783 1331.73679 38.23 < 2e-16 ***
## partner_attractiveness_longterm -0.13794 0.03212 1034.46492 -4.29 0.0000191 ***
## includedhorm_contra -0.10183 0.04063 1047.89725 -2.51 0.01236 *
## partner_attractiveness_shortterm -0.04417 0.02954 1028.62952 -1.50 0.13516
## fertile 0.16187 0.04788 860.21505 3.38 0.00076 ***
## menstruationpre -0.09721 0.02073 767.38149 -4.69 0.0000033 ***
## menstruationyes -0.07489 0.02166 759.38646 -3.46 0.00057 ***
## fertile_mean -0.12263 0.21177 1465.22044 -0.58 0.56261
## partner_attractiveness_longterm:includedhorm_contra -0.02050 0.04095 1035.80899 -0.50 0.61674
## includedhorm_contra:partner_attractiveness_shortterm 0.02597 0.03853 1019.57604 0.67 0.50045
## partner_attractiveness_longterm:partner_attractiveness_shortterm 0.03627 0.02358 1024.78110 1.54 0.12438
## partner_attractiveness_shortterm:fertile -0.08443 0.04217 790.95278 -2.00 0.04560 *
## partner_attractiveness_longterm:fertile 0.07618 0.04566 866.86778 1.67 0.09556 .
## includedhorm_contra:fertile -0.16461 0.06066 833.93321 -2.71 0.00679 **
## includedhorm_contra:menstruationpre 0.07747 0.02662 764.07430 2.91 0.00371 **
## includedhorm_contra:menstruationyes 0.09261 0.02815 783.40463 3.29 0.00105 **
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm 0.00954 0.03683 1032.96251 0.26 0.79572
## partner_attractiveness_longterm:partner_attractiveness_shortterm:fertile -0.02614 0.03424 867.52608 -0.76 0.44546
## includedhorm_contra:partner_attractiveness_shortterm:fertile 0.05515 0.05409 785.86542 1.02 0.30830
## partner_attractiveness_longterm:includedhorm_contra:fertile -0.00920 0.05839 841.98512 -0.16 0.87488
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm:fertile 0.03313 0.05296 847.11183 0.63 0.53176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<ol start="4" style="list-style-type: lower-alpha">
<li>same for the model not including ST:LT:fertile and ST:LT:HC:fertile</li>
</ol>
<pre class="r"><code>anova(epd1_m0_r, epd1_m_r2)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0_r: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m0_r: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m0_r: partner_attractiveness_longterm * included * fertile + included *
## epd1_m0_r: (menstruation + fertile) + fertile_mean + (1 + menstruation +
## epd1_m0_r: fertile | person)
## epd1_m_r2: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm *
## epd1_m_r2: included + partner_attractiveness_longterm * partner_attractiveness_shortterm +
## epd1_m_r2: (partner_attractiveness_shortterm + partner_attractiveness_longterm) *
## epd1_m_r2: included * fertile + included * (menstruation + fertile) +
## epd1_m_r2: fertile_mean + (1 + menstruation + fertile | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0_r 27 48119 48341 -24033 48065
## epd1_m_r2 29 48119 48357 -24031 48061 4.41 2 0.11</code></pre>
<pre class="r"><code>summary(epd1_m_r2)</code></pre>
<pre><code>## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: extra_pair ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + (partner_attractiveness_shortterm + partner_attractiveness_longterm) * included * fertile + included * (menstruation + fertile) + fertile_mean + (1 + menstruation + fertile | person)
## Data: diary
##
## REML criterion at convergence: 48155
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.750 -0.550 -0.137 0.389 7.904
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## person (Intercept) 0.3098 0.557
## menstruationpre 0.0440 0.210 -0.33
## menstruationyes 0.0672 0.259 -0.26 0.77
## fertile 0.2897 0.538 -0.17 0.44 0.57
## Residual 0.3028 0.550
## Number of obs: 26680, groups: person, 1054
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.8278 0.0474 1329.0098 38.55 < 2e-16 ***
## partner_attractiveness_longterm -0.1368 0.0313 1043.1936 -4.37 0.0000137 ***
## includedhorm_contra -0.0992 0.0391 1034.3098 -2.54 0.01136 *
## partner_attractiveness_shortterm -0.0436 0.0294 1030.1859 -1.49 0.13765
## fertile 0.1521 0.0461 831.2153 3.30 0.00100 **
## menstruationpre -0.0973 0.0207 767.4569 -4.69 0.0000032 ***
## menstruationyes -0.0750 0.0216 759.1252 -3.46 0.00056 ***
## fertile_mean -0.1258 0.2117 1466.5677 -0.59 0.55224
## partner_attractiveness_longterm:includedhorm_contra -0.0219 0.0401 1043.9424 -0.55 0.58494
## includedhorm_contra:partner_attractiveness_shortterm 0.0257 0.0384 1022.4230 0.67 0.50444
## partner_attractiveness_longterm:partner_attractiveness_shortterm 0.0387 0.0179 1038.5746 2.17 0.03049 *
## partner_attractiveness_shortterm:fertile -0.0800 0.0417 800.5668 -1.92 0.05552 .
## partner_attractiveness_longterm:fertile 0.0873 0.0432 874.0042 2.02 0.04375 *
## includedhorm_contra:fertile -0.1534 0.0587 811.6509 -2.61 0.00911 **
## includedhorm_contra:menstruationpre 0.0775 0.0266 764.1456 2.91 0.00369 **
## includedhorm_contra:menstruationyes 0.0927 0.0281 783.2011 3.29 0.00103 **
## includedhorm_contra:partner_attractiveness_shortterm:fertile 0.0512 0.0536 788.6214 0.96 0.33976
## partner_attractiveness_longterm:includedhorm_contra:fertile -0.0205 0.0564 845.5383 -0.36 0.71625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
</div>
<div id="model-comparison-extra-pair-desire" class="section level2">
<h2>Model comparison extra-pair desire</h2>
<p>Steve Gangestad also suggested that the moderation test should have focused on the extra-pair desire subscale, rather than the more heavily aggregated one based on EP desire and behaviour. Certainly, arguments can be made for both approaches.</p>
<pre class="r"><code>epd1 = lmer(extra_pair_desire ~ included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)
epd1_m0_old = lmer(extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)
epd1_m0 = lmer(extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm * included + partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)
epd1_m = lmer(extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)
epd1_m2 = lmer(extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + (partner_attractiveness_shortterm + partner_attractiveness_longterm) * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 | person), data = diary)
summary(epd1_m)</code></pre>
<pre><code>## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + (1 | person)
## Data: diary
##
## REML criterion at convergence: 52316
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.848 -0.493 -0.130 0.332 7.433
##
## Random effects:
## Groups Name Variance Std.Dev.
## person (Intercept) 0.376 0.613
## Residual 0.368 0.606
## Number of obs: 26680, groups: person, 1054
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.79098 0.05296 1277.04150 33.82 < 2e-16 ***
## partner_attractiveness_longterm -0.16115 0.03667 1096.82418 -4.39 0.0000121624 ***
## includedhorm_contra -0.18140 0.04451 1217.54640 -4.08 0.0000489058 ***
## partner_attractiveness_shortterm -0.09460 0.03379 1085.79702 -2.80 0.00521 **
## fertile 0.22580 0.03916 25878.17600 5.77 0.0000000082 ***
## menstruationpre -0.11161 0.01854 25878.74227 -6.02 0.0000000018 ***
## menstruationyes -0.07626 0.01749 25971.35112 -4.36 0.0000130554 ***
## fertile_mean -0.28968 0.23449 1395.47166 -1.24 0.21692
## partner_attractiveness_longterm:includedhorm_contra -0.01246 0.04676 1092.83008 -0.27 0.78985
## includedhorm_contra:partner_attractiveness_shortterm 0.02878 0.04410 1075.37170 0.65 0.51404
## partner_attractiveness_longterm:partner_attractiveness_shortterm 0.05312 0.02699 1095.57679 1.97 0.04935 *
## partner_attractiveness_shortterm:fertile -0.08793 0.03433 25974.13541 -2.56 0.01044 *
## partner_attractiveness_longterm:fertile 0.10315 0.03789 26065.54307 2.72 0.00649 **
## includedhorm_contra:fertile -0.21557 0.04938 25973.84928 -4.37 0.0000127466 ***
## includedhorm_contra:menstruationpre 0.07490 0.02379 25875.17192 3.15 0.00165 **
## includedhorm_contra:menstruationyes 0.08843 0.02292 25952.21221 3.86 0.00011 ***
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm 0.01625 0.04215 1093.20520 0.39 0.69995
## partner_attractiveness_longterm:partner_attractiveness_shortterm:fertile 0.02909 0.02805 26004.91919 1.04 0.29981
## includedhorm_contra:partner_attractiveness_shortterm:fertile 0.01481 0.04401 25935.65044 0.34 0.73645
## partner_attractiveness_longterm:includedhorm_contra:fertile -0.04588 0.04804 26024.27143 -0.96 0.33954
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm:fertile 0.00863 0.04342 25949.02856 0.20 0.84242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre class="r"><code>anova(epd1_m0_old, epd1_m)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0_old: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m0_old: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m0_old: partner_attractiveness_shortterm + included * (menstruation +
## epd1_m0_old: fertile) + fertile_mean + (1 | person)
## epd1_m: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m: partner_attractiveness_shortterm + partner_attractiveness_shortterm *
## epd1_m: partner_attractiveness_longterm * included * fertile + included *
## epd1_m: (menstruation + fertile) + fertile_mean + (1 | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0_old 16 52259 52390 -26113 52227
## epd1_m 23 52252 52440 -26103 52206 20.9 7 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre class="r"><code>anova(epd1_m0, epd1_m)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0: extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m0: included + partner_attractiveness_longterm * included * fertile +
## epd1_m0: included * (menstruation + fertile) + fertile_mean + (1 |
## epd1_m0: person)
## epd1_m: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m: partner_attractiveness_shortterm + partner_attractiveness_shortterm *
## epd1_m: partner_attractiveness_longterm * included * fertile + included *
## epd1_m: (menstruation + fertile) + fertile_mean + (1 | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0 19 52260 52416 -26111 52222
## epd1_m 23 52252 52440 -26103 52206 16.3 4 0.0027 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre class="r"><code>anova(epd1_m0, epd1_m2)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0: extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m0: included + partner_attractiveness_longterm * included * fertile +
## epd1_m0: included * (menstruation + fertile) + fertile_mean + (1 |
## epd1_m0: person)
## epd1_m2: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m2: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m2: partner_attractiveness_shortterm + (partner_attractiveness_shortterm +
## epd1_m2: partner_attractiveness_longterm) * included * fertile + included *
## epd1_m2: (menstruation + fertile) + fertile_mean + (1 | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0 19 52260 52416 -26111 52222
## epd1_m2 20 52249 52412 -26104 52209 13.7 1 0.00021 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre class="r"><code>effs = allEffects(epd1_m)
effs = data.frame(effs$`partner_attractiveness_longterm:included:partner_attractiveness_shortterm:fertile`) %>%
filter(round(partner_attractiveness_longterm,1) %in% c(-3, -0.4,0.8),round(partner_attractiveness_shortterm,1) %in% c(-2,0.4,2))
ggplot(effs, aes(fertile, fit, ymin = lower, ymax = upper, color = included)) +
facet_grid(partner_attractiveness_shortterm ~ partner_attractiveness_longterm) +
geom_smooth(stat='identity') +
scale_color_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
scale_fill_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
ggtitle("Moderation", "top-to-bottom: short-term,\nleft-to-right: long-term attractiveness of the partner")+
ylab("Extra-pair desire")</code></pre>
<p><img src="3_fertility_robustness2_files/figure-html/rob-14-1.png" width="768" /></p>
<pre class="r"><code>epd1_r = lmer(extra_pair_desire ~ included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)
epd1_m0_r = lmer(extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm * included + partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)
epd1_m_r = lmer(extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)
epd1_m_r2 = lmer(extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + (partner_attractiveness_shortterm + partner_attractiveness_longterm) * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary)
anova(epd1_m0_r, epd1_m_r)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m0_r: extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm *
## epd1_m0_r: included + partner_attractiveness_longterm * included * fertile +
## epd1_m0_r: included * (menstruation + fertile) + fertile_mean + (1 +
## epd1_m0_r: menstruation + fertile | person)
## epd1_m_r: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m_r: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m_r: partner_attractiveness_shortterm + partner_attractiveness_shortterm *
## epd1_m_r: partner_attractiveness_longterm * included * fertile + included *
## epd1_m_r: (menstruation + fertile) + fertile_mean + (1 + menstruation +
## epd1_m_r: fertile | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m0_r 28 51822 52051 -25883 51766
## epd1_m_r 32 51823 52085 -25879 51759 7.32 4 0.12</code></pre>
<pre class="r"><code>anova(epd1_m, epd1_m_r)</code></pre>
<pre><code>## Data: diary
## Models:
## epd1_m: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m: partner_attractiveness_shortterm + partner_attractiveness_shortterm *
## epd1_m: partner_attractiveness_longterm * included * fertile + included *
## epd1_m: (menstruation + fertile) + fertile_mean + (1 | person)
## epd1_m_r: extra_pair_desire ~ partner_attractiveness_longterm * included +
## epd1_m_r: partner_attractiveness_shortterm * included + partner_attractiveness_longterm *
## epd1_m_r: partner_attractiveness_shortterm + partner_attractiveness_shortterm *
## epd1_m_r: partner_attractiveness_longterm * included * fertile + included *
## epd1_m_r: (menstruation + fertile) + fertile_mean + (1 + menstruation +
## epd1_m_r: fertile | person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## epd1_m 23 52252 52440 -26103 52206
## epd1_m_r 32 51823 52085 -25879 51759 447 9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre class="r"><code>effs = allEffects(epd1_m_r)
effs = data.frame(effs$`partner_attractiveness_longterm:included:partner_attractiveness_shortterm:fertile`) %>%
filter(round(partner_attractiveness_longterm,1) %in% c(-3, -0.4,0.8),round(partner_attractiveness_shortterm,1) %in% c(-2,0.4,2))
ggplot(effs, aes(fertile, fit, ymin = lower, ymax = upper, color = included)) +
facet_grid(partner_attractiveness_shortterm ~ partner_attractiveness_longterm) +
geom_smooth(stat='identity') +
scale_color_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
scale_fill_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
ggtitle("Moderation", "top-to-bottom: short-term,\nleft-to-right: long-term attractiveness of the partner")+
ylab("Extra-pair desire")</code></pre>
<p><img src="3_fertility_robustness2_files/figure-html/rob-14-2.png" width="768" /></p>
<p>Here, we see that this association might turn significant according to the likelihood ratio test, depending on how you set things up (not with random slopes). But this is now really a post-hoc analysis, so for this reason and others spelled out below, I also used a Bayesian approach with LOO below.</p>
</div>
<div id="model-comparison-with-approximative-leave-one-out-cross-validation-in-brms" class="section level2">
<h2>Model comparison with approximative leave-one-out cross-validation in brms</h2>
<p>Here, I fit the same model as above, but as a Bayesian multilevel regression in brms with uninformative priors. I can then compare the models using approximative leave-one-out cross-validation (Vehtari, Gelman, Gabry, Yao, 2018). <a href="http://rpubs.com/rubenarslan/ordinal_floor">Previously</a>, I’ve found that doing so allows me to recover true moderators but exclude moderators that appear significant only because they move the outcome variable closer to the floor (and hence make it impossible for another predictor to explain variance that doesn’t exist). The more proper model here would be to model the data with one row per subject-item combination and use an ordinal model for the outcome. However, this is computationally very expensive and my simulations showed that the comparison using loo also yields the same answer under these circumstances. Also, modeling it per item would mean to model additional sources of uncertainty to be able to generalize to new items. Although I favour this now, it would make the model less comparable to the other models I ran, and I don’t want to get into the question of whether generalizing to new items is what we want with SG.</p>
<pre class="r"><code>library(brms)
epd1_m_r_b0 <- brm(extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm * included + partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary, future = T, iter = 2000, file = "epd1_m_r_b0")
epd1_m_r_b <- brm(extra_pair_desire ~ partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + ( 1 + menstruation + fertile | person), data = diary, future = T, iter = 2000, file = "epd1_m_r_b")</code></pre>
<p>LOO does not favour the model with the ST:HC:fertile interaction.</p>
<pre class="r"><code>comp <- loo(epd1_m_r_b0, epd1_m_r_b)
comp</code></pre>
<pre><code>## LOOIC SE
## epd1_m_r_b0 49359.64 411.8
## epd1_m_r_b 49362.29 411.8
## epd1_m_r_b0 - epd1_m_r_b -2.65 6.5</code></pre>
<p>The 95% credible interval also includes zero for the central interaction term, but, again, the estimates for the interaction term could be affected by the floor effect, which I plan to model in more detail in the future.</p>
<pre class="r"><code>summary(epd1_m_r_b)</code></pre>
<pre><code>## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: extra_pair_desire ~ partner_attractiveness_longterm * included + partner_attractiveness_shortterm * included + partner_attractiveness_longterm * partner_attractiveness_shortterm + partner_attractiveness_shortterm * partner_attractiveness_longterm * included * fertile + included * (menstruation + fertile) + fertile_mean + (1 + menstruation + fertile | person)
## Data: diary (Number of observations: 26680)
## Samples: 4 chains, each with iter = 1500; warmup = 750; thin = 1;
## total post-warmup samples = 3000
##
## Group-Level Effects:
## ~person (Number of levels: 1054)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0.64 0.02 0.61 0.67 676 1.01
## sd(menstruationpre) 0.25 0.02 0.22 0.28 841 1.01
## sd(menstruationyes) 0.29 0.02 0.26 0.32 1288 1.00
## sd(fertile) 0.63 0.03 0.56 0.70 838 1.00
## cor(Intercept,menstruationpre) -0.42 0.05 -0.52 -0.32 1672 1.00
## cor(Intercept,menstruationyes) -0.29 0.05 -0.38 -0.19 1498 1.00
## cor(menstruationpre,menstruationyes) 0.82 0.05 0.71 0.92 237 1.02
## cor(Intercept,fertile) -0.18 0.05 -0.28 -0.07 1819 1.00
## cor(menstruationpre,fertile) 0.45 0.07 0.31 0.57 310 1.01
## cor(menstruationyes,fertile) 0.61 0.06 0.49 0.71 218 1.01
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 1.80 0.05 1.70 1.91 612 1.00
## partner_attractiveness_longterm -0.17 0.03 -0.23 -0.10 517 1.00
## includedhorm_contra -0.19 0.05 -0.28 -0.10 250 1.01
## partner_attractiveness_shortterm -0.09 0.03 -0.15 -0.03 485 1.00
## fertile 0.21 0.05 0.11 0.31 1470 1.00
## menstruationpre -0.12 0.02 -0.17 -0.08 1254 1.00
## menstruationyes -0.08 0.02 -0.13 -0.04 1272 1.00
## fertile_mean -0.33 0.24 -0.80 0.13 867 1.00
## partner_attractiveness_longterm:includedhorm_contra -0.01 0.04 -0.09 0.08 430 1.01
## includedhorm_contra:partner_attractiveness_shortterm 0.03 0.04 -0.06 0.11 457 1.00
## partner_attractiveness_longterm:partner_attractiveness_shortterm 0.06 0.03 0.00 0.11 512 1.01
## partner_attractiveness_shortterm:fertile -0.09 0.05 -0.18 0.01 2005 1.00
## partner_attractiveness_longterm:fertile 0.11 0.05 0.01 0.21 1738 1.00
## includedhorm_contra:fertile -0.20 0.07 -0.33 -0.07 1438 1.00
## includedhorm_contra:menstruationpre 0.09 0.03 0.03 0.15 1195 1.00
## includedhorm_contra:menstruationyes 0.09 0.03 0.03 0.16 1331 1.00
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm 0.02 0.04 -0.06 0.10 524 1.01
## partner_attractiveness_longterm:partner_attractiveness_shortterm:fertile 0.02 0.04 -0.06 0.09 2102 1.00
## includedhorm_contra:partner_attractiveness_shortterm:fertile 0.02 0.06 -0.10 0.13 2074 1.00
## partner_attractiveness_longterm:includedhorm_contra:fertile -0.05 0.06 -0.18 0.07 1699 1.00
## partner_attractiveness_longterm:includedhorm_contra:partner_attractiveness_shortterm:fertile 0.02 0.06 -0.10 0.13 2243 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.59 0.00 0.58 0.59 3000 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).</code></pre>
<p>Finally, what does a plot of the marginal effects show?</p>
<pre class="r"><code>newd <- expand.grid(included = c("cycling", "horm_contra"), fertile = unique(diary$fertile), partner_attractiveness_shortterm = c(-2,0.4,2),
partner_attractiveness_longterm = c(-3, -0.4, 0.8), menstruation = "no", fertile_mean = 0) %>% na.omit()
newd[, c("Estimate", "Std.Error", "lower", "upper")] <- fitted(epd1_m_r_b, newd, re_formula = NA)
ggplot(newd, aes(fertile, Estimate, ymin = lower, ymax = upper, color = included)) +
facet_grid(partner_attractiveness_shortterm ~ partner_attractiveness_longterm) +
geom_smooth(stat='identity') +
scale_color_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
scale_fill_manual(values = c("cycling" = 'red', 'horm_contra' = 'black'), guide = F) +
ggtitle("Moderation", "top-to-bottom: short-term,\nleft-to-right: long-term attractiveness of the partner")+
ylab("Extra-pair desire")</code></pre>
<p><img src="3_fertility_robustness2_files/figure-html/unnamed-chunk-16-1.png" width="768" /></p>
</div>
<div id="summary" class="section level2">
<h2>Summary</h2>
<p>In summary, we think we were consistent in saying that our data don’t provide evidence for the predicted moderation effect (as we understood the prediction). However, they also do not rule out moderator effects of a smaller size than were previously reported. If this wasn’t clear, we’re happy to have corrected this misinterpretation. We wrote</p>
<blockquote>
<p>Although further tests should be conducted, the good genes ovulatory shift hypothesis could be wrong, given that we could not replicate previously reported moderators.</p>
</blockquote>
<p>We really meant this, not just as a paean to uncertainty, but we actually acted on it. We’ve collected data for a second study to test this and other theories again, with more data and better measures. We’re also involved in a larger project which aims to collect data across multiple labs in a Registered Report framework. These discussions may help us frame the predictions and models in this new project better. We may also improve measurement of our moderators, the ones in this study were based on those used in previous studies and maybe suboptimal in some ways.</p>
</div>
</div>
</div>
<script>
// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
$('tr.header').parent('thead').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
bootstrapStylePandocTables();
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>