We construct z \sim 6 - 7 , 8 , and 9 faint Lyman break galaxy samples ( 334 , 61 , and 37 galaxies , respectively ) with accurate size measurements with the software glafic from the complete Hubble Frontier Fields cluster and parallel fields data . These are the largest samples hitherto and reach down to the faint ends of recently obtained deep luminosity functions . At faint magnitudes , however , these samples are highly incomplete for galaxies with large sizes , implying that derivation of the luminosity function sensitively depends on the intrinsic size–luminosity relation . We thus conduct simultaneous maximum-likelihood estimation of luminosity function and size–luminosity relation parameters from the observed distribution of galaxies on the size–luminosity plane with the help of a completeness map as a function of size and luminosity . At z \sim 6 - 7 , we find that the intrinsic size–luminosity relation expressed as r _ { \textrm { e } } \propto L ^ { \beta } has a notably steeper slope of \beta = 0.46 ^ { +0.08 } _ { -0.09 } than those at lower redshifts , which in turn implies that the luminosity function has a relatively shallow faint-end slope of \alpha = -1.86 ^ { +0.17 } _ { -0.18 } . This steep \beta can be reproduced by a simple analytical model in which smaller galaxies have lower specific angular momenta . The \beta and \alpha values for the z \sim 8 and 9 samples are consistent with those for z \sim 6 - 7 but with larger errors . For all three samples , there is a large , positive covariance between \beta and \alpha , implying that the simultaneous determination of these two parameters is important . We also provide new strong lens mass models of Abell S1063 and Abell 370 , as well as updated mass models of Abell 2744 and MACS J0416.1 - 2403 .