Package {twoPhaseGAS}

This code a function that generates a sample Y, Z, G given some initial parameters.

Description

This code a function that generates a sample Y, Z, G given some initial parameters.

Usage

DataGeneration_TPD(
  Beta0 = 2,
  Beta1 = 0.5,
  Disp = 1,
  N = 5000,
  LD.r = 0.75,
  P_g = 0.2,
  P_z = 0.3,
  family = gaussian(),
  tao = 2/5
)

Arguments

Value

A dataframe with complete data Y, G, Z, S, where G and Z come from the same haplotype determined by P_g, P_z and LD.r; Y is generated from E[Y|G] = family$linkinv(Beta0 + Beta1 x G). S is a 3 level variable determined by Y, Beta1, Disp and tao (only returned when family is gaussian). the function iterates across generated datasets until Z and Y are associated at a suggestive genome wide threshold of p<=1e-05.

Examples

data = DataGeneration_TPD()

Likelihood ratio test for twoPhaseSPML models

Description

Compares two nested twoPhaseSPML models via a likelihood ratio test (LRT), mirroring the interface of lmtest::lrtest().

Usage

lrtest(object, object2, ...)

## S3 method for class 'twoPhaseSPML'
lrtest(object, object2, ...)

Arguments

Details

The test statistic is

LRT = -2 \bigl(\ell_0 - \ell_1\bigr),

where \ell_0 and \ell_1 are the maximised log-likelihoods of the smaller and larger models, respectively. Under the null hypothesis that the additional parameters are zero, LRT follows a \chi^2 distribution with degrees of freedom equal to the difference in the number of estimated regression coefficients between the two models.

The function determines model order automatically: whichever of object and object2 has fewer regression coefficients is treated as the null model, regardless of the order in which they are supplied.

Value

An object of class "lrtest.twoPhaseSPML", which is a data.frame with one row per model and the following columns:

#Df

number of estimated regression coefficients.

LogLik

maximised log-likelihood.

Df

change in degrees of freedom relative to the previous row (NA for the first row).

Chisq

LRT statistic (NA for the first row).

Pr(>Chisq)

two-sided p-value from the \chi^2 distribution (NA for the first row).

The object also carries a "heading" attribute (a character vector) and is printed via print.lrtest.twoPhaseSPML.

See Also

twoPhaseSPML, print.lrtest.twoPhaseSPML

Examples

data.table::setDTthreads(1)
data <- DataGeneration_TPD(N=200)
set.seed(1)
R <- rep(0L, nrow(data))
R[sample(nrow(data), 175)] <- 1L
data[R == 0L, c("G1")] <- NA

# Null model (no G1)
res_Ho <- twoPhaseSPML(
  formula = Y ~ Z,
  miscov  = ~ G1,
  auxvar  = ~ Z,
  data   = data
)

# Alternative model (with G1)
res_Ha <- twoPhaseSPML(
  formula = Y ~ Z + G1,
  miscov  = ~ G1,
  auxvar  = ~ Z,
  data   = data
)

lrtest(res_Ho, res_Ha)

Genetic-algorithm optimal phase 2 design for a two-phase genetic association study.

Description

Uses a genetic algorithm (GA) to search over individual-level phase 2 selection indicators R_i \in \{0,1\} for the subset of size n that maximises the anticipated Fisher information for the regression parameters of interest. Wraps kofnGA. When a parallel cluster is supplied, multiple independent GA runs are executed in parallel and the best solution across all runs is returned.

Usage

optimTP.GA(
  formula,
  miscov,
  auxvar,
  family,
  n,
  data,
  beta,
  p_gz,
  disp = NULL,
  ga.popsize,
  ga.propelit,
  ga.proptourney,
  ga.ngen,
  ga.mutrate,
  ga.initpop = NULL,
  optimMeasure,
  K.idx = NULL,
  seed = 1,
  verbose = 0,
  cluster = NULL
)

Arguments

Details

The fitness function evaluates the anticipated Fisher information matrix (FIM) for a candidate selection indicator R:

\mathrm{FIM}(R) = I_3/N + \sum_{i: R_i=1}(I_{2,i} - I_{3,i})/N,

where I_3 is the empirical score outer-product estimator (the numerically stable approximation to the full Louis (1982) observed FIM valid for large N, since I_1 \approx I_2 in the genotype expansion context), and I_{2,i}, I_{3,i} are the per-individual complete-data score outer product and observed-data score outer product, respectively. Under H_0 (beta[G] = 0) the score-test Schur complement is optimised; under H_a the full inverse FIM is optimised.

Value

An object of class "kofnGA" (see kofnGA) with the following elements:

bestsol

Integer vector of length n giving the row indices (in data) of the selected phase 2 individuals.

bestobj

Numeric value of the optimality criterion at bestsol.

pop

Final population matrix (ga.popsize rows, n columns), sorted by fitness.

obj

Objective values for each row of pop.

old

List with per-generation best, obj, and avg records.

time.elapsed

Named numeric vector with the pre-processing (preproc) and GA (GA) times in seconds.

See Also

kofnGA, twoPhaseDesign, optimTP.LM

Examples

data <- DataGeneration_TPD(N = 500)
data <- data[, c("Y", "Z")]  # phase-1 data: no G column
p_gz <- data.frame(Z = c(0, 0, 0, 1, 1, 1, 2, 2, 2), # theoretical joint distribution
        G1 = c(0, 1, 2, 0, 1, 2, 0, 1, 2), # between G and Z (LD = 0.75, MAF_G = 0.2, MAF_Z = 0.3)
        q = c(0.486, 0.004, 0, 0.143, 0.276, 0.001, 0.011, 0.04, 0.039)) 

res_ga <- optimTP.GA(
  formula      = Y ~ Z + G1,
  miscov       = ~ G1,
  auxvar       = ~ Z,
  family       = gaussian(),
  n            = 100,
  data         = data,
  beta         = c(2, 0.5, 0.1),
  p_gz         = p_gz,
  ga.popsize   = 40,
  ga.propelit  = 0.1,
  ga.proptourney = 0.1,
  ga.ngen      = 50,
  ga.mutrate   = 0.01,
  optimMeasure = "Par-spec",
  K.idx        = 2L
)
cat("Best objective:", res_ga$bestobj, "\n")
cat("Selected indices (first 10):", head(res_ga$bestsol, 10), "\n")

## Parallel processing (2 workers)
cl <- parallel::makeCluster(2L)
res_par <- optimTP.GA(
  formula      = Y ~ Z + G1,
  miscov       = ~ G1,
  auxvar       = ~ Z,
  family       = gaussian(),
  n            = 100,
  data         = data,
  beta         = c(2, 0.5, 0.1),
  p_gz         = p_gz,
  ga.popsize   = 40,
  ga.propelit  = 0.1,
  ga.proptourney = 0.1,
  ga.ngen      = 50,
  ga.mutrate   = 0.01,
  optimMeasure = "Par-spec",
  K.idx        = 2L,
  cluster      = cl
)
parallel::stopCluster(cl)
cat("Best objective (parallel):", res_par$bestobj, "\n")

Lagrange-multiplier optimal phase 2 design for a two-phase genetic association study.

Description

Computes stratum-specific sampling fractions that maximise the anticipated Fisher information for the regression parameters of interest subject to a fixed total phase 2 sample size. The optimisation is performed by the COBYLA algorithm (nloptr) with a derivative-free fallback (dfoptim::hjkb) if COBYLA fails to converge.

Usage

optimTP.LM(
  formula,
  miscov,
  auxvar,
  strata,
  family,
  n,
  data,
  beta,
  p_gz,
  disp = NULL,
  optimMeasure,
  K.idx = NULL,
  min.nk = NULL,
  logical.sub = NULL
)

Arguments

Details

The function computes, for each stratum defined by strata, the stratum-specific sampling fraction \pi_k = n_k / N_k that maximises the chosen information criterion. The optimisation is constrained so that (i) the total selected sample equals n, (ii) no stratum is over-sampled (n_k \le N_k), and (iii) the min.nk lower bounds are respected.

The baseline Fisher information uses the empirical score outer-product estimator I_3/N (see optimTP.GA for details).

Value

A data.frame with one row per stratum (as defined by strata) and the following columns from xtabs(strata, data):

(strata variables)

The stratification variable values.

Freq

Stratum size N_k.

prR_cond_optim

Optimal conditional sampling probability n_k / n (proportional allocation weight).

convergence

Convergence code from the solver (0 = success; codes from nloptr::cobyla or dfoptim::hjkb).

See Also

twoPhaseDesign, optimTP.GA

Examples

data <- DataGeneration_TPD(N = 1000)
data <- data[, c("Y", "Z")]  # phase-1 data: no G column
p_gz <- data.frame(Z = c(0, 0, 0, 1, 1, 1, 2, 2, 2), # theoretical joint distribution
        G1 = c(0, 1, 2, 0, 1, 2, 0, 1, 2), # between G and Z (LD = 0.75, MAF_G = 0.2, MAF_Z = 0.3)
        q = c(0.486, 0.004, 0, 0.143, 0.276, 0.001, 0.011, 0.04, 0.039)) 
        
data$S <- as.numeric(cut(data$Y,
            breaks = quantile(data$Y, probs = c(0, 0.4, 0.8, 1)),
            include.lowest = TRUE))

res_lm <- optimTP.LM(
  formula     = Y ~ Z + G1,
  miscov      = ~ G1,
  auxvar      = ~ Z,
  strata      = ~ Z + S,
  family      = gaussian(),
  n           = 200,
  data        = data,
  beta        = c(2, 0.5, 0.1),
  p_gz        = p_gz,
  optimMeasure = "Par-spec",
  K.idx       = 2L
)
head(res_lm)

Print method for lrtest.twoPhaseSPML objects

Description

Print method for lrtest.twoPhaseSPML objects

Usage

## S3 method for class 'lrtest.twoPhaseSPML'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

x is returned invisibly.

See Also

lrtest.twoPhaseSPML

Print method for twoPhaseSPML objects

Description

Prints a concise, glm-style summary of a fitted twoPhaseSPML model.

Usage

## S3 method for class 'twoPhaseSPML'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Details

Under the null model (miscov absent from formula), score statistics (Sobs and Sexp) are displayed. Under the alternative model, Wald statistics (Wobs and Wexp) together with standard errors and p-values for the missing-by-design covariate(s) are displayed. The coefficient table always shows the estimated regression coefficients; standard errors and p-values are shown only when var_theta is available (alternative model).

Value

x is returned invisibly.

See Also

twoPhaseSPML

Examples

data.table::setDTthreads(1)
data <- DataGeneration_TPD(N=200)
set.seed(1)
R <- rep(0, nrow(data))
R[sample(nrow(data), 175)] <- 1
data[R == 0, c("G1")] <- NA

res_Ho <- twoPhaseSPML(
  formula = Y ~ Z,
  miscov  = ~ G1,
  auxvar  = ~ Z,
  data   = data
)
print(res_Ho)

Internal function to optimize over a range of maf, LD and betas.

Description

Internal function to optimize over a range of maf, LD and betas.

Usage

twoPhase(
  beta = c(1, 1, 1),
  maf_G = 0.1,
  LD = 0.3,
  data = NA,
  n2 = NA,
  design_formula = Y ~ G + Z,
  family = gaussian(),
  useGeneticAlgorithm = FALSE
)

Arguments

Value

sth

Examples

data = twoPhase()

Compute sample allocations for a two phase study design.

Description

Compute sample allocations for a two phase study design.

Usage

twoPhaseDesign(
  beta,
  maf_G,
  LD,
  data,
  n2,
  design_formula = Y ~ G + Z,
  family = gaussian(),
  S,
  perc_Y = c(1/5, 4/5),
  p_gz,
  design = c("RDS", "LM", "GA", "PPS", "BAL", "COM", "TZL"),
  ndraws = 10,
  optimCriterion = c("Par-spec", "A-opt", "D-opt"),
  overallMethod = c("med-max", "cumm")
)

Arguments

Value

sth

Select samples for phase 2 under heuristic designs.

Description

Select samples for phase 2 under heuristic designs.

Usage

twoPhaseHeuristic(
  design = c("pps", "bal", "comb"),
  ndraws = 1,
  data = NA,
  n2 = NA,
  family = gaussian(),
  S = NULL,
  perc_Y = c(1/5, 4/5)
)

Arguments

Value

sth

Performs inference on two-phase studies data via semiparametric maximum likelihood.

Description

Performs inference on two-phase studies data via semiparametric maximum likelihood.

Usage

twoPhaseSPML(
  formula,
  miscov,
  auxvar,
  family = gaussian,
  data,
  start.values = NULL,
  control = list()
)

Arguments

Details

The function implements the semiparametric maximum likelihood (SPML) estimator for two-phase genetic association studies described by Espin-Garcia et al. The EM algorithm alternates between an E-step that computes observation weights based on the current estimate of the joint distribution of the missing covariate and the auxiliary variable, and an M-step that updates the regression coefficients and the non-parametric distribution via iteratively re-weighted least squares and a closed-form update, respectively.

When miscov is absent from formula, inference is performed under H_0 and score statistics are returned. When miscov is included in formula, inference is performed under H_a and Wald statistics are returned.

Value

An object of class "twoPhaseSPML", which is a list containing:

theta

named numeric vector of estimated regression coefficients.

qGZ

data frame with the estimated joint distribution of miscov and auxvar.

dispersion

estimated dispersion parameter (variance for Gaussian; 1 for other families).

iter

integer, number of EM iterations performed.

ll

log-likelihood at convergence.

df

degrees of freedom for the test statistic.

Sobs

(null model only) observed score test statistic.

Sexp

(null model only) expected score test statistic.

var_theta

(alternative model only) named numeric vector of marginal variances for theta.

Wobs

(alternative model only) observed Wald test statistic.

Wexp

(alternative model only) expected Wald test statistic.

qG

data frame with the estimated marginal distribution of miscov.

formula

the formula used.

miscov

the miscov formula used.

family

the family object used.

null.model

logical, TRUE when the model was fitted under H_0.

Examples

data.table::setDTthreads(1)
data <- DataGeneration_TPD(N=200)
set.seed(1)
R <- rep(0, nrow(data))
R[sample(nrow(data), 175)] <- 1   # random phase 2 subsample of 500
data[R == 0, c("G1")] <- NA

# Null model (score test)
res_Ho <- twoPhaseSPML(
  formula = Y ~ Z,
  miscov  = ~ G1,
  auxvar  = ~ Z,
  data   = data
)
print(res_Ho)

# Alternative model (Wald test)
res_Ha <- twoPhaseSPML(
  formula = Y ~ Z + G1,
  miscov  = ~ G1,
  auxvar  = ~ Z,
  data   = data
)
print(res_Ha)

Auxiliary for Controlling the EM algorithm in twoPhaseSPML.

Description

Auxiliary for Controlling the EM algorithm in twoPhaseSPML.

Usage

twoPhaseSPML.control(EMmaxit = 1000, EMtol = 1e-05, trace = FALSE)

Arguments

Details

The function implements the semiparametric maximum likelihood (SPML)

Value

A list with components named as the arguments.