We construct an updated and extended compilation of growth rate data consisting of 34 points and including corrections for model dependence. In order to maximize the independence of the datapoints we also construct a subsample of this compilation (`Gold' growth dataset) which consists of 18 datapoints. We test the consistency of this dataset with the best fit Planck15/$\Lambda$CDM parameters in the context of General Relativity (GR) using the evolution equation for the growth factor $\delta(a)$ with a $w$CDM background. We find tension at the $\sim 3 \sigma$ level between the best fit parameters $w$ (the dark energy equation of state), $\Omega_{0m}$ (the matter density parameter) and $\sigma_8$ (the matter power spectrum normalization on scales $8h^{-1}$Mpc) and the corresponding Planck15/$\Lambda$CDM parameters. We show that the tension disappears if we allow for evolution of the effective Newton's constant, parametrized as $G_{eff}(a)/G_N = 1 + g_a(1-a)^n-g_a(1-a)^{2n}$ with $n\ge2$ where $g_a$, $n$ are parameters of the model, $a$ is the scale factor and $z = 1/a-1$ is the redshift. This parametrization satisfies three criteria: a. $G_{eff} > 0$, b. Consistency with Big Bang Nucleosynthesis ($G_{eff}(a\ll 1)/G_N=1$), c. Consistency with solar system tests ($G_{eff}(a=1)/G_N=1$ and $G_{eff}'(a=1)/G_N=0$). We show that the best fit form of $G_{eff}(z)$ obtained from the growth data corresponds to weakening gravity at recent redshifts (decreasing function of $z$) and we demonstrate that this behavior is not consistent with any scalar-tensor Lagrangian with a real scalar field. Finally, we use MGCAMB to find the best fit $G_{eff}(z)$ obtained from the Planck CMB power spectrum on large angular scales and show that it is a mildly increasing function of $z$, in $3\sigma$ tension with the corresponding decreasing best fit $G_{eff}(z)$ obtained from the growth data.

Dark energy equation of state $w(z)$ parametrizations with two parameters and given monotonicity are generically either convex or concave functions. This makes them suitable for fitting either freezing or thawing quintessence models but not both simultaneously. Fitting a dataset based on a freezing model with an unsuitable (concave when increasing) $w(z)$ parametrization (like CPL) can lead to significant misleading features like crossing of the phantom divide line, incorrect $w(z=0)$, incorrect slope \etc that are not present in the underlying cosmological model. To demonstrate this fact we generate scattered cosmological data both at the level of $w(z)$ and the luminosity distance $D_L(z)$ based on either thawing or freezing quintessence models and fit them using parametrizations of convex and of concave type. We then compare statistically significant features of the best fit $w(z)$ with actual features of the underlying model. We thus verify that the use of unsuitable parametrizations can lead to misleading conclusions. In order to avoid these problems it is important to either use both convex and concave parametrizations and select the one with the best $\chi^2$ or use principal component analysis thus splitting the redshift range into independent bins. In the latter case however, significant information about the slope of $w(z)$ at high redshifts is lost. Finally, we propose a new family of parametrizations (nCPL) $w(z)=w_0+w_a (\frac{z}{1+z})^n$ which generalizes the CPL and interpolates between thawing and freezing parametrizations as the parameter $n$ increases to values larger than 1.