Individual variations in drug efficacy and side effects pose serious problems in medicine. These variations are influenced by factors such as drug-metabolizing enzymes, drug transporters, and drug targets (e.g., receptors). For many medications, these factors can be attributed to genetic polymorphisms [1, 2]. Indeed, these genomic biomarkers are sometimes used to improve drug responses and reduce side effects by controlling the medication or dose according to the patient’s genotype [3, 4].
However, only a few of these genomic biomarkers have been validated. For this reason, many pharmacogenomics (PGx) studies have been launched around the world. The purpose of these studies is to identify genes that affect drug-metabolizing enzymes, drug transporters, and drug targets. Therefore, pharmacokinetics (PK) studies that include analyses of single-nucleotide polymorphisms (SNPs) as genomic markers in candidate-gene or genome-wide studies can be used to identify these genes. The availability of powerful array-based SNP-typing platforms has facilitated genome-wide studies, which have become a standard strategy. Such platforms make available to researchers genotype data for 100,000–4,300,000 SNPs.
In PK studies, it is common to apply compartmental models, which are often nonlinear models that include several PK parameters, in order to describe the profiles of drug concentrations in blood . Because drug concentrations in blood are usually related to drug efficacy and side effects, via their interactions with drug-metabolizing enzymes, drug transporters, and drug targets, differences in PK parameters indicate differences in effectiveness and toxicity. Therefore, one object of PGx studies is to identify genes associated with PK parameters.
Because drug-concentration data is measured from multiple subjects in PGx studies, inter-individual variability in PK parameters should be considered. Such data is referred to as population PK data. If the impact of inter-individual variability in model parameters is ignored, no statistically valid inference is possible. The mixed-effects model, which includes both fixed and random effects, is one method that accounts for inter-individual variability. Inter-individual variability in model parameters is modeled as random effects with strong parametric assumptions about the random-effects distribution. These models have often been used in analysis of longitudinal data , and they represent a useful method for accounting for inter-individual variability. Moreover, the nonlinear mixed-effects model (NLMM), an extension of the mixed-effects model to nonlinear functions, is often used to analyze population PK data [7, 8, 9, 10, 11].
The association between PK parameters and SNPs are typically analyzed using an NLMM [12, 13, 14, 15] in conjunction with population PK data. However, applying an NLMM to large-scale data can be problematic for the following reasons:
Computation time: NLMMs can be computationally intensive, because these models must compute the marginal log-likelihood by integrating out random effects . In NLMMs, inferences about model parameters are based on the marginal log-likelihood function, which includes a multiple integral with respect to the unobservable random effects. Because the regression functions are non-linear, the integral in the marginal log-likelihood function has no closed form, and it is necessary to compute the integral. To address this issue, various methods have been proposed to compute the integral approximation. However, these methods are computationally inefficient.
Convergence of iterative calculations: For instance, a major statistical software package, the SAS NLMIXED Procedure (SAS Institute, Inc., Cary, North Carolina) with adaptive Gauss–Hermite quadrature, is now used to approximate the maximum marginal log-likelihood [16, 17]. These computations are based on iterative calculations; for complex models, however, these calculations may not converge [18, 19]. If iterative calculations do not converge, we derive no information from valuable data.
Random-effects misspecification: Random-effects misspecification leads to bias in parameter estimates of the regression coefficients, and slightly inflates type-I error rates of tests for the regression coefficients in generalized mixed-effects models [20, 21, 22, 23] and NLMMs . Therefore, careful model building and checking are needed for each of the 100,000–4,300,000 analyses, but in practice this may be difficult to apply.
In conclusion, it seems that these three problems occur in association with a strong assumption of random effects.
Therefore, we consider a new method that uses a potentially misspecified model to avoid the strong assumptions of the random-effects distribution. Misspecified models are useful and powerful tools for studying the behavior of estimators under model misspecification. Model misspecification means that an incorrect working model is used for estimation. In this paper, we consider a “true” model that includes inter-individual variability in model parameters as fixed-effect parameter vectorsfor each subject (), and a “working” model that misspecifies the presence of inter-individual variability in model parameters as a common parameter vector . In this paper, we describe a new interpretation of the estimator as a weighted average of the individual parameter vectors . The proposed method allows for computation that is faster than NLMM by a factor of 100, performs stable computations, and is robust for various structures of individual variations, because it is based on misspecified fixed-effect models instead of random-effect models.
A general theory of misspecified models has been proposed by White  for maximum-likelihood methods, and by Yi and Reid  for estimating equations. In both papers, the authors demonstrate the asymptotic normality of estimators of working-model parameters under mild conditions. The method we propose uses generalized estimating equations (GEE) of a working model based on Yi and Reid’s result. GEE has been widely used in regression analyses of the generalized linear models with correlated response, such as repeated-measurement data [26, 27, 28]. Under mild conditions, the estimator from GEE of a misspecified model is consistent and asymptotically normal. Therefore, the proposed method may be applied to correlated-response data with inter-individual variability in model parameters, which includes a wide range of applications. In this paper, however, we focus on the problem of estimating PK parameters in the presence of inter-individual variability in PGx studies, i.e., the motivating example.
The proposed method focuses only on estimating fixed effect parameter vectors, because one object of PGx studies is to identify genes associated with PK parameters. Other parameters are nuisance parameters. Therefore, a misspecified model that gives an estimator for a weighted average of fixed effect parameter vectors is “intentionally” used. As a result, the proposed estimator is different from the estimator of NLMMs. Marginal (or population-averaged) models and mixed (or subject-specific) models can lead to a different estimator in non-linear settings . The proposed method relates to marginal models.
For each SNP, there are three genotypes for each locus: the “aa”, “Aa”, and “AA” genotypes, where “a” is the major allele and “A” is the minor allele. There are often a considerable number of genes for which the frequency of the minor homozygous genotype, “AA”, is very small, because a SNP is defined as a mutation involving a single DNA base substitution that is observed with a frequency of at least 1% in a population. Therefore, for valid statistical inference, a small–sample correction is needed. To address the small–sample size problem, we propose a Wald-type test and an asymptotic -test for determining the effects of a genetic polymorphism on PK parameters (see Section 3.5).
In Section 2, motivating data and issues of NLMM are introduced. In Section 3, misspecified models and the proposed method are presented, and some of the proposed method’s theoretical properties are discussed. In Section 4, we study the performance of the proposed method using simulations. In Section 5, we present the application of the method to published experimental data. We present our concluding remarks in Section 6.
2 Motivating example
The motivation for this paper stems from an PGx study data  on gemcitabine (2’,2’-difluorodeoxycytidine), which is a nucleoside anticancer drug. The study was designed to screen for genes related to the PK of gemcitabine. The participants consisted of 233 gemcitabine-naive cancer patients (mainly with pancreatic carcinoma). For the PK analysis, heparinized blood samples were taken before administration and at 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, and 4.5 hours after the beginning of the administration. The dose was adjusted according to the surface area of the body of each subject. A total of 109,365 gene-centric SNPs were genotyped using the Sentrix Human-1 Genotyping BeadChip (Illumina Inc., San Diego, CA).
Because the main object of this PGx study was to screen PK-related genes, the SNP genotype effect on PK parameters was modeled using a compartmental model. Compartmental models, which are derived from differential equations that describe drug kinetics, are nonlinear models with several PK parameters. It is common to apply such models in order to describe the profiles of drug concentrations in blood.
In general, analyses of genome-wide data use appropriate statistical methods to investigate the association between an outcome variable and a set of SNPs 
. Based on the results, favorable SNPs that strongly associate with the outcome variable are screened with appropriate criteria (e.g., the Bonferroni adjustment, false-discovery rate, etc.). The appropriate statistical methods are determined by the nature of the outcome variable and the study design (e.g., trend tests for odds ratios are used in case-control studies), and these analyses are commonly performed one by one for each SNP. In genome-wide PGx studies, to identify genes that associate with PK parameters, analyses of the associations between PK parameters and SNPs are applied to population PK data. We consider that the SNP genotype effect reflects the difference in average PK parameters between different genotypes. Moreover, because the PK data include multiple individuals, the data are population PK data, and we must consider the impact of inter-individual variability in PK parameters.
Now, we introduce notations and describe the data structure. Suppose we have a plasma drug concentration dataset with subjects and a genotype dataset with SNPs. For each subject (), there is a random ()-dimensional vector , a covariate ()-dimensional matrix , and an ()-dimensional vector representing time after starting measurement. The covariate matrix includes genotype data of the -th SNP (). Because similar analysis is repeated for each of the SNPs, we shall write instead of for simplicity. Note that the superscript “” indicates the transpose of a matrix or a vector.
In the literature to date, several studies have analyzed a genetic polymorphism in relation to population PK data using NLMMs [12, 13, 14, 15]. In NLMMs, it is often assumed that arises from the non-linear model,
where is a compartmental model function that is non-linear in its PK parameters , is a error vector, and are vectors of fixed effects and random effects, is a design matrix for random effects, and is a covariance matrix. The
are assumed to have a multivariate normal distribution with mean vectorand a covariance matrix . NLMM incorporates unmeasured random effects into the compartmental model function to account for inter-individual variability in the PK parameters .
Because two-compartment models have been widely used in gemcitabine PK analyses [32, 33], we fitted a two-compartment constant intravenous-infusion model (Figure 1). For the gemcitabine PGx data, the function form of is as follows:
because the PK parameters are restricted to be positive, and in practice are empirically log-transformed 
. Commonly, the three genotypes are considered in evaluating the relationship between SNPs and the PK parameters. We use the dummy variablesand for the covariate matrix . Let , , or denote that the -th subject has the genotype “aa”, “Aa”, or “AA”, respectively. We assume that the effect of a SNP on the PK parameter can be described by the following relationship;
where is an intercept parameter for each PK parameter, and , , , and are the effect parameters of a SNP genotype “Aa” and “AA” for each PK parameter, respectively.
In order to evaluate the effect parameters of a SNP genotype, a test for the effect for the single genotype (e.g., vs.
) and a multiple degrees-of-freedom test (e.g.,vs. ) can be considered. When the null hypotheses, , of these tests are rejected at an appropriate significance level, it can be concluded that the SNP affects the profiles of drug concentrations in blood.
However, applying an NLMM to large-scale data, such as genome-wide PGx studies, can be problematic. (i) NLMMs can be computationally intensive, because these models must compute the marginal log-likelihood by integrating out random effects . (ii) These computations are based on iterative calculations, and may not converge in complex models [18, 19]. (iii) Random-effects misspecification leads to bias in parameter estimates of the regression coefficients and slightly inflates type-I error rates of tests for the regression coefficients in an NLMM .
Therefore, we consider an alternative approach that relates to a marginal modeling approach that avoids the specification of random effects. The approach potentially results in misspecification of the model for the parameters of interest. Therefore, we evaluated the estimator based on the proposed approach via a misspecified model.
3 Estimation and inference
3.1 Misspecified models
Misspecified models are useful and powerful tools for studying a behavior of estimators under model misspecification. Model misspecification means that an incorrect working model is used for estimation. The proposed method “intentionally” uses an incorrect working model that gives an estimator for a weighted average of fixed effect parameter vectors.
The general theory of misspecified models has been proposed by White  for maximum likelihood methods, and by Yi and Reid  for estimating equations. White showed that the maximum likelihood estimator of a misspecified model converges to a constant vector
which minimizes the Kullback–Leibler divergence. A similar property holds in the case of estimating equations. Under mild conditions, Yi and Reid showed thatis asymptotically normally distributed with a mean vector and a covariance matrix that can be consistently estimated by the so-called sandwich estimator. However, because generally do not have a simple analytical form, we need to evaluate the properties of .
The literature includes several attempts to uncover the relation between the parameters of a true model and the estimators of model parameters from an incorrect working model. For example, misspecified models under non-proportional hazards models with a time-varying effect parameter have been discussed in semi-parametric survival models [35, 36]. Xu and O’Quigley  evaluate an asymptotic property of the estimator from a misspecified proportional-hazards model that replaces with a constant . They showed that the estimator
converges in probability to a constantthat is approximated by a weighted average of over time, , where
is the conditional variance of a stochastic process, and is a possibly time-dependent covariate. Xu and O’Quigley showed that the estimator can be interpreted as a weighted average of true parameters even when an incorrect working model is used.
In this paper, we consider GEE of a misspecified model. We assume a true model with inter-individual variability in model parameters as fixed-effect parameter vectors , and a working model that misspecifies the presence of inter-individual variability in model parameters as a common parameter vector . We demonstrate a new interpretation of the estimator as a weighted average of the individual parameter vectors .
3.2 Assumptions about the true model
To describe the true structure of the observations, we assume that the true cumulative distribution of is with density , where is a ()-dimensional vector of effect parameters with inter-individual variability as fixed effects, is a scale parameter, and is a variance model parameter vector.
The expectation of the observation is modeled as , where is a PK function that is nonlinear in its PK parameters , is an individual PK parameter ()-dimensional vector, and is modeled in linear form as .
Furthermore, we assume that the variance is modeled as , where is a known variance function that has the variance model parameter vector .
3.3 An estimator of a weighted average effect by GEE
Under the true distribution in Section 3.2, an average effect may be obtained heuristically by replacingwith a constant and then fitting to data. We consider GEE of a working model that has parameters . GEE is well known to be inadequate when the mean structure is misspecified. However, we will show that the estimator can be interpreted as a weighted average effect under the true model, in Section 3.4.
We define GEE of a potentially misspecified model as
is a working correlation ()-dimensional matrix that can depend on a parameter vector , and is a parameter ()-dimensional vector that is common between individuals. The true distribution has the effect parameter vector , which has inter-individual variability; however, this working model assumes no inter-individual variability. In equation (2), the individual PK parameter ()-dimensional vector is modeled as . Note that indicates a diagonal matrix with diagonal elements in parentheses.
Here, we denote the solution to equation (2) as . Yi and Reid  showed following theorems.
Theorem 1. Under the true model (the distribution function of is ), the estimator converges in probability to a constant vector as , where is a constant vector that satisfies the equation
Theorem 2. Under the true model, is asymptotically normal with a mean vector and a covariance matrix as , where
is a sandwich variance,
The solution to equation (3) can be interpreted as a weighted average of the individual parameter vectors . In equation (3) from Theorem 1, minimizes the distance between the true model and the misspecified model. For example, when is the score function, White  showed that minimizes the Kullback–Leibler divergence between the true model and the misspecified model.
However, do not have a simple analytical form. To evaluate the properties of , we consider a first-order Taylor expansion of the expectation of equation (2) around . As a result, we get
which is a weighted average of the individual parameter vectors with weights , where is an inverse matrix of the model-based variance of the -th subject. The derivation of equation (4) is shown in Appendix A.
Therefore, according to Theorem 1 and equation (4), the estimator from the working model can asymptotically estimate a weighted average of the individual parameter vectors . For instance, if , then ; moreover, if and , then . The model corresponding to (2) is a misspecified model. Nevertheless, the estimator converges to and can be interpreted as population-weighted average parameters. Furthermore, this approach is robust for various structures of inter-individual variability (see Section 4.1), because it does not require a strong assumption of a random-effect distribution.
3.5 A Wald-type test and an asymptotic -test
As indicated in Section 1, a small–sample correction is needed for valid statistical inference in genome-wide PGx studies that analyze SNP genotyping data. Therefore, we consider a Wald-type test based on distributions and an asymptotic -test instead of the asymptotic Wald chi-square tests.
We considered linear hypotheses of the form versus , where is a contrast-coefficient ()-dimensional vector.
We proposed a Wald-type test statistic;
is the estimator of the covariance matrix ,
is asymptotically normally distributed with mean 0 and variance 1 under the null hypothesis, and assuming that
follows a chi-square distribution withdegrees of freedom (d.f.), the test statistic is asymptotically -distributed with d.f., which must be estimated.
We applied the moment estimator of the d.f., as proposed by Fay and Graubard , where is a block-diagonal matrix, , , and is a block-diagonal matrix. We discuss a derivation of in Appendix B.
We also considered approximating the distribution of multiple degrees-of-freedom tests of versus , where is a contrast-coefficient ()-dimensional matrix and is the number of contrast-coefficient vectors that one wishes to test.
We proposed an asymptotic -test statistic:
The statistic is asymptotically -distributed with a numerator d.f. and a denominator d.f. that must be estimated. We applied the moment estimator of the denominator d.f.,
3.6 Bias correction for
The sandwich-variance estimator is biased downward under small–sample size conditions, as shown by Mancl and DeRouen  (see also [41, 42, 43]). As indicated in Section 1, a bias correction is also needed in genome-wide PGx studies.
To calculate , a product of the residual vector is used to estimate . However, using a first-order Taylor expansion of and around , and , we have , where , , and can be obtained by replacing by in the expression , , and ; ; and is an (
)-dimensional identity matrix. Replacingin equation (5) by gives the bias-corrected sandwich-variance estimator,
3.7 Example application: Compartment models and the effects of SNPs
In this study, we applied our methods to a genome-wide PGx study , introduced above in Section 2.
We introduce the true distribution of and a generalized estimating equation for the data. It is common to use a constant coefficient of variation (CV) model in PK data analysis [9, 45, 46, 47]. Under the constant CV model, the expectation and variance of the log-transformed random vector, , is modeled as , , and , where
is the log-transformed PK function, and is given by equation (1). Furthermore, the working correlation matrix is modeled as, .
Along with NLMM as shown in Section 2, in order to evaluate the association between PK parameters and SNPs, the individual PK parameters in GEE of a working model is modeled as
where the covariate matrix is =12
and the parameter vector is
Simulations were conducted to study the performance of the proposed Wald-type test and asymptotic -test for population PK data. The simulation conditions for population PK data were determined by reference to an actual genome-wide PGx study .
For simplicity, observed responses , which should have individual variations, were generated from the NLMM as follows:
where is a random-effect vector of the -th subject for each PK parameter, for which conditions are shown in Section 4.1 and Section 4.2. Further, we assumed , where is a log-transformed two-compartment constant intravenous-infusion PK function in equation (7) setting to 1400 mg and to 0.5 hours. The intercept terms of the log-transformed PK parameters were set to , , , and
; the standard deviationwas set to 0.27; and values of the remaining parameters are shown in Section 4.1 and Section 4.2. Note that these parameters were set based on a preliminary NLMM analysis of gemcitabine data without covariates; we assumed that the random-effect vector is normally distributed with a mean vector and a diagonal covariance matrix , because the random effects of the elimination parameter were too small. In the simulations, we changed the blood sampling points to 0.1, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, and 4.5 hours after drug administration. The covariate matrices in equation (8) were generated non-randomly by taking
In actual studies, the sample size of each genotype group is not controlled, but depends on allele frequency. Generally, these studies are likely to have unequal sample sizes for different genotypes and a minor-allele frequency (MAF) of less than 0.5, most commonly around 0.2 . When MAF is small, there are too few subjects homozygous for the minor allele. Therefore, the MAF was set to 0.25 or 0.50. If we let , , and denote the sample size of each genotype, then . We assumed that the population was in Hardy–Weinberg equilibrium, with the total sample size ; the sample size for each group was set to , , and for , and , , and for . The total sample size is not realistic for genome-wide PGx studies, but is sufficient to evaluate statistical performance.
For each data configuration, 1000 simulations were generated.
For each simulation, the generalized estimating equation of a misspecified model in equation (2) was fitted assuming the two-compartment constant intravenous-infusion model as shown in Section 3.7, and the NLMM was fitted assuming a normal random-effects model with adaptive Gauss–Hermite quadrature. We used a diagonal covariance matrix for the random-effect vector , which is normally distributed with a mean vector , and a covariance matrix for the NLMM, which is a commonly-used method (see Section 2).
In order to assess the statistical performance of tests for the effect of a SNP on PK parameters (e.g., vs. and vs.