I am trying to compute p-values for the slope obtained using xarray.polyfit() by extracting the covariance matrix and using a t-test. However, I am getting 1s for all p-values, which seems incorrect. The same datasets give valid p-values when using scipy.stats.linregress() for specific points but slow to apply to the whole dataset. def compute_trends(data : xr.Dataset): trends = xr.Dataset() # data = data.sortby("time") n = data.dims["time"] # Number of time points # Perform linear regression using xarray.polyfit() ds = data.polyfit(dim="time", deg=1, cov=True) # Extract slope and covariance matrix slope = ds["polyfit_coefficients"].sel(degree=1) slope_variance = ds["polyfit_covariance"].sel(cov_i=0, cov_j=0) # Variance of the slope # Compute standard error, t-statistic, and p-values stderr = np.sqrt(slope_variance) t_stat = slope / stderr p_values = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n - 2)) # Store in dataset trends["slope"] = slope trends["p_value"] = xr.DataArray(p_values, coords=slope.coords, dims=slope.dims) return trends I suspect an issue with how I’m extracting the covariance matrix (polyfit_covariance) or computing stderr. Am I correctly extracting the slope variance from polyfit_covariance? Continue reading...
