Based on the operator product expansion, the perturbative and nonperturbative contributions to the polarized Bjorken sum rule (BSR) can be separated conveniently, and the nonperturbative one can be fitted via a proper comparison with the experimental data. In the paper, we first give a detailed study on the pQCD corrections to the leading-twist part of BSR. Basing on the accurate pQCD prediction of BSR, we then give a novel fit of the non-perturbative high-twist contributions by comparing with JLab data. Previous pQCD corrections to the leading-twist part derived under conventional scale-setting approach still show strong renormalization scale dependence. The principle of maximum conformality (PMC) provides a systematic and strict way to eliminate conventional renormalization scale-setting ambiguity by determining the accurate alpha _s-running behavior of the process with the help of renormalization group equation. Our calculation confirms the PMC prediction satisfies the standard renormalization group invariance, e.g. its fixed-order prediction does scheme-and-scale independent. In low Q^2-region, the effective momentum of the process is small and in order to derive a reliable prediction, we adopt four low-energy alpha _s models to do the analysis, i.e. the model based on the analytic perturbative theory (APT), the Webber model (WEB), the massive pQCD model (MPT) and the model under continuum QCD theory (CON). Our predictions show that even though the high-twist terms are generally power suppressed in high Q^2-region, they shall have sizable contributions in low and intermediate Q^2 domain. Based on the more accurate scheme-and-scale independent pQCD prediction, our newly fitted results for the high-twist corrections at Q^2=1;mathrm{GeV}^2 are, f_2^{p-n}|_{mathrm{APT}}=-0.120pm 0.013, f_2^{p-n}|_mathrm{WEB}=-0.081pm 0.013, f_2^{p-n}|_{mathrm{MPT}}=-0.128pm 0.013 and f_2^{p-n}|_{mathrm{CON}}=-0.139pm 0.013; mu _6|_mathrm{APT}=0.003pm 0.000, mu _6|_{mathrm{WEB}}=0.001pm 0.000, mu _6|_mathrm{MPT}=0.003pm 0.000 and mu _6|_{mathrm{CON}}=0.002pm 0.000, respectively, where the errors are squared averages of those from the statistical and systematic errors from the measured data.