Last data update: 2014.03.03
R: EM-PAVA function
EM_PAVA_Func R Documentation
EM-PAVA function
Description
This function is used to estimate the genotype-specific distribution of time-to-event outcomes using EM-PAVA algorithm (Qin et al. 2014).
Usage
EM_PAVA_Func (q, x, delta, timeval, p, ep = 1e-05, maxiter = 400)
Arguments
q
matrix of 2 columns, where the first and second columns are the probabilities of belonging to the carrier p0G
and non-carrier groups 1 - p0G
, respectively.
x
observed event time or censoring time.
delta
indicator of event.
timeval
grid points at which the distribution function values are estimated.
p
number of groups.
ep
convergence criterion. Here, ep = 1e-5
is used as the default value.
maxiter
maximum number of EM iterations.
Details
Technical details can be found in Qin et al. (2014).
Value
Returns a list of prediction values for classes
Fest
estimated values at the points of timeval
.
Fest.all
estimated values of cumulative distribution function on both carrier and non-carrier groups.
References
Qin, J., Garcia, T., Ma, Y., Tang, M., Marder, K. & Wang, Y. (2014). Combining isotonic regression and EM algorithm to predict genetic risk under monotonicity constraint. The Annals of Applied Statistics 8 (2), 1182-1208.
See Also
p0G_Func()
, Sieve_NPMLE_Switch()
Examples
data("Simulated_data");
OY = Simulated_data[,2];
ind = order(OY);
ODelta = Simulated_data[,3];
Op0G = Simulated_data[,4];
Y = OY[ind];
Delta = ODelta[ind];
p0G = Op0G[ind];
Grid = seq(0.01, 3.65, 0.01);
fix_t1 = c(0.288, 0.693, 1.390);
fix_t2 = c(0.779, 1.860, 3.650);
EMpava_result = EM_PAVA_Func ( q = rbind(p0G,1-p0G), x = Y, delta = Delta,
timeval = Grid, p = 2, ep = 1e-4 );
all = sort(c(Grid, Y));
F_carr_func = function(x){ EMpava_result$Fest.all[1, which.max(all[all <= x]) ] };
F_non_func = function(x){ EMpava_result$Fest.all[2, which.max(all[all <= x]) ] };
PAVA_F1.hat_fix_t = apply( matrix(fix_t1, ncol=1), 1, F_carr_func );
PAVA_F2.hat_fix_t = apply( matrix(fix_t2, ncol=1), 1, F_non_func );
PAVA_F.hat_fix_t = data.frame( fix_t1 = fix_t1, PAVA_F1.hat = PAVA_F1.hat_fix_t,
fix_t2 = fix_t2, PAVA_F2.hat = PAVA_F2.hat_fix_t );
print(PAVA_F.hat_fix_t);
# plot estimated curves
F_carr = apply( matrix(Grid, ncol=1), 1, F_carr_func );
F_non = apply( matrix(Grid, ncol=1), 1, F_non_func );
plot( Grid, F_carr, type = 's', lty = 1,
xlab = "Y", ylab = "Estimated Cumulative Distribution Function",
ylim = c(0,1), col = 'blue' );
lines(Grid, F_non, type='s', lty=2, col='red');
legend("topleft", legend=c("Carrier group", "Non-Carrier group"),
lty=c(1,2), col=c("blue", "red") );
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(GSSE)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GSSE/EM_PAVA_Func.Rd_%03d_medium.png", width=480, height=480)
> ### Name: EM_PAVA_Func
> ### Title: EM-PAVA function
> ### Aliases: EM_PAVA_Func
>
> ### ** Examples
>
>
> data("Simulated_data");
>
> OY = Simulated_data[,2];
> ind = order(OY);
> ODelta = Simulated_data[,3];
> Op0G = Simulated_data[,4];
>
> Y = OY[ind];
> Delta = ODelta[ind];
> p0G = Op0G[ind];
>
> Grid = seq(0.01, 3.65, 0.01);
> fix_t1 = c(0.288, 0.693, 1.390);
> fix_t2 = c(0.779, 1.860, 3.650);
>
> EMpava_result = EM_PAVA_Func ( q = rbind(p0G,1-p0G), x = Y, delta = Delta,
+ timeval = Grid, p = 2, ep = 1e-4 );
>
> all = sort(c(Grid, Y));
>
> F_carr_func = function(x){ EMpava_result$Fest.all[1, which.max(all[all <= x]) ] };
> F_non_func = function(x){ EMpava_result$Fest.all[2, which.max(all[all <= x]) ] };
>
> PAVA_F1.hat_fix_t = apply( matrix(fix_t1, ncol=1), 1, F_carr_func );
> PAVA_F2.hat_fix_t = apply( matrix(fix_t2, ncol=1), 1, F_non_func );
>
> PAVA_F.hat_fix_t = data.frame( fix_t1 = fix_t1, PAVA_F1.hat = PAVA_F1.hat_fix_t,
+ fix_t2 = fix_t2, PAVA_F2.hat = PAVA_F2.hat_fix_t );
>
> print(PAVA_F.hat_fix_t);
fix_t1 PAVA_F1.hat fix_t2 PAVA_F2.hat
1 0.288 0.2696055 0.779 0.2488954
2 0.693 0.5084504 1.860 0.4177679
3 1.390 0.8137684 3.650 0.6965440
>
> # plot estimated curves
>
> F_carr = apply( matrix(Grid, ncol=1), 1, F_carr_func );
> F_non = apply( matrix(Grid, ncol=1), 1, F_non_func );
>
> plot( Grid, F_carr, type = 's', lty = 1,
+ xlab = "Y", ylab = "Estimated Cumulative Distribution Function",
+ ylim = c(0,1), col = 'blue' );
> lines(Grid, F_non, type='s', lty=2, col='red');
> legend("topleft", legend=c("Carrier group", "Non-Carrier group"),
+ lty=c(1,2), col=c("blue", "red") );
>
>
>
>
>
>
> dev.off()
null device
1
>