Area under the curve (AUC)-directed vancomycin therapy is recommended, but Bayesian AUC estimation in critically ill children is difficult due to inadequate methods for estimating kidney function. We prospectively enrolled 50 critically ill children receiving IV vancomycin for suspected infection and divided them into model training (n = 30) and testing (n = 20) groups. We performed nonparametric population PK modeling in the training group using Pmetrics, evaluating novel urinary and plasma kidney biomarkers as covariates on vancomycin clearance. In this group, a two-compartment model best described the data. During covariate testing, cystatin C-based estimated glomerular filtration rate (eGFR) and urinary neutrophil gelatinase-associated lipocalin (NGAL; full model) improved model likelihood when included as covariates on clearance. We then used multiple-model optimization to define the optimal sampling times to estimate AUC24 for each subject in the model testing group and compared the Bayesian posterior AUC24 to AUC24 calculated using noncompartmental analysis from all measured concentrations for each subject. Our full model provided accurate and precise estimates of vancomycin AUC (bias 2.3%, imprecision 6.2%). However, AUC prediction was similar when using reduced models with only cystatin C-based eGFR (bias 1.8%, imprecision 7.0%) or creatinine-based eGFR (bias −2.4%, imprecision 6.2%) as covariates on clearance. All three model(s) facilitated accurate and precise estimation of vancomycin AUC in critically ill children.