Statistical analysis
To determine if the range shift (centroid distance: distance between the
current and future centroids) is associated with the amount of current
suitable habitat and elevational distribution, I ran two separate GLMMs
in R using the ‘glmmTMB’ package (Magnusson et al., 2017). To calculate
whether the range shift (centroid distance) is associated with the
current suitability distribution along the elevational gradient, I added
the centroid distance as the response variable and the interaction
between mean elevation (current) and migratory status as the explanatory
variable. I used the interaction of migratory status with the response
variable to account for how migratory species can adjust to the new
habitat faster than the non-migratory species and assess if Bangladeshi
migratory butterflies follow a similar trend. To calculate whether the
range shift is associated with the amount of suitable habitat (current),
I ran a similar model considering the range shift as the response
variable and the interactive effect of the amount of suitable habitat
(current) and species migratory status as the explanatory variable.
In both cases, I used the ‘nbinom2’ family to account for overdispersion
and butterfly families as a random factor. Finally, if any interaction
was significant, I conducted post hoc pairwise comparison tests using
the ‘emmeans’ R package (Lenth et al., 2019).
To assess the vulnerability status of the threatened and non-threatened
species, I ran two separate GLMMs (for centroid distance and elevation
shift) in R using the ‘glmmTMB’ package (Magnusson et al., 2017). To
test whether threatened species experienced more range shift or
elevation shift, I added centroid distance and elevation shift
periodically as the response variables and the threat status as the
explanatory variable. I used the ‘nbinom2’ family to account for
overdispersion and butterfly families as a random factor. Given that
there were negative values in the elevation shift (species for which the
mean elevation declined), before fitting the model, I transformed the
values by adding (1+ maximum negative values) for each species.