[Genabel-commits] r1332 - pkg/PredictABEL/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 6 16:31:43 CEST 2013
Author: lckarssen
Date: 2013-09-06 16:31:39 +0200 (Fri, 06 Sep 2013)
New Revision: 1332
Modified:
pkg/PredictABEL/R/PredictABEL.R
Log:
In PredictABEL: minor change in the layout of the help text.
Modified: pkg/PredictABEL/R/PredictABEL.R
===================================================================
--- pkg/PredictABEL/R/PredictABEL.R 2013-09-06 13:22:24 UTC (rev 1331)
+++ pkg/PredictABEL/R/PredictABEL.R 2013-09-06 14:31:39 UTC (rev 1332)
@@ -12,13 +12,14 @@
#' discrimination box plot, predictiveness curve, and several risk distributions.
#'
#'
-#' These functions can be applied to predicted risks that are obtained using
-#' logistic regression analysis, to weighted or unweighted risk scores, for
-#' which the functions are included in this package. The functions can also be
-#' used to assess risks or risk scores that are constructed using other methods, e.g., Cox Proportional
-#' Hazards regression analysis, which are not included in the current version.
-#' Risks obtained from other methods can be imported into R for assessment
-#' of the predictive performance.
+#' These functions can be applied to predicted risks that are obtained
+#' using logistic regression analysis, to weighted or unweighted risk
+#' scores, for which the functions are included in this package. The
+#' functions can also be used to assess risks or risk scores that are
+#' constructed using other methods, e.g., Cox Proportional Hazards
+#' regression analysis, which are not included in the current version.
+#' Risks obtained from other methods can be imported into R for
+#' assessment of the predictive performance.
#'
#'
#' The functions to construct the risk models using logistic regression analyses
@@ -474,7 +475,7 @@
n1 <- matrix (nrow=dim(p)[2], ncol=12)
for (i in 1:dim(p)[2])
- {
+ {
s<-table(p[,i],o)
if (dim(s)[1]==1) {s <- rbind(s,c(0,0));s <- rbind(s,c(0,0))}
if (dim(s)[1]==2) {s <- rbind(s,c(0,0))}
@@ -490,40 +491,40 @@
d<-oddsratio.wald(c)$measure
e<-oddsratio.wald(c)$data
- m1[i,1]<- colnames(p)[i]
- m1[i,2]<- ( b[1,2])
- m1[i,3]<- ( round((b[1,2]/b[4,2])*100 ,1))
- m1[i,4]<- ( b[2,2])
- m1[i,5]<- ( round((b[2,2]/b[4,2])*100 ,1))
- m1[i,6]<- ( b[3,2])
- m1[i,7]<- ( round((b[3,2]/b[4,2])*100 ,1))
- m1[i,8]<- ( b[1,1])
- m1[i,9]<- ( round((b[1,1]/b[4,1])*100 ,1))
- m1[i,10]<- ( b[2,1])
- m1[i,11]<- ( round((b[2,1]/b[4,1])*100 ,1))
- m1[i,12]<- ( b[3,1])
- m1[i,13]<- ( round((b[3,1]/b[4,1])*100 ,1))
- m1[i,14]<- ( round(a[2,1] ,2))
- m1[i,15]<- ( round(a[2,2] ,2))
- m1[i,16]<- ( round(a[2,3] ,2))
- m1[i,17]<- ( round(a[3,1] ,2))
- m1[i,18]<- ( round(a[3,2] ,2))
- m1[i,19]<- ( round(a[3,3] ,2))
+ m1[i,1]<- colnames(p)[i]
+ m1[i,2]<- ( b[1,2])
+ m1[i,3]<- ( round((b[1,2]/b[4,2])*100 ,1))
+ m1[i,4]<- ( b[2,2])
+ m1[i,5]<- ( round((b[2,2]/b[4,2])*100 ,1))
+ m1[i,6]<- ( b[3,2])
+ m1[i,7]<- ( round((b[3,2]/b[4,2])*100 ,1))
+ m1[i,8]<- ( b[1,1])
+ m1[i,9]<- ( round((b[1,1]/b[4,1])*100 ,1))
+ m1[i,10]<- ( b[2,1])
+ m1[i,11]<- ( round((b[2,1]/b[4,1])*100 ,1))
+ m1[i,12]<- ( b[3,1])
+ m1[i,13]<- ( round((b[3,1]/b[4,1])*100 ,1))
+ m1[i,14]<- ( round(a[2,1] ,2))
+ m1[i,15]<- ( round(a[2,2] ,2))
+ m1[i,16]<- ( round(a[2,3] ,2))
+ m1[i,17]<- ( round(a[3,1] ,2))
+ m1[i,18]<- ( round(a[3,2] ,2))
+ m1[i,19]<- ( round(a[3,3] ,2))
- n1[i,1]<- colnames(p)[i]
- n1[i,2]<- (e[1,2])
- n1[i,3]<- round((e[1,2]/e[3,2])*100,1)
- n1[i,4]<- e[2,2]
- n1[i,5]<- round((e[2,2]/e[3,2])*100,1)
- n1[i,6]<- e[1,1]
- n1[i,7]<- round((e[1,1]/e[3,1])*100,1)
- n1[i,8]<- e[2,1]
- n1[i,9]<- round((e[2,1]/e[3,1])*100,1)
- n1[i,10]<- round(d[2,1] ,2)
- n1[i,11]<- round(d[2,2] ,2)
- n1[i,12]<- round(d[2,3] ,2)
+ n1[i,1]<- colnames(p)[i]
+ n1[i,2]<- (e[1,2])
+ n1[i,3]<- round((e[1,2]/e[3,2])*100,1)
+ n1[i,4]<- e[2,2]
+ n1[i,5]<- round((e[2,2]/e[3,2])*100,1)
+ n1[i,6]<- e[1,1]
+ n1[i,7]<- round((e[1,1]/e[3,1])*100,1)
+ n1[i,8]<- e[2,1]
+ n1[i,9]<- round((e[2,1]/e[3,1])*100,1)
+ n1[i,10]<- round(d[2,1] ,2)
+ n1[i,11]<- round(d[2,2] ,2)
+ n1[i,12]<- round(d[2,3] ,2)
- }
+ }
m1 <- as.table(m1)
dimnames(m1)[[1]] <- c(1:dim(p)[2])
@@ -532,7 +533,7 @@
"CI-Low","CI-high")
names(dimnames(m1)) <- c("", " Genotype frequencies for Cases & Controls, and OR(95% CI)")
if (!missing(filenameGeno))
- write.table(m1,file=filenameGeno, row.names = FALSE,sep = "\t")
+ write.table(m1,file=filenameGeno, row.names = FALSE,sep = "\t")
n1 <- as.table(n1)
dimnames(n1)[[1]] <- c(1:dim(p)[2])
@@ -543,7 +544,7 @@
write.table(n1,file=filenameAllele, row.names = FALSE,sep = "\t")
p<- list(Genotype =m1, Allelic=n1)
- return(p)
+ return(p)
}
#' Function to obtain multivariate odds ratios from a logistic regression model.
#' The function estimates multivariate (adjusted) odds ratios (ORs) with
@@ -624,7 +625,7 @@
round(Upper.CI ,4),"p-value"=round((sum.coef)[, 4],4))
if (!missing(filename))
- write.table(tab,file=filename, row.names=TRUE,sep = "\t")
+ write.table(tab,file=filename, row.names=TRUE,sep = "\t")
B <- mean((riskModel$y) * (1-predict(riskModel, type="response"))^2 +
(1-riskModel$y) * (predict(riskModel, type="response"))^2)
@@ -735,15 +736,15 @@
abline(b,col = 8)
if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
+ if (!missing(fileplot))
savePlot(filename = fileplot,
- type =plottype,
- device = dev.cur())
+ type =plottype,
+ device = dev.cur())
tab<- cbind(riskScore,predRisk)
tab <- as.table(tab)
dimnames(tab)[[2]] <- c("Risk Score", "Predicted risk ")
if (!missing(filename))
- write.table(tab,file=filename, row.names=TRUE,sep = "\t")
+ write.table(tab,file=filename, row.names=TRUE,sep = "\t")
}
#' Function to plot posterior risks against prior risks.
#'
@@ -843,45 +844,45 @@
risk2 <- posteriorrisk
if( plotAll)
- {
- plot(risk1,risk2,xlab= xlabel, ylab=ylabel,main=plottitle,
- cex.lab=1.2, cex.axis=1.1, las=1,pty='s',xlim=rangeaxis ,
- ylim=rangeaxis,pch=20)
- abline(a=0,b=1, lwd=1,col=8)
- if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
- savePlot(filename = fileplot,
- type =plottype,
- device = dev.cur()
- )
+ {
+ plot(risk1,risk2,xlab= xlabel, ylab=ylabel,main=plottitle,
+ cex.lab=1.2, cex.axis=1.1, las=1,pty='s',xlim=rangeaxis ,
+ ylim=rangeaxis,pch=20)
+ abline(a=0,b=1, lwd=1,col=8)
+ if (missing(plottype)) {plottype<- "jpg"}
+ if (!missing(fileplot))
+ savePlot(filename = fileplot,
+ type =plottype,
+ device = dev.cur()
+ )
}
else
{
- op <- par(mfrow=c(1,2),pty="s" )
- plot(risk1,risk2,xlab= xlabel, ylab=ylabel,
- col = (1-(data[,cOutcome]))*1,cex.lab=1.2, cex.axis=1, las=1,pty='s',
- xlim=rangeaxis , ylim=rangeaxis,pch="*")
- abline(a=0,b=1, lwd=1,col=4)
- title(labels[1],cex.main=1)
+ op <- par(mfrow=c(1,2),pty="s" )
+ plot(risk1,risk2,xlab= xlabel, ylab=ylabel,
+ col = (1-(data[,cOutcome]))*1,cex.lab=1.2, cex.axis=1, las=1,pty='s',
+ xlim=rangeaxis , ylim=rangeaxis,pch="*")
+ abline(a=0,b=1, lwd=1,col=4)
+ title(labels[1],cex.main=1)
- plot(risk1,risk2, xlab= xlabel, ylab=ylabel,
- col =(data[,cOutcome])*1, cex.lab=1.2, cex.axis=1, las=1,pty='s',
- xlim=rangeaxis , ylim=rangeaxis,pch="*")
- abline(a=0,b=1, lwd=1,col=4)
- title(labels[2], cex.main=1)
- par(op)
- if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
- savePlot(filename = fileplot,
- type =plottype,
- device = dev.cur())
+ plot(risk1,risk2, xlab= xlabel, ylab=ylabel,
+ col =(data[,cOutcome])*1, cex.lab=1.2, cex.axis=1, las=1,pty='s',
+ xlim=rangeaxis , ylim=rangeaxis,pch="*")
+ abline(a=0,b=1, lwd=1,col=4)
+ title(labels[2], cex.main=1)
+ par(op)
+ if (missing(plottype)) {plottype<- "jpg"}
+ if (!missing(fileplot))
+ savePlot(filename = fileplot,
+ type =plottype,
+ device = dev.cur())
}
- tab<- cbind(risk1,risk2,data[,cOutcome])
- tab <- as.table(tab)
+ tab<- cbind(risk1,risk2,data[,cOutcome])
+ tab <- as.table(tab)
dimnames(tab)[[2]] <- c("Predicted risk 1","Predicted risk 2", "outcome")
- if (!missing(filename))
- write.table(tab,file=filename, row.names=TRUE,sep = "\t")
+ if (!missing(filename))
+ write.table(tab,file=filename, row.names=TRUE,sep = "\t")
}
#' Function for calibration plot and Hosmer-Lemeshow goodness of fit test.
#' The function produces a calibration plot and provides Hosmer-Lemeshow
@@ -970,8 +971,8 @@
p=predRisk
y=data[,cOutcome]
if (length(unique(y))!=2) {
- stop(" The specified outcome is not a binary variable.\n")
- }
+ stop(" The specified outcome is not a binary variable.\n")
+ }
else{
matres <-matrix(NA,nrow=groups,ncol=5)
@@ -994,13 +995,13 @@
pval<- 1-pchisq(chisqr,df)
lines(x =c(0,1),y=c(0,1))
if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
+ if (!missing(fileplot))
savePlot(filename = fileplot,
- type =plottype,
- device = dev.cur())
+ type =plottype,
+ device = dev.cur())
if (!missing(filename))
- write.table(matres,file=filename, row.names=TRUE,sep = "\t",dec=",")
+ write.table(matres,file=filename, row.names=TRUE,sep = "\t",dec=",")
out <- list(Table_HLtest=matres,Chi_square = round(chisqr,3), df=df,
p_value =round(pval,4))
return(out)
@@ -1085,8 +1086,8 @@
if (missing(ylabel)) {ylabel<- "Sensitivity"}
if (class(predrisk) == "numeric") {predrisk<- cbind(predrisk)}
a<-c(1:dim(predrisk)[2])
- for(i in 1:dim(predrisk)[2])
- {
+ for(i in 1:dim(predrisk)[2])
+ {
rAllele <- rcorr.cens(predrisk[,i], data[,cOutcome], outx=FALSE)
pred <- prediction(predrisk[,i], data[,cOutcome])
perf <- performance(pred,"tpr","fpr")
@@ -1107,16 +1108,16 @@
cat("AUC [95% CI] for the model",i, ": ", round(rAllele[1],3),
"[", round(rAllele[1]-1.96/2*rAllele[3],3)," - ",
round(rAllele[1]+1.96/2*rAllele[3],3), "] \n")
- }
- if (!missing(labels)){
- legend( "bottomright",legend= labels, col=c(17:(16+dim(predrisk)[2])),
+ }
+ if (!missing(labels)){
+ legend( "bottomright",legend= labels, col=c(17:(16+dim(predrisk)[2])),
lty=c(1:(dim(predrisk)[2])),lwd =2,cex=1)
}
if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
+ if (!missing(fileplot))
savePlot(filename = fileplot,
- type =plottype,
- device = dev.cur())
+ type =plottype,
+ device = dev.cur())
}
#' Function for predictiveness curve. The function creates a plot of cumulative percentage
@@ -1207,7 +1208,7 @@
n <- length(x)
y <- (1:n)/n
z <- y >= ylim[1] & y <= ylim[2]
-
+
evalCall(plot,
argu = list(x = y[z], y = x[z], type = "n",
xlab = "", ylab = "", las = 1, mgp = c(0, 0.6, 0),
@@ -1215,7 +1216,7 @@
ylim = rangeyaxis, main = plottitle),
checkdef = TRUE,
checkpar = TRUE)
-
+
evalCall(lines,
argu = list(x = y[z], y = x[z], col=16+i, lty=i, lwd=2),
checkdef = TRUE,
@@ -1232,7 +1233,7 @@
n <- length(x)
y <- (1:n)/n
z <- y >= ylim[1] & y <= ylim[2]
-
+
evalCall(lines,
argu = list(x = y[z], y = x[z], col=16+i,lty=i,lwd=2,
add= TRUE), checkdef = TRUE, checkpar = TRUE)
@@ -1347,9 +1348,9 @@
b<-round(seq(rangexaxis[1],rangexaxis[2],2*interval),2)
for(i in 1:length(a))
- {
- if(any(a[i]==b)){a[i]<-a[i] }
- else{a[i] <- ""}
+ {
+ if(any(a[i]==b)){a[i]<-a[i] }
+ else{a[i] <- ""}
}
dimnames(p)[[1]] <- a
@@ -1359,10 +1360,10 @@
cex.lab = 1.2,cex.axis=1.1, las=1,space = c(0,.10))
if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
+ if (!missing(fileplot))
savePlot(filename = fileplot,
- type =plottype,
- device = dev.cur())
+ type =plottype,
+ device = dev.cur())
}
#' Function for reclassification table and statistics.
@@ -1474,8 +1475,8 @@
x<-improveProb(x1=as.numeric(c11)*(1/(length(levels(c11)))),
x2=as.numeric(c22)*(1/(length(levels(c22)))), y=data[,cOutcome])
-
+
y<-improveProb(x1=predrisk1, x2=predrisk2, y=data[,cOutcome])
cat("\n NRI(Categorical) [95% CI]:", round(x$nri,4),"[",round(x$nri-1.96*x$se.nri,4),"-",
@@ -1570,7 +1571,7 @@
mean(risk[data[,cOutcome]==0]),3))
if (missing(plottype)) {plottype<- "jpg"}
- if (!missing(fileplot))
+ if (!missing(fileplot))
savePlot(filename = fileplot,type =plottype,device = dev.cur())
return(p)
}
@@ -1633,7 +1634,7 @@
out<- list(riskModel1=riskmodel1, riskModel2=riskmodel2)
return(out)
}
-#' Function to construct a simulated dataset containing individual genotype data,
+#' Function to construct a simulated dataset containing individual genotype data,
#' genetic risks and disease status for a hypothetical population.
#' Construct a dataset that contains individual genotype data, genetic risk,
#' and disease status for a hypothetical population.
@@ -1735,10 +1736,10 @@
#' van Duijn CM. Predictive testing for complex diseases using multiple genes:
#' fact or fiction? Genet Med. 2006;8:395-400.
#'
-#' Kundu S, Karssen LC, Janssens AC: Analytical and simulation methods for
-#' estimating the potential predictive ability of genetic profiling: a comparison
+#' Kundu S, Karssen LC, Janssens AC: Analytical and simulation methods for
+#' estimating the potential predictive ability of genetic profiling: a comparison
#' of methods and results. Eur J Hum Genet. 2012 May 30.
-#'
+#'
#' van Zitteren M, van der Net JB, Kundu S, Freedman AN, van Duijn CM,
#' Janssens AC. Genome-based prediction of breast cancer risk in the general
#' population: a modeling study based on meta-analyses of genetic associations.
@@ -1775,7 +1776,7 @@
#' # Obtain the AUC and produce ROC curve
#' plotROC(data=Data, cOutcome=4, predrisk=Data[,3])
#'
-"simulatedDataset" <- function(ORfreq, poprisk, popsize, filename)
+"simulatedDataset" <- function(ORfreq, poprisk, popsize, filename)
{
if (missing(poprisk)) {stop("Population disease risk is not specified")}
if (missing(popsize)) {stop("Total number of individuals is not mentioned")}
@@ -1794,103 +1795,103 @@
}
reconstruct.2x3table <- function(OR1,OR2,p1,p2,d,s){
- a <- 1
- eOR <- 0
- while (eOR<=OR2){
- b <- p2*s*(1-d)
- snew <- s-a-b
- p1new <-p1/(1-p2)
- dnew <- (d-(a/s))/((d-(a/s))+ ((1-d)-b/s))
- c <- (OR1*p1new*snew*(1-dnew)*dnew*snew)/((1-p1new)*snew*(1-dnew)+OR1*p1new*snew*(1-dnew))
- dd <- p1new*((1-d)-b/s)*s
- e <- (d-(a/s))*s-c
- f <- ((1-d)-b/s)*s-dd
- eOR <- (a*f)/(b*e)
- tabel <- cbind(a,b,c,dd,e,f,g,OR1,OR2)
- a <- a+1
- tabel
- }
- tabel
+ a <- 1
+ eOR <- 0
+ while (eOR<=OR2){
+ b <- p2*s*(1-d)
+ snew <- s-a-b
+ p1new <-p1/(1-p2)
+ dnew <- (d-(a/s))/((d-(a/s))+ ((1-d)-b/s))
+ c <- (OR1*p1new*snew*(1-dnew)*dnew*snew)/((1-p1new)*snew*(1-dnew)+OR1*p1new*snew*(1-dnew))
+ dd <- p1new*((1-d)-b/s)*s
+ e <- (d-(a/s))*s-c
+ f <- ((1-d)-b/s)*s-dd
+ eOR <- (a*f)/(b*e)
+ tabel <- cbind(a,b,c,dd,e,f,g,OR1,OR2)
+ a <- a+1
+ tabel
+ }
+ tabel
}
reconstruct.2x3tableHWE <- function(OR,p,d,s){
OR1 <- OR
OR2 <- OR^2
- p1 <- 2*p*(1-p)
- p2 <- p*p
+ p1 <- 2*p*(1-p)
+ p2 <- p*p
- a <- 1
- eOR <- 0
- while (eOR<=OR2){
- b <- p2*s*(1-d)
- snew <- s-a-b
- p1new <-p1/(1-p2)
- dnew <- (d-(a/s))/((d-(a/s))+ ((1-d)-b/s))
- c <- (OR1*p1new*snew*(1-dnew)*dnew*snew)/((1-p1new)*snew*(1-dnew)+OR1*p1new*snew*(1-dnew))
- dd <- p1new*((1-d)-b/s)*s
- e <- (d-(a/s))*s-c
- f <- ((1-d)-b/s)*s-dd
- eOR <- (a*f)/(b*e)
- tabel <- cbind(a,b,c,dd,e,f,g,OR1,OR2)
- a <- a+1
- tabel
- }
- tabel
+ a <- 1
+ eOR <- 0
+ while (eOR<=OR2){
+ b <- p2*s*(1-d)
+ snew <- s-a-b
+ p1new <-p1/(1-p2)
+ dnew <- (d-(a/s))/((d-(a/s))+ ((1-d)-b/s))
+ c <- (OR1*p1new*snew*(1-dnew)*dnew*snew)/((1-p1new)*snew*(1-dnew)+OR1*p1new*snew*(1-dnew))
+ dd <- p1new*((1-d)-b/s)*s
+ e <- (d-(a/s))*s-c
+ f <- ((1-d)-b/s)*s-dd
+ eOR <- (a*f)/(b*e)
+ tabel <- cbind(a,b,c,dd,e,f,g,OR1,OR2)
+ a <- a+1
+ tabel
+ }
+ tabel
}
-adjust.postp <- function (pd, LR){
- odds.diff <- 0
- prior.odds <- pd/(1-pd)
- for (i in (1:100000)) {
- Postp <- (prior.odds*LR)/(1+(prior.odds*LR))
- odds.diff <- (pd-mean(Postp))/ (1-(pd-mean(Postp)))
- prior.odds <- prior.odds+odds.diff
- if (odds.diff < .0001) break
- }
- Postp
+adjust.postp <- function (pd, LR){
+ odds.diff <- 0
+ prior.odds <- pd/(1-pd)
+ for (i in (1:100000)) {
+ Postp <- (prior.odds*LR)/(1+(prior.odds*LR))
+ odds.diff <- (pd-mean(Postp))/ (1-(pd-mean(Postp)))
+ prior.odds <- prior.odds+odds.diff
+ if (odds.diff < .0001) break
+ }
+ Postp
}
func.data <- function(p,d,OR,s,g){
Data <- matrix (NA,s,4+g)
- Data[,1] <- rep(0,s)
- Data[,2] <- rep(1,s)
- Data[,3] <- rep(0,s)
- i <- 0
- while (i < g){
+ Data[,1] <- rep(0,s)
+ Data[,2] <- rep(1,s)
+ Data[,3] <- rep(0,s)
+ i <- 0
+ while (i < g){
i <- i+1
cells2x3 <- rep(NA,9)
cells2x3 <- if(p[i,2]==0) {reconstruct.2x2table(p=p[i,1],d,OR=OR[i,1],s)} else {if(p[i,2]==1) {reconstruct.2x3tableHWE(OR=OR[i,1],p=p[i,1],d,s)}
- else {reconstruct.2x3table(OR1=OR[i,1],OR2=OR[i,2],p1=p[i,1],p2=p[i,2],d,s)}}
- LREE <- ((cells2x3[1]/d*s)/(cells2x3[2]/(1-d)*s))
- LREe <- ((cells2x3[3]/d*s)/(cells2x3[4]/(1-d)*s))
- LRee <- ((cells2x3[5]/d*s)/(cells2x3[6]/(1-d)*s))
-
+ else {reconstruct.2x3table(OR1=OR[i,1],OR2=OR[i,2],p1=p[i,1],p2=p[i,2],d,s)}}
+ LREE <- ((cells2x3[1]/d*s)/(cells2x3[2]/(1-d)*s))
+ LREe <- ((cells2x3[3]/d*s)/(cells2x3[4]/(1-d)*s))
+ LRee <- ((cells2x3[5]/d*s)/(cells2x3[6]/(1-d)*s))
+
Gene <- if(p[i,2]==0){c(rep(0,((1-p[i,1]-p[i,2])*s)),rep(1,p[i,1]*s),rep(2,p[i,2]*s))}
else {if(p[i,2]==1) {c(rep(0,(((1-p[i,1])^2)*s)),rep(1,2*p[i,1]*(1-p[i,1])*s),rep(2,p[i,1]*p[i,1]*s))}
- else {c(rep(0,((1-p[i,1]-p[i,2])*s)),rep(1,p[i,1]*s),rep(2,p[i,2]*s))}}
- Filler <- s-length(Gene)
- Gene <- sample(c(Gene,rep(0,Filler)),s,replace=FALSE)
+ else {c(rep(0,((1-p[i,1]-p[i,2])*s)),rep(1,p[i,1]*s),rep(2,p[i,2]*s))}}
+ Filler <- s-length(Gene)
+ Gene <- sample(c(Gene,rep(0,Filler)),s,replace=FALSE)
Data[,4+i] <- Gene
GeneLR <- ifelse(Gene==0,LRee,ifelse(Gene==1,LREe,LREE))
-
+
Data[,1] <- Data[,1]+Gene
- Data[,2] <- Data[,2]*GeneLR
-
-# cat(i,"")
- }
-
- Data[,3] <- adjust.postp(pd=d, LR=Data[,2])
- Data[,4] <- ifelse(runif(s)<=(Data[,3]), 1, 0)
+ Data[,2] <- Data[,2]*GeneLR
+
+# cat(i,"")
+ }
+
+ Data[,3] <- adjust.postp(pd=d, LR=Data[,2])
+ Data[,4] <- ifelse(runif(s)<=(Data[,3]), 1, 0)
Data <- as.data.frame(Data)
Data
}
-
- simulatedData <- func.data (p=ORfreq[,c(3,4)],d=poprisk,OR=ORfreq[,c(1,2)],s=popsize,g=nrow(ORfreq))
+ simulatedData <- func.data (p=ORfreq[,c(3,4)],d=poprisk,OR=ORfreq[,c(1,2)],s=popsize,g=nrow(ORfreq))
-if (!missing(filename))
- {write.table( simulatedData,file=filename, row.names=TRUE,sep = "\t") }
+if (!missing(filename))
+ {write.table( simulatedData,file=filename, row.names=TRUE,sep = "\t") }
+
return(simulatedData)
}
More information about the Genabel-commits
mailing list