2.3 Data analysis
For this study, we performed all analyses in R 4.2.0 (R Development Core Team, 2022). We improved the residual normality and reduced the heteroscedasticity of our data by log-transforming (natural logarithm) female snout-vent length. To answer our central questions regarding the determinants of intraspecific variation in female body size across geographical gradients, we constructed linear mixed models using the lmer function implemented in the lme4 package (Bates et al., 2015). For each model, we used the lizard population origin as a random intercept; because this allows us to account for the non-independence of lizards within populations (Bolker et al., 2009).
First, we asked whether the geographical patterns of body size at the intraspecific level of lizards vary with altitude; to achieve this, we constructed a simple linear model of ln-transformed body sizes of lizards as the response variable with altitude (binned into three categories) as the predictor variable. We then determined the significance using the Anova function from the car package (Fox and Weisberg, 2019). Further, we constructed a post hoc test for our model using the emmeans function from theemmeans package (Lenth, 2019) to specifically test for differences in body size between our three levels of altitude.
we then fit univariate linear models of each environmental factor (annual mean temperature, annual mean precipitation, temperature seasonality, and precipitation seasonality) as the response variable with the altitude as the predictor variable in each case to understand the relationship between climatic conditions and altitude in our study.
Next, we asked whether female body size varied with resource availability and/or in response to the changing climatic and seasonal conditions across altitudes. However, we could not achieve a single model for all our predictor variables due to the significant collinearity between net primary productivity and annual mean temperature (r  = 0.70; p <0.0001) and net primary productivity and annual mean precipitation (r  = 0.96;p <0.0001), which resulted in all linear mixed models failing to converge. Thus, we fit two alternative models to explore climate-body size relationships: (1) with net primary productivity, temperature seasonality, and precipitation seasonality; and (2) with annual mean temperature, annual mean precipitation, temperature seasonality, precipitation seasonality and altitude as our predictor variables. In all our fitted models, we included altitudes as a covariate because climatic and seasonal changes that influence the life-history traits of species significantly vary across geographic gradients such as altitudes (Hille and Cooper, 2015; Laiolo and Obeso, 2015).
We binned elevation in our models because of the discontinuous elevational distribution of the lizards (Fig. 1). Due to the discontinuous nature of the elevational data, modeling it as continuous led to model failure. Therefore, we follow the approach of Bhat et al. (2020) to categorize the lizard populations in our data into three altitudinal levels (low altitudes: <1,000m; mid altitudes: 1,000–2,000 m; high altitudes: >2,000 m asl).