Analyzing data from paleoclimate archives such as tree rings or lake sediments offers the opportunity of inferring information on past climate variability. Often, such data sets are univariate and a proper reconstruction of the systems higher-dimensional phase space can be crucial for further analyses. In this study, we systematically compare the methods of time delay embedding and differential embedding for phase space reconstruction. Differential embedding relates the systems higher-dimensional coordinates to the derivatives of the measured time series. For implementation, this requires robust and efficient algorithms to estimate derivatives from noisy and possibly non-uniformly sampled data. For this purpose, we consider several approaches: (i) central differences adapted to irregular sampling, (ii) a generalized version of discrete Legendre coordinates and (iii) the concept of Moving Taylor Bayesian Regression. We evaluate the performance of differential and time delay embedding by studying two paradigmatic model systems - the Lorenz and the Rossler system. More precisely, we compare geometric properties of the reconstructed attractors to those of the original attractors by applying recurrence network analysis. Finally, we demonstrate the potential and the limitations of using the different phase space reconstruction methods in combination with windowed recurrence network analysis for inferring information about past climate variability. This is done by analyzing two well-studied paleoclimate data sets from Ecuador and Mexico. We find that studying the robustness of the results when varying the analysis parameters is an unavoidable step in order to make well-grounded statements on climate variability and to judge whether a data set is suitable for this kind of analysis.