I have finished my coalescent simulation analysis for all my populations, and now want to plot the results using ggplot, but I don't know how to parameterise it so that it matches the expectations for such a plot type. The above plot is kind of the expectation, but when I use the same parameters in ggplot, I produce the bottom one.
It's the more "jagged" lines that trail off the plot that I don't know how to create. The labels and axis are easy to modify.
A sample of the code I used for the top is here:
plot(CA_L_fixed_rho2_8cores$left_time_boundary/mu*gen, (1/CA_L_fixed_rho2_8cores$lambda)/(2*mu), log="x",ylim=c(0,100000), type="n", xlab="Years Ago", ylab="Effective Population Size")
lines(CA_L_fixed_rho2_8cores$left_time_boundary/mu*gen, (1/CA_L_fixed_rho2_8cores$lambda)/(2*mu), type="s", col=as.character(Pop_Col_Matrix[grep("CA_L", Pop_Col_Matrix$Population),2]))
lines(CA_R_fixed_rho2_8cores$left_time_boundary/mu*gen, (1/CA_R_fixed_rho2_8cores$lambda)/(2*mu), type="s", col=as.character(Pop_Col_Matrix[grep("CA_R", Pop_Col_Matrix$Population),2]))
And the ggplot is here:
ggplot(output, aes(x = left_time_boundary/mu*gen, y = (1/lambda_00)/(2*mu), linetype = Ecotype, colour = Population)) +
geom_line() +
ylim(0,100000)+
scale_x_log10(limits = c(100,10000)) +
theme_classic()
EDIT** Added 1 population of example data.
time_index left_time_boundary right_time_boundary lambda_00 Population
0 0.00000e+00 1.21183e-06 1.58059e+00 CA_L
1 1.21183e-06 2.45515e-06 7.37336e+01 CA_L
2 2.45515e-06 3.73162e-06 6.25560e+01 CA_L
3 3.73162e-06 5.04307e-06 2.34338e+02 CA_L
4 5.04307e-06 6.39146e-06 7.09668e+02 CA_L
5 6.39146e-06 7.77894e-06 1.30547e+03 CA_L
6 7.77894e-06 9.20785e-06 1.78073e+03 CA_L
7 9.20785e-06 1.06807e-05 2.02606e+03 CA_L
8 1.06807e-05 1.22004e-05 2.07207e+03 CA_L
9 1.22004e-05 1.37699e-05 2.06643e+03 CA_L
10 1.37699e-05 1.53926e-05 2.57822e+03 CA_L
11 1.53926e-05 1.70722e-05 2.57822e+03 CA_L
12 1.70722e-05 1.88129e-05 2.77558e+03 CA_L
13 1.88129e-05 2.06194e-05 2.77558e+03 CA_L
14 2.06194e-05 2.24967e-05 2.68765e+03 CA_L
15 2.24967e-05 2.44506e-05 2.68765e+03 CA_L
16 2.44506e-05 2.64877e-05 2.68153e+03 CA_L
17 2.64877e-05 2.86154e-05 2.68153e+03 CA_L
18 2.86154e-05 3.08421e-05 2.54312e+03 CA_L
19 3.08421e-05 3.31774e-05 2.54312e+03 CA_L
20 3.31774e-05 3.56325e-05 2.43018e+03 CA_L
21 3.56325e-05 3.82205e-05 2.43018e+03 CA_L
22 3.82205e-05 4.09563e-05 2.29070e+03 CA_L
23 4.09563e-05 4.38581e-05 2.29070e+03 CA_L
24 4.38581e-05 4.69472e-05 2.14761e+03 CA_L
25 4.69472e-05 5.02496e-05 2.14761e+03 CA_L
26 5.02496e-05 5.37967e-05 2.00975e+03 CA_L
27 5.37967e-05 5.76280e-05 2.00975e+03 CA_L
28 5.76280e-05 6.17928e-05 1.87895e+03 CA_L
29 6.17928e-05 6.63548e-05 1.87895e+03 CA_L
30 6.63548e-05 7.13978e-05 1.74989e+03 CA_L
31 7.13978e-05 7.70355e-05 1.74989e+03 CA_L
32 7.70355e-05 8.34270e-05 1.62090e+03 CA_L
33 8.34270e-05 9.08054e-05 1.62090e+03 CA_L
34 9.08054e-05 9.95322e-05 1.48383e+03 CA_L
35 9.95322e-05 1.10213e-04 1.48383e+03 CA_L
36 1.10213e-04 1.23983e-04 1.33870e+03 CA_L
37 1.23983e-04 1.43390e-04 1.33870e+03 CA_L
38 1.43390e-04 1.76568e-04 1.11796e+03 CA_L
39 1.76568e-04 Inf 1.11796e+03 CA_L
Those two plots look like they were generated from different data sets. So it's a little bit hard to figure out what you want exactly. But I think what you're after is essentially cropping so it's "zoomed" in. In that case us
coord_cartesian(ylim=c(), xlim=c())
to set the upper and lower bounds.https://stackoverflow.com/questions/25685185/limit-ggplot2-axes-without-removing-data-outside-limits-zoom
Hi Amar, it's the same data. I've added one population worth of data so you can see for yourself. It's not zooming that I'm after, it's to make the lines less "smoothed" and more abrupt when they go back in time as this is the norm for this type of plot in publications. Among other differences that are obvious from the plots.
Miles,
I'm attempting to do something similar but have two species for which I've generated psmc plots. I just want to plot them on the same axes to save some space in my figures and to allow direct comparison.
I don't typically use R but can if needed. My main question is this: how did you get the data structure you posted above?
How did you get this information from the psmc file? Did you get this information from the psmc file? If not, from where did it come? If I can get my psmc data into a format like this, I could do what I need to but I don't know how to get that step completed.
Any help would be appreciated.
You're likely to get more help if you post your own question. Piggybacking on another one makes it tough for other users to find your question.