We consider the important problem of estimating parameter sensitivities for stochastic models of reaction networks that describe the dynamics as a continuous-time Markov process over a discrete lattice. These sensitivity values are useful for understanding network properties, validating their design and identifying the pivotal model parameters. Many methods for sensitivity estimation have been developed, but their computational feasibility suffers from the critical bottleneck of requiring time-consuming Monte Carlo simulations of the exact reaction dynamics. To circumvent this problem one needs to devise methods that speed up the computations while suffering acceptable and quantifiable loss of accuracy. We develop such a method by first deriving a novel integral representation of parameter sensitivity and then demonstrating that this integral may be approximated by any convergent tau-leap method. Our method is easy to implement, works with any tau-leap simulation scheme and its accuracy is proved to be similar to that of the underlying tau-leap scheme. We demonstrate the efficiency of our methods through numerical examples. We also compare our method with the tau-leap versions of certain finite-difference schemes that are commonly used for sensitivity estimations.