Paper 27: Linkage between fecal androgen and glucocorticoid metabolites, spermaturia, body weight and onset of puberty in male African lions (Panthera leo)

References

Putman SB, Brown JL, Saffoe C, Franklin AD, Pukazhenthi BS (2019) Linkage between fecal androgen and glucocorticoid metabolites, spermaturia, body weight and onset of puberty in male African lions (Panthera leo). PLoS ONE 14(7):e0217986. https://doi.org/10.1371/journal.pone.0217986

Disclosure

This reproducibility project was conducted to the best of our ability, with careful attention to statistical methods and assumptions. The research team comprises four senior biostatisticians (three of whom are accredited), with 20 to 30 years of experience in statistical modelling and analysis of healthcare data. While statistical assumptions play a crucial role in analysis, their evaluation is inherently subjective, and contextual knowledge can influence judgements about the importance of assumption violations. Differences in interpretation may arise among statisticians and researchers, leading to reasonable disagreements about methodological choices.

Our approach aimed to reproduce published analyses as faithfully as possible, using the details provided in the original papers. We acknowledge that other statisticians may have differing success in reproducing results due to variations in data handling and implicit methodological choices not fully described in publications. However, we maintain that research articles should contain sufficient detail for any qualified statistician to reproduce the analyses independently.

Methods used in our reproducibility analyses

There were two parts to our study. First, 100 articles published in PLOS ONE were randomly selected from the health domain and sent for post-publication peer review by statisticians. Of these, 95 included linear regression analyses and were therefore assessed for reporting quality. The statisticians evaluated what was reported, including regression coefficients, 95% confidence intervals, and p-values, as well as whether model assumptions were described and how those assumptions were evaluated. This report provides a brief summary of the initial statistical review.

The second part of the study involved reproducing linear regression analyses for papers with available data to assess both computational and inferential reproducibility. All papers were initially assessed for data availability and the statistical software used. From those with accessible data, the first 20 papers (from the original random sample) were evaluated for computational reproducibility. Within each paper, individual linear regression models were identified and assigned a unique number. A maximum of three models per paper were selected for assessment. When more than three models were reported, priority was given to the final model or the primary models of interest as identified by the authors; any remaining models were selected at random.

To assess computational reproducibility, differences between the original and reproduced results were evaluated using absolute discrepancies and rounding error thresholds, tailored to the number of decimal places reported in each paper. Results for each reported statistic, e.g., regression coefficient, were categorised as Reproduced, Incorrect Rounding, or Not Reproduced, depending on how closely they matched the original values. Each paper was then classified as Reproduced, Mostly Reproduced, Partially Reproduced, or Not Reproduced. The mostly reproduced category included cases with minor rounding or typographical errors, whereas partially reproduced indicated substantial errors were observed, but some results were reproduced.

For models deemed at least partially computationally reproducible, inferential reproducibility was further assessed by examining whether statistical assumptions were met and by conducting sensitivity analyses, including bootstrapping where appropriate. We examined changes in standardized regression coefficients, which reflect the change in the outcome (in standard deviation units) for a one standard deviation increase in the predictor. Meaningful differences were defined as a relative change of 10% or more, or absolute differences of 0.1 (moderate) and 0.2 (substantial). When non-linear relationships were identified, inferential reproducibility was assessed by comparing model fit measures, including R², Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC). When the Gaussian distribution was not appropriate for the dependent variable, alternative distributions were considered, and model fit was evaluated using AIC and BIC.

Results from the reproduction of the Putman et al. (2019) paper are presented below. An overall summary of results is presented first, followed by model-specific results organised within tab panels. Within each panel, the Original results tab displays the linear regression outputs extracted from the published paper. The Reproduced results tab presents estimates derived from the authors’ shared data, along with a comprehensive assessment of linear regression assumptions. The Differences tab compares the original and reproduced models to assess computational reproducibility. Finally, the Sensitivity analysis tab evaluates inferential reproducibility by examining whether identified assumption violations meaningfully affected the results.

Summary from statistical review

This paper explored metabolites, body weight, and puberty onset in male African lions. The primary analysis in this paper is based on repeated measures of ANCOVA. However, the article goes between using repeated measures and averaging results in a confusing manner. Results are reported in the text rather than in tables; Degrees of Freedom was reported for repeated measures but not for linear regression. ANCOVA/ linear regression models examined associations between captivity and weight. No assumptions were checked, collinearity was not mentioned, and modelling was unclear.

Data availability and software used

The authors provided data in a working Excel file in the supporting information. The file was poorly structured for use in statistical software, requiring each analysis dataset to be identified and extracted manually. Importing the data was also difficult because the Excel sheets contained duplicate column names and inconsistent naming conventions. SAS was reported as the statistical software used for the analyses.

Regression sample

This paper used a mix of repeated-measures ANOVA and linear regression; however, it was often unclear which analysis corresponded to which results, making it difficult to determine the total number of regression models fitted. The authors reported using linear regression to examine weight gain over 30 months in captive and wild male lions. Based on the reporting, it appeared that four linear regression models were fitted: one for wild lions and three for captive male lions. These included one overall model and two additional models stratified by follow-up duration (up to 24 months and beyond 24 months).

Visualisations were provided for two of the models in the main paper. For these models, the regression equations and 𝑅² values were reported in the supporting information. For the remaining two captive models, only the slope estimates were reported, without accompanying intercepts or model fit statistics. To align with the level of reporting in the original paper, reproducibility was assessed for the overall wild-lion model and for the two captive-lion models for which slope estimates were reported.

Computational reproducibility results

This paper was not computationally reproducible, as only one of the three models could be reproduced. The regression model for wild male weight was reproduced. However, interpretation was complicated by inconsistencies in the coding of the time variable: months ranged from 0 to 30 in the dataset but needed to be recoded to 1 to 31 to obtain the same intercept results. The two models examining captive weight could not be reproduced; the timing was unclear, and the model beyond 24 months had missing data, with only three points for the regression.

Inferential reproducibility results

The model for wild lion weights was not inferentially reproduced due to concerns about linearity. A quadratic term for month was therefore considered. Relative to the linear specification, the quadratic model provided a modestly improved fit (AIC lower by 6.2; BIC lower by 4.9) and improved diagnostic behaviour, with clearer residual patterns and attenuation of departures from linearity and normality. The quadratic term captured mild concave-down curvature in the month–weight relationship rather than a strictly linear increase. On balance, the quadratic model was considered better specified and was retained as the final model.

Recommended Changes

• Provide data in recognised wide or long format with a unique animal identifier (ID).
• Provide tables in the Supporting Information that present all analyses conducted in the paper, including full model outputs such as regression coefficients and all variables used for adjustment.
• Consider linearity of lion weight with time.
• Evaluate the assumptions of the linear regression models by examining residuals, identifying influential outliers, and assessing multicollinearity among predictors. If any assumptions are violated, address them using appropriate methods.
• Include a data dictionary.

Model 1

Original results

Model results for Wild Male Weight (KG)

Term	B	SE	Lower	Upper	t	p-value
Intercept	9.5168
Month	3.8759
SE = Standard error; Lower = lower confidence interval; Upper = upper confidence interval.

Fit statistics for Wild Male Weight (KG)

R	R2	R2Adj	AIC	RMSE	F	DF1	DF2	p-value
	0.9451
R2 Adj = Adjusted R2; AIC = Akaike Information Criterion; RMSE = The Root Mean Squared Error; DF1 = Degrees of freedom for the model; DF2 = Degrees of freedom for the residuals.

ANOVA table for Wild Male Weight (KG)

Term	SS	DF	MS	F	p-value
Month
Residuals
SS = Sum of Squares; DF = Degrees of freedom; MS = Mean Square.

Reproduced results

Model results for Wild Male Weight (KG)

Term	B	SE	Lower	Upper	t	p-value
Intercept	9.517	3.558	2.174	16.859	2.675	0.0132
Month	3.876	0.191	3.482	4.269	20.332	<0.001
SE = Standard error; Lower = lower confidence interval; Upper = upper confidence interval.

Fit statistics for Wild Male Weight (KG)

R	R2	R2Adj	AIC	RMSE	F	DF1	DF2	p-value
0.972	0.945	0.943	192.805	8.789	413.390	1	24	<0.001
R2 Adj = Adjusted R2; AIC = Akaike Information Criterion; RMSE = The Root Mean Squared Error; DF1 = Degrees of freedom for the model; DF2 = Degrees of freedom for the residuals.

ANOVA table for Wild Male Weight (KG)

Term	SS	DF	MS	F	p-value
Month	34,591.164	1	34,591.164	413.390	<0.001
Residuals	2,008.245	24	83.677
SS = Sum of Squares; DF = Degrees of freedom; MS = Mean Square; Calculated using type III SS.

Visualisation of regression model

The blue line shows the best line of fit with shading representing 95% confidence intervals, while holding all other covariates constant. The dots show partial residuals, which reflect the observed data adjusted for all other predictors except the one being plotted.

Checking residuals plots for patterns

Blue line showing quadratic fit for residuals

Testing residuals for non linear relationships

Term	Statistic	p-value	Results
Month	−2.920	0.0077	Linearity violation
Tukey test	−2.920	0.0035	Linearity violation
Specification test for predictors using quadratic tests, for fitted values curvature is tested through Tukey's one-degree-of-freedom test for nonadditivity.

Checking univariate relationships with the dependent variable using scatterplots

Blue line shows linear relationship, red line indicates relationship inferred by GAM modelling

Linearity results

• Plots and tests show an unaccounted non-linear relationship may be present.

Testing for homoscedasticity

Statistic	p-value	Parameter	Method
0.420	0.5170	1	studentized Breusch-Pagan test

Homoscedasticity results

• The studentized Breusch-Pagan test supports homoscedasticity.
• There is no distinct funnelling pattern observed, supporting homoscedasticity of residuals.

Model descriptives including cook’s distance and leverage to understand outliers

Term	N	Mean	SD	Median	Min	Max	Skewness	Kurtosis
Wild Male Weight (KG)	26	71.978	38.262	75.081	6.947	126.900	−0.228	−1.419
Month	26	16.115	9.597	16.000	1.000	31.000	−0.013	−1.448
.fitted	26	71.978	37.197	71.531	13.393	129.668	−0.013	−1.448
.resid	26	−0.000	8.963	−1.520	−11.788	23.318	1.326	0.958
.leverage	26	0.077	0.032	0.070	0.039	0.138	0.455	−1.260
.sigma	26	9.144	0.374	9.281	7.914	9.344	−2.379	4.645
.cooksd	26	0.033	0.045	0.013	0.000	0.149	1.457	0.961
.std.resid	26	−0.009	1.011	−0.176	−1.358	2.604	1.294	0.894
dfb.1_	26	0.001	0.148	0.003	−0.310	0.360	−0.011	0.173
dfb.Mnth	26	0.001	0.164	−0.001	−0.361	0.258	−0.460	−0.258
dffit	26	−0.029	0.274	−0.058	−0.461	0.631	0.891	0.364
cov.r	26	1.094	0.163	1.143	0.585	1.223	−2.125	3.584
* categorical variable

Cooks threshold

Cook’s distance measures the overall change in fit, if the ith observation is removed. Potential influential observations are identified by \(\text{Cook's Distance}_i > \frac{4}{n}\), where n is the number of observations. In practice a threshold of 0.5 to 1 is often used to identify influential observations.

DFFIT threshold

DFFIT measures how many standard deviations the fitted values will change when the ith observation is removed. Potential influential observations \(\left| \text{DFFITS}_i \right| > \frac{2\sqrt{p}}{\sqrt{n}}\) where p is the number of predictors (including the intercept) and n is the number of observations. In practice this can result in a large number of points identified, a practical cut-off of 1 was used to flag observations with meaningful impact.

DFBETA threshold

DFBETAS quantify the influence of the ith observation on the jth regression coefficient as the change in that coefficient when the observation is omitted, expressed in units of the coefficient’s estimated standard error. There is a DFBETA for each parameter in the model. Potential influential observations \(|\text{DFBETA}_{ij}| > \frac{2}{\sqrt{n}}\), where n is the number of observations. In larger datasets, this threshold can flag a high number of observations with only minor influence on the model. A practical cut-off of 1 was used to flag observations with meaningful impact.

Influence plot

Observations with high leverage (horizontal) and large residuals (vertical, typically at ±2 or ±3 studentized residuals) are concerning, as they may disproportionately influence the model. This combination is reflected by large bubbles with high Cook’s distance indicated by darker shadings of blue.

COVRATIO plot

COVRATIO measures the overall change in the precision (covariance matrix) of the estimated regression coefficients when the ith observation is removed. Values close to 1 indicate little influence on the model’s precision. Values below 1 suggest that an observation inflates the variances and reduces precision, resulting in wider confidence intervals, whereas values above 1 suggest deflated variances and narrower confidence intervals. A commonly cited guideline is \(\left|\mathrm{COVRATIO}_i - 1\right| > \frac{3p}{n}\), where p is the number of parameters and n is the number of observations. A practical cut-off between 0.9 to 1.1 was used to flag observations with meaningful impact on precision, although there is no agreed universal alternative cut-off.

Observations of interest identified by the influence plot

ID	StudRes	Leverage	CookD	dfb.1_	dfb.Mnth	dffit	cov.r
1	−0.752	0.138	0.046	−0.300	0.255	−0.300	1.203
31	−1.071	0.135	0.089	0.195	−0.357	−0.423	1.142
20	2.854	0.045	0.148	0.085	0.236	0.620	0.622
19	3.010	0.042	0.149	0.145	0.185	0.631	0.585
StudRes = studentized residuals; CookD = Cook's Distance a combined measure of leverage and influence. DFBETAS (dfb.*) measures how much a specific regression coefficient changes (in standard errors) when an observation is removed; DFFITS measures how much the fitted (predicted) value for an observation changes (in standard deviations) when that observation is removed; cov.r = coefficient covariance ratio which measures how much the overall variance (precision) of the coefficients changes when that observation is removed.

Results for outliers and influential points

Two observations had studentized residuals near or greater than 3. Both obsevations had low leverage and small Cook’s distance, with DFBETAS and DFFITS within conventional ranges. The COVRATIO indicated observations that may affect confidence intervals widths.

Checking for normality of the residuals using a Q–Q plot

Normality of residuals using Shapiro-Wilk and Kolmogorov-Smirnov tests

Statistic	p-value	Method
0.235	0.0957	Exact one-sample Kolmogorov-Smirnov test

Statistic	p-value	Method
0.844	0.0011	Shapiro-Wilk normality test

Normality results

• The Kolmogorov-Smirnov supports residuals being normally distributed.
• The Shapiro-Wilk normality test indicates residuals may not be normally distributed.
• QQ-plot indicates the residuals are not normally distributed.

Assessing independence with the Durbin–Watson test for autocorrelation

AutoCorrelation	Statistic	p-value
0.279	1.380	0.0820

Independence results

• The Durbin–Watson test suggests there are no auto-correlation issues.
• The study design was not independent. Although the model used aggregated data, measurements were taken repeatedly from the same population over time; therefore, residuals should be carefully assessed for autocorrelation and lag effects.

Assumption conclusions

Normality tests and inspection of the Q–Q plot suggested that the residuals may not be normally distributed. Outlier diagnostics indicated that point estimates were unlikely to be substantially influenced by individual observations; however, the width of the confidence intervals could be affected and warrants further investigation. Residual plots also suggested potential non-linearity in weight over time, which should be assessed by fitting a quadratic term. Although no statistically significant autocorrelation was detected, this assumption should be interpreted cautiously, as measurements were aggregated at the population level but derived from repeated assessments of the same population over time.

Differences

Forest plot showing original and reproduced coefficients and 95% confidence intervals for Wild Male Weight (KG)

Change in regression coefficients

term	O_B	R_B	Change.B	reproduce.B
Intercept	9.5168	9.5168	0.0000	Reproduced
Month	3.8759	3.8759	0.0000	Reproduced
O_B = original B; R_B = reproduced B; Change.B = change in R_B - O_B; Reproduce.B = B reproduced.

Change in R²

O_R2	R_R2	Change.R2	Reproduce.R2
0.945	0.9451	0.0000	Reproduced
O_R2 = original R2; R_R2 = reproduced R2; Change.R2 = change in R_R2 - O_R2

Conclusion computational reproducibility

This model was computationally reproducible, with all reported statistics that were assessed being reproducible.

Sensitivity analysis

Methods

The model was successfully reproduced; however, residual diagnostics suggested potential non-linearity. To assess the impact of this departure from linearity, non-linear terms were introduced for month and the resulting model fit (AIC and BIC) were compared with the original specification. In these models, all variables were standardized. The non-linear predictor, which was centered and scaled prior to creation of squared and cubic terms. Visualisation was performed using scatterplots on the original scale to illustrate the fitted relationships from the linear and quadratic.

Linearity conclusion

A quadratic term for Month was retained. Compared with the linear specification, the quadratic model provided a moderately better fit (AIC lower by 6.2; BIC lower by 4.9) and improved diagnostic behaviour, with clearer residual patterns and attenuation of the previously observed departures from linearity and normality. With month standardized, a one–standard deviation increase in month was associated with a 0.92 standard deviation increase in weight at the mean month, with a small negative quadratic term (β = −0.14) indicating mild concave-down curvature. This specification is also conceptually plausible, as mild curvature over time would be expected. On balance, the quadratic model was considered a better fit and was retained. Therefore, this model was not considered inferentially reproducible.

Characteristic	Linear			quadratic
	Beta	95% CI1	p-value	Beta	95% CI1	p-value
(Intercept)	-0.0117	-0.1085, 0.0851	0.805	0.1333	0.0002, 0.2664	0.050
z_Month	0.9210	0.8275, 1.0145	<0.001	0.9226	0.8408, 1.0043	<0.001
z_Month2				-0.1354	-0.2312, -0.0395	0.008
1 CI = Confidence Interval

Model_type	R2adj	sigma	df	logLik	AIC	BIC	deviance	df.residual	nobs
Linear	0.943	0.239	1.000	1.353	3.293	7.067	1.372	24	26
Quadratic	0.956	0.209	2.000	5.454	−2.908	2.125	1.001	23	26

Although predictors were standardized to facilitate comparison of regression coefficients, the fitted relationship is plotted on the original scale to aid interpretation of curvature.

Model 2

Original results

Model results for Captive male Weight (KG) 24 months

Term	B	SE	Lower	Upper	t	p-value
Intercept
Month	7.27
SE = Standard error; Lower = lower confidence interval; Upper = upper confidence interval.

Fit statistics for Captive male Weight (KG) 24 months

R	R2	R2Adj	AIC	RMSE	F	DF1	DF2	p-value
R2 Adj = Adjusted R2; AIC = Akaike Information Criterion; RMSE = The Root Mean Squared Error; DF1 = Degrees of freedom for the model; DF2 = Degrees of freedom for the residuals.

ANOVA table for Captive male Weight (KG) 24 months

Term	SS	DF	MS	F	p-value
Month
Residuals
SS = Sum of Squares; DF = Degrees of freedom; MS = Mean Square.

Reproduced results

Model results Captive male Weight (KG) 24 months

Term	B	SE	Lower	Upper	t	p-value
Intercept	−13.076	3.765	−20.905	−5.246	−3.473	0.0023
Month	7.644	0.275	7.073	8.215	27.837	<0.001
SE = Standard error; Lower = lower confidence interval; Upper = upper confidence interval.

Model fit for Captive male Weight (KG) 24 months

R	R2	R2Adj	AIC	RMSE	F	DF1	DF2	p-value
0.987	0.974	0.972	168.876	8.347	774.917	1	21	<0.001
R2 Adj = Adjusted R2; AIC = Akaike Information Criterion; RMSE = The Root Mean Squared Error; DF1 = Degrees of freedom for the model; DF2 = Degrees of freedom for the residuals.

ANOVA table for Captive male Weight (KG) 24 months

Term	SS	DF	MS	F	p-value
Month	59,126.499	1	59,126.499	774.917	<0.001
Residuals	1,602.309	21	76.300
SS = Sum of Squares; DF = Degrees of freedom; MS = Mean Square; Calculated using type III SS.

Visualisation of regression model

The blue line shows the best line of fit with shading representing 95% confidence intervals, while holding all other covariates constant. The dots show partial residuals, which reflect the observed data adjusted for all other predictors except the one being plotted.

Differences

Forest plot showing Original and Reproduced coefficients and 95% confidence intervals for Captive male Weight (KG) 24 months

Change in regression coefficients

term	O_B	R_B	Change.B	reproduce.B
Intercept		−13.0757
Month	7.27	7.6437	0.3737	Not Reproduced
O_B = original B; R_B = reproduced B; Change.B = change in R_B - O_B; Reproduce.B = B reproduced.

Conclusion computational reproducibility

This model was not computationally reproducible.

Sensitivity analysis

As this model was not computationally reproducible, inferential reproducibility was not considered, since the original analyses could not be reproduced and therefore, statistical assumptions could not be meaningfully compared or interpreted.

Model 3

Original results

Model results for Captive male Weight (KG) 25-30 months

Term	B	SE	Lower	Upper	t	p-value
Intercept
Month	2.66
SE = Standard error; Lower = lower confidence interval; Upper = upper confidence interval.

Fit Statistics Captive male Weight (KG) 25-30 months

R	R2	R2Adj	AIC	RMSE	F	DF1	DF2	p-value
R2 Adj = Adjusted R2; AIC = Akaike Information Criterion; RMSE = The Root Mean Squared Error; DF1 = Degrees of freedom for the model; DF2 = Degrees of freedom for the residuals.

ANOVA Table for Captive male Weight (KG) 25-30 months

Term	SS	DF	MS	F	p-value
Month
Residuals
SS = Sum of Squares; DF = Degrees of freedom; MS = Mean Square.

Reproduced results

Model results for Captive male Weight (KG) 25-30 months

Term	B	SE	Lower	Upper	t	p-value
Intercept	47.461	23.394	−249.785	344.706	2.029	0.2915
Month	3.553	0.843	−7.162	14.267	4.213	0.1484
SE = Standard error; Lower = lower confidence interval; Upper = upper confidence interval.

Fit statistics for Captive male Weight (KG) 25-30 months

R	R2	R2Adj	AIC	RMSE	F	DF1	DF2	p-value
0.973	0.947	0.893	17.812	1.733	17.750	1	1	0.1484
R2 Adj = Adjusted R2; AIC = Akaike Information Criterion; RMSE = The Root Mean Squared Error; DF1 = Degrees of freedom for the model; DF2 = Degrees of freedom for the residuals.

ANOVA Table for Captive male Weight (KG) 25-30 months

Term	SS	DF	MS	F	p-value
Month	159.868	1	159.868	17.750	0.1484
Residuals	9.007	1	9.007
SS = Sum of Squares; DF = Degrees of freedom; MS = Mean Square; Calculated using type III SS.

Visualisation of regression model

The blue line shows the best line of fit with shading representing 95% confidence intervals, while holding all other covariates constant. The dots show partial residuals, which reflect the observed data adjusted for all other predictors except the one being plotted.

Differences

Forest plot showing original and reproduced coefficients and 95% confidence intervals for Captive male Weight (KG) 25-30 months

Change in regression coefficients

term	O_B	R_B	Change.B	reproduce.B
Intercept		47.4605
Month	2.66	3.5526	0.8926	Not Reproduced
O_B = original B; R_B = reproduced B; Change.B = change in R_B - O_B; Reproduce.B = B reproduced.

Conclusion computational reproducibility

This model was not computationally reproducible.

Sensitivity analysis

As this model was not computationally reproducible, inferential reproducibility was not considered, since the original analyses could not be reproduced and therefore, statistical assumptions could not be meaningfully compared or interpreted.