Is it possible to make a scatterplot that fits this model?
mymodel <- glm(cbind(Success,Failure) ~ Temperature + Pesticide +
Temperature*Pesticide, data=mydata, family=binomial(link="logit"))
I have hatching success data from an experiment exposing an invertebrate to a pesticide at different temperatures. I want to make a scatterplot of the datapoints with a regression line from the glm (one line per pesticide treatment (three), in different colours. If i plot the raw data (0/1 response for each egg, instead of the proportion of an egg clutch that successfully hatched) I get the datapoints, but then i dont manage to add the regression, and also not sure if it is correct.
I use R, and prefer ggplot2.
This didn´t work:
library(ggplot2)
myplot <- ggplot(mydata, aes(y=cbind(Success, Failure), x=Temperature, group=Pesticide, color=Pesticide)) +
theme_light()+
scale_color_manual(values=c("#009E73","#F0E442","#FF9999")) +
labs(x="Temperature (°C)", y="Hatching succeess")+
geom_point() +
stat_smooth(here I will put my model)
Thank you for any tips!
output of subset using dput()
structure(list(Temperature = c(10, 10, 10, 10, 10, 10, 10, 10,
10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
10, 10, 10, 10, 10, 10, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 15, 15, 15, 15, 15, 15, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 25, 25, 25, 25, 25, 25, 25,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
25, 25, 25, 25), Pesticide = c("CONTROL", "CONTROL", "CONTROL",
"CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL",
"CONTROL", "LOW", "LOW", "LOW", "LOW", "LOW", "LOW", "LOW", "LOW",
"LOW", "LOW", "LOW", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "CONTROL", "CONTROL", "LOW",
"LOW", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "CONTROL", "CONTROL", "LOW", "HIGH", "CONTROL", "HIGH",
"LOW", "HIGH", "LOW", "CONTROL", "HIGH", "CONTROL", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "LOW", "CONTROL", "LOW", "LOW",
"LOW", "LOW", "LOW", "LOW", "CONTROL", "CONTROL", "CONTROL",
"CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL",
"CONTROL", "CONTROL", "LOW", "CONTROL", "CONTROL", "CONTROL",
"CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL",
"CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL", "CONTROL",
"CONTROL", "LOW", "HIGH", "CONTROL", "LOW", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH",
"HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "HIGH", "LOW", "LOW",
"LOW", "LOW", "LOW", "HIGH", "LOW", "LOW", "HIGH", "HIGH", "CONTROL",
"CONTROL", "HIGH", "LOW", "LOW", "HIGH", "HIGH", "CONTROL", "HIGH",
"CONTROL", "CONTROL", "HIGH", "HIGH", "HIGH", "HIGH", "CONTROL",
"LOW", "LOW", "CONTROL", "LOW", "CONTROL", "CONTROL", "LOW",
"LOW", "CONTROL", "CONTROL", "LOW", "CONTROL", "LOW", "CONTROL",
"CONTROL", "LOW", "CONTROL", "LOW", "LOW", "CONTROL", "LOW",
"CONTROL", "CONTROL", "LOW", "LOW", "CONTROL", "LOW", "CONTROL",
"CONTROL", "LOW", "LOW", "CONTROL", "HIGH", "LOW", "CONTROL",
"LOW", "LOW", "LOW", "LOW", "CONTROL", "LOW", "CONTROL", "CONTROL",
"LOW", "LOW"), Batch_ID = c("156_c", "156_d", "153_g", "153_h",
"153_e", "153_f", "156_m", "156_o", "156_n", "153_d", "166_h",
"166_i", "165_b", "165_a", "168_d", "166_d", "166_c", "166_b",
"168_a", "164_a", "166_a", "178_i", "177_h", "177_i", "180_g",
"177_e", "179_k", "179_l", "177_d_2", "178_k", "179_j", "179_i",
"177_f", "177_g", "179_h", "177_e_2", "179_m", "177_d", "178_l",
"178_m", "180_o", "179_f", "177_c2", "179_g", "177_c", "180_n",
"180_m", "177_b", "180_i", "179_d", "179_e", "179_d2", "179_e2",
"180_k", "180_l", "179_b", "179_c", "173_b", "180_j", "173_a",
"177_a", "178_j", "180_h", "179_a", "180_c", "180_d", "180_e",
"180_f", "178_h", "178_f", "178_d", "180_b", "178_g", "178_a",
"178_c", "178_e", "178_b", "180_a", "192_g2", "189_o", "201_g",
"204_h", "216_a", "215_n", "215_o", "215_m", "215_k", "215_l",
"216_h", "216_j", "192_p", "192_s", "201_u", "215_f2", "192_q",
"215_f3", "201_t", "216_g", "201_s", "192_o", "216_f2", "192_f2",
"216_k", "216_l", "216_e2", "213_o", "213_n", "216_i", "204_n",
"190_b", "202_c", "202_d", "202_a", "202_b", "203_b", "203_d",
"189_i", "189_j", "189_k", "189_h", "191_a", "191_b", "192_a",
"192_b", "187_a", "185_b", "186_b", "203_a", "186_c", "185_c",
"185_d", "185_e", "188_c", "188_d", "185_a", "188_b", "190_a",
"189_b", "189_d", "189_a", "188_a", "187_b", "186_a", "192_c",
"239_j", "249_o", "228_g", "238_b", "250_h", "252_o", "249_k",
"249_m", "249_l", "249_n", "249_j", "250_o2", "250_o", "251_o",
"249_i", "250_n", "250_m", "249_e", "249_g", "250_g", "251_h2",
"238_i", "238_j", "238_k", "238_l", "238_m", "252_n", "240_n",
"238_h", "251_b2", "251_c", "227_k", "227_m", "252_h", "237_i",
"240_g", "251_d", "251_e", "226_a", "238_c", "227_f", "227_c",
"240_a", "246_a", "246_b", "247_a", "223_a", "233_a", "233_b",
"261_i_2", "273_j2", "261_g", "260_a", "276_a", "273_k", "261_i",
"261_j", "273_m", "261_h", "273_l", "264_d", "264_e", "274_d",
"261_h", "275_h2", "273_j", "263_e", "273_h", "262_f", "262_g",
"273_d", "273_e", "261_d", "275_h", "261_b", "263_d", "274_b",
"273_c", "262_c", "287_a2", "275_e2", "263_b", "275_b2", "275_c2",
"275_d2", "275_f2", "262_d", "274_c", "264_a", "261_a", "275_g",
"275_f"), No_eggs = c(45, 15, 14, 24, 15, 33, 21, 18, 11, 13,
11, 21, 15, 20, 14, 30, 23, 15, 28, 6, 10, 14, 17, 18, 74, 22,
12, 10, 19, 42, 20, 30, 12, 34, 35, 32, 6, 14, 23, 33, 35, 14,
17, 36, 25, 18, 7, 30, 35, 23, 14, 15, 29, 21, 13, 15, 15, 10,
17, 9, 38, 24, 4, 43, 16, 25, 5, 9, 18, 24, 6, 26, 9, 30, 3,
6, 9, 19, 14, 34, 11, 21, 4, 19, 12, 34, 19, 32, 20, 19, 16,
25, 27, 14, 21, 19, 9, 22, 34, 31, 10, 28, 19, 22, 19, 12, 37,
14, 23, 16, 24, 7, 12, 11, 15, 14, 5, 11, 18, 15, 3, 8, 13, 5,
8, 17, 14, 6, 13, 13, 10, 13, 12, 14, 10, 13, 9, 13, 12, 17,
13, 11, 11, 4, 10, 22, 20, 10, 16, 14, 19, 20, 8, 15, 23, 14,
18, 6, 14, 18, 10, 18, 24, 3, 10, 19, 9, 6, 14, 8, 8, 20, 10,
9, 7, 12, 14, 12, 18, 14, 11, 9, 7, 9, 7, 7, 11, 11, 11, 6, 11,
17, 5, 11, 24, 16, 5, 8, 8, 16, 7, 7, 3, 11, 17, 8, 5, 3, 4,
11, 5, 10, 11, 12, 6, 10, 3, 3, 6, 3, 12, 18, 4, 1, 10, 6, 17,
5, 6, 10, 5, 3, 1, 2, 3, 2), Success = c(44, 15, 2, 6, 13, 17,
8, 6, 11, 0, 6, 18, 15, 20, 14, 25, 17, 10, 27, 3, 7, 14, 8,
14, 64, 22, 11, 9, 18, 0, 19, 30, 6, 31, 33, 12, 0, 14, 0, 0,
34, 13, 1, 29, 22, 16, 6, 30, 34, 22, 14, 15, 29, 0, 0, 15, 15,
9, 17, 8, 38, 21, 3, 34, 12, 25, 4, 9, 0, 13, 5, 25, 8, 24, 3,
6, 9, 18, 14, 31, 11, 11, 3, 13, 7, 34, 13, 11, 6, 12, 2, 19,
24, 0, 14, 17, 8, 6, 34, 31, 8, 26, 19, 22, 2, 12, 36, 13, 11,
15, 17, 5, 12, 11, 13, 14, 5, 11, 18, 14, 3, 8, 13, 5, 8, 17,
14, 6, 13, 12, 9, 13, 12, 14, 9, 12, 8, 13, 12, 17, 13, 11, 10,
0, 9, 19, 20, 1, 3, 12, 2, 0, 0, 15, 7, 13, 10, 3, 14, 18, 7,
17, 2, 1, 0, 19, 8, 6, 14, 8, 8, 8, 0, 1, 2, 1, 0, 0, 0, 1, 11,
6, 7, 9, 7, 7, 11, 7, 11, 5, 11, 16, 4, 0, 1, 0, 0, 8, 0, 0,
0, 0, 2, 2, 4, 0, 0, 2, 0, 8, 1, 0, 0, 0, 0, 7, 3, 0, 3, 2, 2,
13, 3, 0, 8, 5, 3, 5, 6, 9, 2, 0, 1, 2, 0, 0), Failure = c(1,
0, 12, 18, 2, 16, 13, 12, 0, 13, 5, 3, 0, 0, 0, 5, 6, 5, 1, 3,
3, 0, 9, 4, 10, 0, 1, 1, 1, 42, 1, 0, 6, 3, 2, 20, 6, 0, 23,
33, 1, 1, 16, 7, 3, 2, 1, 0, 1, 1, 0, 0, 0, 21, 13, 0, 0, 1,
0, 1, 0, 3, 1, 9, 4, 0, 1, 0, 18, 11, 1, 1, 1, 6, 0, 0, 0, 1,
0, 3, 0, 10, 1, 6, 5, 0, 6, 21, 14, 7, 14, 6, 3, 14, 7, 2, 1,
16, 0, 0, 2, 2, 0, 0, 17, 0, 1, 1, 12, 1, 7, 2, 0, 0, 2, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0,
0, 0, 0, 0, 1, 4, 1, 3, 0, 9, 13, 2, 17, 20, 8, 0, 16, 1, 8,
3, 0, 0, 3, 1, 22, 2, 10, 0, 1, 0, 0, 0, 0, 12, 10, 8, 5, 11,
14, 12, 18, 13, 0, 3, 0, 0, 0, 0, 0, 4, 0, 1, 0, 1, 1, 11, 23,
16, 5, 0, 8, 16, 7, 7, 1, 9, 13, 8, 5, 1, 4, 3, 4, 10, 11, 12,
6, 3, 0, 3, 3, 1, 10, 5, 1, 1, 2, 1, 14, 0, 0, 1, 3, 3, 0, 0,
3, 2)), row.names = c(NA, -236L), class = "data.frame")
I think you want something like this?
library(ggplot2)
## create sensible non-default ordering for the factor
mydata <- transform(mydata,
Pesticide = factor(Pesticide,
levels =c("CONTROL", "LOW", "HIGH")))
ggplot(mydata,
aes(Temperature, Success/No_eggs,
colour = Pesticide)) +
geom_point(alpha = 0.6, aes(size = No_eggs)) +
geom_smooth(method = "glm",
aes(weight = No_eggs),
method.args = list(family = binomial),
formula = y~x ## default: avoids message
) +
theme_light()+
scale_color_manual(values=c("#009E73","#F0E442","#FF9999")) +
labs(x="Temperature (°C)", y="Hatching succeess")
Notes:
weights mapping in the geom_smooth (this corresponds to the weights argument in glm(), which you use for a binomial model when you use a proportion as the response — see ?glm