OBJECTIVES: Digoxin is a well-known probe for the activity of P-glycoprotein. The objective of this work was to apply different methods for covariate selection in non-linear mixed-effect models to study the relationship between the pharmacokinetic parameters of digoxin and the genotype for two major exons located on the multi-drug-resistance 1 (MDR1) gene coding for P-glycoprotein. METHODS: Thirty-two healthy volunteers were recruited in three pharmacokinetic drug interaction studies. The data after a single oral administration of digoxin alone were pooled. All subjects were genotyped for the MDR1 C3435T and G2677T/A genotypes. The concentration-time profile of digoxin was established using 12-16 blood samples taken between 15 min and 72 h after administration. We modelled the pharmacokinetics of digoxin using non-linear mixed-effect models. Parameter estimation was performed using the stochastic approximation EM method (SAEM). We used three methods to select the covariate model: selection from a full model using Wald tests, forward inclusion using the log-likelihood ratio test and model selection using the Bayesian Information Criterion. RESULTS: The three covariate inclusion methods led to the same final model. Carriers of two T alleles for the C3435T polymorphism in exon 26 of MDR1 had a lower apparent volume of distribution than carriers of a C allele. The only other covariate effect was a shorter absorption time-lag in women. CONCLUSION: The apparent volume of distribution of digoxin is lower in TT subjects, probably reflecting differences in bioavailability. Non-linear mixed-effect models can be useful for detecting the influence of covariates on pharmacokinetic parameters.