One way I know requires the assumption that the function is nice enough, i.e., analytic on the interval that you have some values available. Then you can write the Taylor series expansion of that function around the point you are trying to approximate the derivative. Evaluate these expressions at given points, then try to find a linear combination that gives the best approximation of the first derivative.
More generally, let's assume we have values at points $x_1 < \ldots < x_n$ which are not very far from each other, and we want to approximate $f^{(k)}(x_0)$ where $x_0$ is inside the interval $(x_1, x_n)$. The Taylor expansion of $f$ around $x_0$ is
$f(x) = \sum_{i=0}^\infty f^{(i)}(x_0) \frac{(x - x_0)^i}{i!}.$
Evaluate at $x_1, \ldots, x_n$ to get $n$ equations:
$f(x_k) = \sum_{i=0}^\infty f^{(i)}(x_0)\frac{(x_k - x_0)^i}{i!}, \quad k = 1, 2, \ldots, n.$
Suppose $h = \max_k\left|x_k - x_0\right|$ (or we can use the rougher estimate $h = x_n - x_1$). The above equation can be written as
$f(x_k) = \sum_{i=0}^{n - 1} f^{(i)}(x_0)\frac{(x_k - x_0)^i}{i!} + O(h^n), \quad k = 1, 2, \ldots, n.$
Ignoring the $O(h^n)$, we get a system of linear equations:
$f(x_k) = \sum_{i=0}^{n - 1} \hat f^{(i)}(x_0)\frac{(x_k - x_0)^i}{i!}, \quad k = 1, 2, \ldots, n$
where $\hat f$ is an approximation of $f$. You can write this as a matrix-vector equation
$ \begin{pmatrix} 1 & \frac{x_1 - x_0}{1!} & \frac{(x_1 - x_0)^2}{2!} & \ldots & \frac{(x_1 - x_0)^{n-1}}{(n-1)!} \\ 1 & \frac{x_2 - x_0}{1!} & \frac{(x_2 - x_0)^2}{2!} & \ldots & \frac{(x_2 - x_0)^{n-1}}{(n-1)!} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \frac{x_n - x_0}{1!} & \frac{(x_n - x_0)^2}{2!} & \ldots & \frac{(x_n - x_0)^{n-1}}{(n-1)!} \\ \end{pmatrix} \begin{pmatrix} \hat f(x_0) \\ \hat f'(x_0) \\ \vdots \\ \hat f^{(n-1)}(x_0) \end{pmatrix} = \begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{pmatrix}. $
Inverting the matrix on the left gives you formulas for approximating $f^{(k)}$ that are at least order $n - k$ accurate.
I think $\hat f$ is actually the Lagrange polynomial that passes through all the given points if you impose $\hat f^{(k)} = 0$ for $k \ge n$. However, the matrix-vector equation is much easier to use (at least in my opinion).