The analytic inverse kinematics solutions you found do depend on those $0$ terms in your transformation matrices. Those values are, as you've implied, based on the $0$ and $90$ degree values for the axes relationships. The kinematic mapping process works for other $a_i$ and $d_i$ values, too, but it becomes difficult to find closed-form inverse solutions when those nice zeroes become other values.
The solution most often implemented is to use numerical solutions for the inverse kinematics after calibration. The closed-form, analytic, inverse kinematics are a good starting estimate for the numerical solver, but you must then account for the as-built errors with some (often iterative) approach to find the actual joint angles given the EE pose. There are many published approaches to the numerical solution methods.
EDIT (Industrial robots): For general-purpose industrial robots the above is rarely done. Instead of calibrating to eliminate as-built errors in the kinematics, the robots are usually calibrated to align the task space coordinate system with the robot's base coordinate system (resulting in a single additional transformation matrix of constant values). The end effector destination locations are then taught (in either task, end-effector, or joint coordinates). Because they use a "teach and repeat" method of programming, the small errors in the calibration are absorbed by the teaching process.
Additional calibrations of the robot may be accomplished by the robot manufacturer. In my experience those calibrations generally set precise values for the link lengths and encoder home, or zero, values - but not joint axis offset angles. It's possible that some robot manufacturers calibrate more precisely, then use numerical methods to find inverse solutions, but I have not seen this absolute precision needed in an industrial setting.
