Für ein Projekt an Blattaustrieb an verschiedenen Baumarten möchte ich eine «Survival» /»Time-to-event» Analyse mit einem CoxPH Modell durchführen.
Ganz grob habe ich untersucht wie verschiedene Umeltfaktoren (Chilling_Hours, Forcing_GDH etc) und chemische Parameter (SA_ng_g etc.) sich auf die Blattaustriebsgeschwindigeit (Days_first_bud_st2) auswirken.
Dazu habe ich folgendes Modell verwendet:
fit_phyto <- coxph(
Surv(Days_first_bud_st2_NAis60, Status) ~ Species*IAA_ng_g+Species*SA_ng_g+Species*Chilling_Hours+Species*Forcing_GDH_O_B_4+frailty(Tree.ID),
data = dtaX
)
Nun möchte ich einzelne signifkante Interaktionen (z.B. Species*Chilling_Hours) einem «marginal effects plot» darstellen.
Untenstehen mal ein Plot in dem ich die Rohdaten «Chilling_Hours vs. Days_first_bud_st2 geplottet habe. Ich möchte in meinem marginal effects plot etwas ähnliches erreichen und «Chilling_Hours vs. «Probability of first budburst» i.e. »Wahrscheinlichkeit für Knospenaustrieb» plotten (oder ist es Probability of budburst not happening?)
Gerne würde ich das ganze mit ggplot plotten.
Kann mir da jemand weiterhelfen? Insbesonder mit dem frailty term klappt es nicht. Notfalls würde ich den aber einfach weglassen
Untenstehend der R-code und Daten für mein coxph Modell
Vielen herzlichen Dank für eure Hilfe.
# CODE MODELL
library(ggplot2)
library(ggsurvfit)
library(ggthemes)
library(dplyr)
library(tidyr)
library(RColorBrewer)
library(survival)
library(survminer)
library(coxme)
Code: Alles auswählen
fit_phyto <- coxph(
Surv(Days_first_bud_st2_NAis60, Status) ~ Species*IAA_ng_g+Species*SA_ng_g+Species*Chilling_Hours+Species*Forcing_GDH_O_B_4+frailty(Tree.ID),
data = dtaX
)
dtaX<-structure(list(Days_first_bud_st2_NAis60 = c(1L, 2L, 2L, 2L,
3L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L,
13L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 15L, 15L,
15L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L,
17L, 17L, 18L, 18L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L,
20L, 21L, 21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 22L, 22L, 22L,
22L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 24L, 24L, 26L,
27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 28L, 28L, 29L, 29L, 29L,
29L, 29L, 29L, 30L, 30L, 30L, 31L, 31L, 31L, 32L, 32L, 32L, 33L,
33L, 33L, 33L, 34L, 34L, 35L, 35L, 35L, 35L, 36L, 37L, 37L, 37L,
37L, 37L, 39L, 39L, 39L, 40L, 40L, 41L, 41L, 43L, 44L, 44L, 45L,
46L, 47L, 47L, 47L, 48L, 48L, 48L, 48L, 48L, 49L, 49L, 51L, 51L,
51L, 51L, 52L, 53L, 53L, 55L, 55L, 57L, 60L, 60L, 60L, 60L, 60L,
60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L,
60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L, 60L
), Status = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Species = structure(c(2L,
4L, 4L, 4L, 2L, 2L, 2L, 2L, 1L, 2L, 4L, 4L, 4L, 4L, 2L, 4L, 1L,
2L, 2L, 2L, 2L, 4L, 4L, 4L, 1L, 1L, 1L, 4L, 4L, 1L, 2L, 4L, 4L,
4L, 4L, 4L, 4L, 1L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 1L, 4L, 4L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 1L, 1L,
1L, 1L, 2L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 4L, 1L, 1L, 2L, 4L, 4L,
4L, 4L, 1L, 3L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 4L, 4L, 2L,
3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 3L, 4L, 1L, 1L, 1L, 2L, 3L, 4L,
4L, 4L, 1L, 1L, 1L, 2L, 3L, 3L, 3L, 3L, 4L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 3L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 2L,
2L, 1L, 1L, 4L, 1L, 3L, 4L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 1L,
2L, 3L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 4L, 3L, 3L, 2L, 3L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 2L, 2L,
2L, 3L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L), levels = c("Beech", "Hornbeam", "Maple",
"Oak"), class = "factor"), IAA_ng_g = c(14.565, 12.295, 15.935,
28.315, 13.775, 11.915, 11.505, 16.185, 7.055, 7.985, 25.835,
31.785, 19.325, 14.015, 6.255, 11.215, 11.575, 14.605, 26.355,
7.155, 10.415, 7.115, 52.175, 20.685, 12.305, 7.835, 6.865, 7.275,
14.205, 9.795, 10.805, 23.565, 50.485, 23.375, 44.455, 17.705,
11.405, 8.845, 7.605, 5.175, 9.305, 8.255, 32.345, 178.605, 18.785,
23.335, 195.435, 39.345, 22.545, 43.275, 30.335, 0.005, 65.625,
21.835, 7.305, 6.035, 5.065, 0.005, 57.995, 6.895, 35.765, 6.375,
29.175, 5.085, 10.015, 5.105, 0.005, 10.945, 15.235, 12.785,
22.055, 47.345, 50.885, 0.005, 5.885, 38.135, 6.225, 8.695, 7.635,
27.525, 16.025, 37.625, NA, 4.645, 8.625, 76.685, 21.755, 22.905,
19.905, 9.905, 0.005, 0.005, 0.005, 17.945, 117.545, 5.895, 5.345,
10.595, 27.435, 16.425, 15.815, 46.505, 5.165, 4.795, 6.965,
8.115, 18.155, 5.185, 6.335, 0.005, 5.685, 6.865, 44.005, 60.785,
94.505, 0.005, 0.005, 0.005, 6.205, 9.535, 6.525, 10.395, 13.915,
50.825, 5.585, 8.575, 6.955, 0.005, 0.005, 0.005, 4.865, 8.895,
7.365, 84.625, 428.865, 0.005, 5.045, 8.695, 0.005, 6.415, 0.005,
7.585, 11.705, 0.005, 6.795, 0.005, 9.275, 5.585, 36.655, 5.635,
5.955, 18.355, 10.595, 0.005, 8.465, 17.765, 5.625, 9.535, 0.005,
6.325, 0.005, 9.735, 5.105, 5.675, 0.005, 4.795, 8.835, 9.125,
0.005, 0.005, 24.845, 19.915, 7.015, 0.005, 8.105, 6.365, 0.005,
0.005, 0.005, 5.315, 0.005, 18.065, 5.905, 0.005, 6.485, 0.005,
8.685, 5.435, 8.285, 4.985, 5.185, 0.005, 6.555, 4.735, 5.855,
6.145, 0.005, 0.005, 11.585, 5.955, 8.375, 9.065, 20.305, 11.605,
6.445, 5.465, 0.005, 7.185, 4.715, 8.995, 9.865, 11.565, 5.645,
5.565, 8.435, 5.465, 4.975, 8.105, 11.145, 13.435, 15.425, 6.205,
5.555, 7.615, 6.915, NA, 13.475, 8.505, 18.875, 6.795, 12.445
), SA_ng_g = c(8127.205, 100610.505, 151579.605, 91827.105, 9712.605,
7373.605, 3211.005, 4800.905, 2253.405, 5583.105, 104505.805,
71635.605, 80039.805, 130087.905, 8362.905, 117421.505, 3000.005,
8000.905, 4247.405, 4777.105, 5078.505, 45277.105, 185819.205,
106966.805, 7634.705, 1761.205, 6250.905, 133932.405, 80941.005,
1828.405, 6703.005, 81327.705, 153563.505, 94199.005, 175187.305,
70723.405, 75926.205, 2250.105, 5983.105, 4066.405, 3247.405,
7600.905, 94489.005, 101121.905, 73277.905, 82610.505, 164115.105,
117505.505, 82479.805, 135294.205, 210496.005, 1721.605, 138562.605,
117924.005, 3819.305, 3898.105, 1441.405, 6553.005, 108407.305,
89243.805, 47331.005, 72610.505, 96024.105, 7456.805, 2341.105,
5573.505, 3896.705, 1993.105, 80470.105, 156783.405, 78344.405,
70186.205, 128095.805, 2248.305, 2002.705, 79103.705, 3378.505,
1851.905, 5483.305, 138705.705, 68095.705, 112580.205, NA, 2273.705,
1191.605, 232896.705, 153643.305, 103460.705, 139717.105, 195114.905,
2835.605, 6554.005, 7115.605, 1931.405, 95453.105, 95963.905,
7132.105, 746.105, 138767.605, 155923.405, 119726.505, 108449.505,
9500.005, 1480.405, 2847.005, 4303.705, 152748.405, 1386.305,
1419.005, 4456.605, 3749.305, 262.705, 71703.005, 182235.905,
72325.105, 7121.305, 6352.705, 1461.805, 8717.105, 354.205, 804.805,
902.905, 773.605, 70051.205, 1887.905, 1558.105, 1832.205, 1897.605,
7414.405, 3936.005, 2198.005, 7226.905, 1288.805, 111252.905,
109157.205, 1222.705, 1698.705, 1995.705, 1720.605, 2557.505,
2061.105, 1699.605, 391.505, 18784.505, 5133.005, 2616.805, 2028.105,
3390.205, 130018.005, 5768.905, 125.105, 123983.505, 2077.505,
3166.505, 1275.605, 556.105, 6134.305, 9362.605, 6372.505, 7740.005,
7756.505, 7749.805, 319.405, 2023.505, 6090.105, 1654.105, 8032.505,
1555.205, 5675.105, 4166.005, 46760.705, 624.105, 1122.505, 3206.805,
1093.605, 1687.905, 2832.605, 6013.105, 7616.405, 2446.905, 4204.405,
1428.705, 10532.105, 1417.305, 1519.805, 1106.205, 330.405, 248.505,
7360.305, 13061.105, 7313.705, 2024.505, 2434.305, 10862.905,
701.905, 1806.605, 2036.805, 4396.805, 1250.505, 125.705, 7217.405,
7383.005, 20606.005, 12168.305, 19871.505, 13866.005, 14668.305,
10054.105, 13998.405, 7083.105, 3926.705, 7686.505, 7772.005,
908.405, 499.505, 1304.705, 707.605, 1207.505, 1167.605, 1066.705,
890.705, 2075.705, 357.005, 202.805, 341.005, NA, 639.605, 830.305,
1776.405, 402.705, 1153.005), Chilling_Hours = c(1473L, 1811L,
150L, 1811L, 1473L, 1473L, 1473L, 1473L, 1811L, 1473L, 1811L,
1811L, 1811L, 1811L, 1473L, 1473L, 1811L, 1337L, 1337L, 1337L,
1337L, 1811L, 6L, 1667L, 1811L, 1811L, 1811L, 6L, 1667L, 1811L,
1096L, 1667L, 150L, 1667L, 6L, 150L, 1667L, 1811L, 1096L, 1337L,
1337L, 1337L, 1473L, 6L, 1337L, 1473L, 6L, 6L, 1473L, 6L, 6L,
1667L, 1667L, 1667L, 1096L, 1096L, 1096L, 1096L, 6L, 6L, 1473L,
1473L, 150L, 1667L, 1667L, 1667L, 1667L, 1096L, 541L, 1096L,
541L, 1096L, 1096L, 1473L, 1473L, 1337L, 1667L, 1667L, 672L,
1337L, 1337L, 150L, 1337L, 6L, 1811L, 6L, 541L, 541L, 541L, 6L,
1473L, 1473L, 1473L, 6L, 1096L, 1337L, 672L, 1811L, 150L, 672L,
150L, 1096L, 6L, 1337L, 6L, 1811L, 6L, 1337L, 1473L, 1473L, 541L,
1473L, 541L, 1096L, 541L, 1337L, 1337L, 1337L, 672L, 1667L, 1473L,
1667L, 1811L, 672L, 6L, 6L, 1337L, 1096L, 1096L, 150L, 672L,
541L, 1473L, 1096L, 672L, 1096L, 1337L, 6L, 6L, 6L, 1096L, 1096L,
1667L, 6L, 541L, 541L, 6L, 6L, 6L, 150L, 1337L, 672L, 6L, 1096L,
1337L, 1337L, 541L, 672L, 6L, 150L, 1096L, 150L, 672L, 672L,
672L, 672L, 672L, 672L, 541L, 672L, 672L, 1096L, 1096L, 541L,
1096L, 150L, 541L, 541L, 672L, 150L, 541L, 541L, 150L, 541L,
150L, 672L, 1096L, 1473L, 6L, 6L, 150L, 150L, 541L, 150L, 672L,
150L, 150L, 541L, 672L, 541L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 150L, 6L, 6L, 6L, 6L, 150L, 541L, 6L, 6L, 150L, 541L,
6L, 6L, 150L, 672L, 6L, 6L, 150L, 541L, 1337L, 1667L), Forcing_GDH_O_B_4 = c(1131.67828,
2993.231743, 0, 2993.231743, 1131.67828, 1131.67828, 1131.67828,
1131.67828, 2993.231743, 1131.67828, 2993.231743, 2993.231743,
2993.231743, 2993.231743, 1131.67828, 1131.67828, 2993.231743,
433.4186002, 433.4186002, 433.4186002, 433.4186002, 2993.231743,
0, 1495.425431, 2993.231743, 2993.231743, 2993.231743, 0, 1495.425431,
2993.231743, 396.2630564, 1495.425431, 0, 1495.425431, 0, 0,
1495.425431, 2993.231743, 396.2630564, 433.4186002, 433.4186002,
433.4186002, 1131.67828, 0, 433.4186002, 1131.67828, 0, 0, 1131.67828,
0, 0, 1495.425431, 1495.425431, 1495.425431, 396.2630564, 396.2630564,
396.2630564, 396.2630564, 0, 0, 1131.67828, 1131.67828, 0, 1495.425431,
1495.425431, 1495.425431, 1495.425431, 396.2630564, 0.963261379,
396.2630564, 0.963261379, 396.2630564, 396.2630564, 1131.67828,
1131.67828, 433.4186002, 1495.425431, 1495.425431, 8318.579175,
433.4186002, 433.4186002, 0, 433.4186002, 0, 2993.231743, 0,
0.963261379, 0.963261379, 0.963261379, 0, 1131.67828, 1131.67828,
1131.67828, 0, 396.2630564, 433.4186002, 8318.579175, 2993.231743,
0, 4.317929707, 0, 396.2630564, 0, 433.4186002, 0, 2993.231743,
0, 433.4186002, 1131.67828, 1131.67828, 8315.224506, 1131.67828,
0.963261379, 396.2630564, 0.963261379, 433.4186002, 433.4186002,
433.4186002, 8318.579175, 1495.425431, 1131.67828, 1495.425431,
2993.231743, 4.317929707, 0, 0, 433.4186002, 396.2630564, 396.2630564,
7918.401951, 8318.579175, 8315.224506, 1131.67828, 396.2630564,
4.317929707, 396.2630564, 433.4186002, 0, 0, 0, 396.2630564,
396.2630564, 1495.425431, 0, 8315.224506, 8315.224506, 0, 0,
0, 0, 433.4186002, 4.317929707, 0, 396.2630564, 433.4186002,
433.4186002, 8315.224506, 8318.579175, 0, 0, 396.2630564, 7918.401951,
8318.579175, 4.317929707, 4.317929707, 4.317929707, 8318.579175,
8318.579175, 0.963261379, 4.317929707, 4.317929707, 396.2630564,
396.2630564, 8315.224506, 396.2630564, 0, 0.963261379, 0.963261379,
4.317929707, 0, 0.963261379, 0.963261379, 7918.401951, 0.963261379,
0, 4.317929707, 396.2630564, 1131.67828, 1904.665701, 5341.045771,
7918.401951, 7918.401951, 8315.224506, 7918.401951, 8318.579175,
0, 0, 0.963261379, 8318.579175, 8315.224506, 1904.665701, 5341.045771,
5341.045771, 1904.665701, 1904.665701, 5341.045771, 1904.665701,
5341.045771, 1904.665701, 5341.045771, 7918.401951, 1904.665701,
5341.045771, 1904.665701, 5341.045771, 7918.401951, 8315.224506,
1904.665701, 5341.045771, 7918.401951, 8315.224506, 1904.665701,
5341.045771, 7918.401951, 8318.579175, 1904.665701, 5341.045771,
7918.401951, 8315.224506, 433.4186002, 1495.425431), Tree.ID = structure(c(8L,
30L, 31L, 31L, 9L, 10L, 13L, 15L, 1L, 12L, 26L, 25L, 33L, 32L,
14L, 31L, 4L, 8L, 9L, 10L, 15L, 34L, 27L, 30L, 2L, 3L, 5L, 25L,
31L, 6L, 8L, 25L, 29L, 33L, 27L, 30L, 34L, 7L, 9L, 12L, 13L,
14L, 34L, 25L, 31L, 30L, 28L, 28L, 25L, 31L, 30L, 1L, 26L, 32L,
10L, 11L, 12L, 14L, 32L, 30L, 32L, 33L, 27L, 2L, 3L, 5L, 7L,
13L, 30L, 31L, 27L, 29L, 30L, 1L, 3L, 25L, 4L, 6L, 10L, 30L,
32L, 28L, 27L, 5L, 22L, 29L, 31L, 25L, 29L, 31L, 4L, 5L, 2L,
4L, 25L, 33L, 11L, 23L, 32L, 31L, 25L, 32L, 2L, 3L, 7L, 17L,
32L, 1L, 6L, 7L, 9L, 17L, 32L, 28L, 28L, 2L, 5L, 6L, 14L, 17L,
22L, 22L, 24L, 32L, 3L, 6L, 4L, 1L, 5L, 10L, 12L, 14L, 23L, 27L,
28L, 3L, 7L, 1L, 3L, 4L, 4L, 6L, 23L, 5L, 10L, 11L, 6L, 7L, 29L,
5L, 17L, 27L, 1L, 7L, 22L, 24L, 8L, 9L, 2L, 2L, 2L, 11L, 17L,
1L, 2L, 3L, 8L, 13L, 4L, 7L, 30L, 17L, 20L, 12L, 22L, 1L, 1L,
2L, 5L, 4L, 5L, 6L, 9L, 3L, 6L, 6L, 21L, 24L, 9L, 10L, 8L, 12L,
13L, 14L, 22L, 3L, 7L, 7L, 20L, 21L, 8L, 8L, 9L, 10L, 11L, 11L,
12L, 12L, 13L, 13L, 13L, 14L, 14L, 17L, 17L, 17L, 17L, 20L, 20L,
20L, 20L, 21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 23L, 24L), levels = c("Beech_1",
"Beech_2", "Beech_3", "Beech_4", "Beech_5", "Beech_6", "Beech_7",
"Hornbeam_1", "Hornbeam_2", "Hornbeam_3", "Hornbeam_4", "Hornbeam_5",
"Hornbeam_6", "Hornbeam_7", "Hornbeam_8", "Maple_1", "Maple_2",
"Maple_3", "Maple_4", "Maple_5", "Maple_6", "Maple_7", "Maple_8",
"Maple_9", "Oak_1", "Oak_10", "Oak_2", "Oak_3", "Oak_4", "Oak_5",
"Oak_6", "Oak_7", "Oak_8", "Oak_9"), class = "factor")), row.names = c(NA,
-231L), class = "data.frame")