Understanding time-dependent effects of charged vacancies on electromechanical responses of materials is at the forefront of research for designing materials exhibiting metal-insulator transition, and memresistive behavior. A Rayleighian approach is used to develop a model for studying the non-linear kinetics of reaction leading to generation of vacancies and electrons by the dissociation of vacancy-electron pairs. Also, diffusion and elastic effects of charged vacancies are considered to model polarization-electric potential and strain-electric potential hysteresis loops. The model captures multi-physics phenomena by introducing couplings among polarization, electric potential, stress, strain, and concentrations of charged (multivalent) vacancies and electrons (treated as classical negatively charged particles), where the concentrations can vary due to association-dissociation reactions. Derivation of coupled time-dependent equations based on the Rayleighian approach is presented. Three limiting cases of the governing equations are considered highlighting effects of 1) non-linear reaction kinetics on the generation of charged vacancies and electrons, 2) the Vegard's law (i.e., the concentration-dependent local strain) on asymmetric strain-electric potential relations, and 3) coupling between a fast component and the slow component of the net polarization on the polarization-electric field relations. The Rayleighian approach discussed in this work should pave the way for developing a multi-scale modeling framework in a thermodynamically consistent manner while capturing multi-physics phenomena in ferroelectric materials.