We present a comprehensive new global QCD analysis of unpolarized parton distribution functions (PDFs) based upon proton, deuteron and $A=3$ data, including the latest inclusive deep-inelastic scattering (DIS) measurements from Jefferson Lab at high Bjorken-$x$. Using the JAM Bayesian Monte Carlo framework, we systematically explore the stability of the PDFs with respect to variations in the cuts on the invariant mass $W$ of the DIS final state, the implementation of target mass and higher twist corrections, as well as on the nuclear wave functions for the $A=2$ and 3 data. We find the $u$ and $d$ quark PDFs (and the $d/u$ ratio) are relatively stable up to $x\approx 0.8$, and able to describe DIS data down to $W^2=3.5$ GeV${}^2$ and $Q^2=m_c^2$. Within the collinear factorization framework, the fitted higher twist corrections to DIS are found to be positive, and largely isospin independent. The description of the nuclear data also requires nonzero isoscalar and isovector nucleon off-shell PDF contributions, which gives specific predictions for the ratio, $R_D$, of deuteron to isoscalar nucleon structure functions.