.Rhistory

for(i in nsimul){
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val
}
return(chisq.pvals)
}
pwr.curve(1.5, 0.2, 500, 1000, 1000)
pwr.curve <- function(OR, p, ncontrol, ncase, nsimul){
p.cases <- OR*p #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul)
for(i in nsimul){
print(i)
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val
}
return(chisq.pvals)
}
pwr.curve(1.5, 0.2, 500, 1000, 1000)
pwr.curve <- function(OR, p, ncontrol, ncase, nsimul){
p.cases <- OR*p #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul)
for(i in 1:nsimul){
print(i)
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val
}
return(chisq.pvals)
}
pwr.curve(1.5, 0.2, 500, 1000, 1000)
for(i in 1000){
print(i)
}
for(i in 1:1000){
print(i)
}
pwr.curve <- function(OR, p, ncontrol, ncase, nsimul){
p.cases <- OR*p #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul)
for(i in 1:nsimul){
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val
}
return(chisq.pvals)
}
pwr.curve(1.5, 0.2, 500, 1000, 1000)
pwr.curve(1.5, 0.2, 500, 1000, 1000)
0.1:1
0.1:0.3
1:10
ORs <- c(0:2, by=0.1)
ORs
ORs <- seq(0, 2, 0.1)
ORs
ORs <- seq(0.1, 2, 0.1)
ORs
p.frqs <- seq(0, 1, 0.1)
p.frqs
p.frqs <- seq(0.1, 0.9, 0.05)
p.frqs
p.frqs <- seq(0.05, 0.9, 0.05)
p.frqs
p.frqs <- seq(0.05, 0.95, 0.05)
p.frqs
pwr.curve(1.5, 0.2, 500, 1000, 1000)
test1 <- pwr.curve(1.5, 0.2, 500, 1000, 1000)
test1<0.05
.05/500000
test1<1e-07
sum(test1<1e-07)
sum(test1<=1e-07)
sum(test1<=1e-07)/length(test1)
test1 <- pwr.curve(.5, 0.2, 500, 1000, 1000)
sum(test1<=0.05)
test1
test1 <- pwr.curve(.05, 0.2, 500, 1000, 1000)
sum(test1<=0.05)
test1
test1 <- pwr.curve(.005, 0.2, 500, 1000, 1000)
test1
test1 <- pwr.curve(.2, 0.2, 500, 1000, 1000)
test1
test1 <- pwr.curve(.5, 0.2, 500, 1000, 1000)
test1
.01*.2
pwr.calc <- function(OR, p, ncontrol, ncase=1000, nsimul=1000){
p.cases <- (OR*(p/(1-p)))/(1+(OR*(p/(1-p)))) #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul) #Initializes a blank vector for storing the chi-squared test p-values on the simulations of case and control data seen below.
for(i in 1:nsimul){ #Iterates through case and control data creation
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val #Extracts the p-value from the chisq.test result and puts it into the chisq.pvals vector.
}
return(chisq.pvals)
}
individual.OR.table <- function(OR){
p.frqs <- seq(0.05, 0.95, 0.05)
sample.sizes <- c(500, 1000, 1500, 2000, 3000)
final.table <- matrix(nrow=length(p.frqs)*length(sample.sizes), ncol=3)
colnames(final.table) <- c("Allele frequency", "Sample Size", "Power")
index <- 1
for(frequency in p.frqs){
for(size in sample.sizes){
power <- sum(pwr.calc(OR, frequency, size)<=(0.05/500000))
final.table[index,] <- c(frequency, size, power)
index <- index+1
}
}
return(final.table)
}
power.plotter <- function(OR){
data <- individual.OR.table(OR)
ordered.data <- data[order(data[,2]),]
plot(ordered.data[1:19,1], ordered.data[1:19, 3], xlab="Allele frequency in population", ylab="Power", ylim=c(0, 100), type="l", col="red", main=paste("Power at OR ", OR))
lines(ordered.data[20:38, 1], ordered.data[20:38, 3], col="green")
lines(ordered.data[39:57, 1], ordered.data[39:57, 3], col="blue")
lines(ordered.data[58:76, 1], ordered.data[58:76, 3], col="orange")
lines(ordered.data[77:95, 1], ordered.data[77:95, 3], col="purple")
legend(0.83, 99, legend=unique(ordered.data[,2]), title="Controls", title.col="black", text.col=c("red", "green", "blue", "orange", "purple"), cex=0.75)
}
for(OR in ORs){
power.plotter(OR)
}
pwr.calc <- function(OR, p, ncontrol, ncase=1000, nsimul=1000){
p.cases <- (OR*(p/(1-p)))/(1+(OR*(p/(1-p)))) #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul) #Initializes a blank vector for storing the chi-squared test p-values on the simulations of case and control data seen below.
for(i in 1:nsimul){ #Iterates through case and control data creation
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val #Extracts the p-value from the chisq.test result and puts it into the chisq.pvals vector.
}
return(chisq.pvals)
}
individual.OR.table <- function(OR){
p.frqs <- seq(0.05, 0.95, 0.05)
sample.sizes <- c(500, 1000, 1500, 2000, 3000)
final.table <- matrix(nrow=length(p.frqs)*length(sample.sizes), ncol=3)
colnames(final.table) <- c("Allele frequency", "Sample Size", "Power")
index <- 1
for(frequency in p.frqs){
for(size in sample.sizes){
power <- sum(pwr.calc(OR, frequency, size)<=(0.05/500000))
final.table[index,] <- c(frequency, size, power)
index <- index+1
}
}
return(final.table)
}
power.plotter <- function(OR){
data <- individual.OR.table(OR)
ordered.data <- data[order(data[,2]),]
plot(ordered.data[1:19,1], ordered.data[1:19, 3], xlab="Allele frequency in population", ylab="Power", ylim=c(0, 100), type="l", col="red", main=paste("Power at OR ", OR))
lines(ordered.data[20:38, 1], ordered.data[20:38, 3], col="green")
lines(ordered.data[39:57, 1], ordered.data[39:57, 3], col="blue")
lines(ordered.data[58:76, 1], ordered.data[58:76, 3], col="orange")
lines(ordered.data[77:95, 1], ordered.data[77:95, 3], col="purple")
legend(0.83, 99, legend=unique(ordered.data[,2]), title="Controls", title.col="black", text.col=c("red", "green", "blue", "orange", "purple"), cex=0.75)
}
ORs <- seq(0.1, 2, 0.1)
p.frqs <- seq(0.05, 0.95, 0.05)
sample.sizes <- c(500, 1000, 1500, 2000, 3000)
for(OR in ORs){
power.plotter(OR)
}
individual.OR.table(1.5)
?expand.grid()
?names()
test1 <- matrix(c(50, 100, 60, 120), nrow=2, ncol=2)
test1
test2 <- 2*test1
test2
chisq.test(test1)
test1 <- matrix(c(50, 100, 300, 475), nrow=2, ncol=2)
test1
chisq.test(test1)
test2 <- 2*test1
chisq.test(test2)
test1
test2
dis.table <- cbind(c(0.1, 0.2, 0.3, 0.4, 0.5), c(4, 2, 1, 3, 5), c(10, 50, 20, 30, 40))
dis.table
ordered.table <- dis.table[order(dis.table[,3], dis.table[,2]),]
ordered.table
dis.table <- cbind(c(0.1, 0.2, 0.3, 0.4, 0.5), c(4, 2, 1, 3, 5), c(10, 20, 20, 10, 20))
dis.table
ordered.table <- dis.table[order(dis.table[,3], dis.table[,2]),]
ordered.table
pwr.calc <- function(OR, p, ncontrol, ncase=1000, nsimul=1000){
p.cases <- (OR*(p/(1-p)))/(1+(OR*(p/(1-p)))) #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul) #Initializes a blank vector for storing the chi-squared test p-values on the simulations of case and control data seen below.
for(i in 1:nsimul){ #Iterates through case and control data creation for nsimul number of times.
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val #Extracts the p-value from the chisq.test result and puts it into the chisq.pvals vector.
}
return((sum(chisq.pvals<=(0.05/500000)))/length(chisq.pvals))
}
individual.OR.table <- function(OR){
p.frqs <- seq(0.05, 0.95, 0.05)
sample.sizes <- c(500, 1000, 1500, 2000, 3000)
final.table <- matrix(nrow=length(p.frqs)*length(sample.sizes), ncol=3)
colnames(final.table) <- c("Allele frequency", "Sample Size", "Power")
index <- 1
for(frequency in p.frqs){
for(size in sample.sizes){
power <- pwr.calc(OR, frequency, size)
final.table[index,] <- c(frequency, size, power)
index <- index+1
}
}
return(final.table)
}
power.plotter <- function(OR){
data <- individual.OR.table(OR)
ordered.data <- data[order(data[,2]),]
plot(ordered.data[1:19,1], ordered.data[1:19, 3], xlab="Allele frequency in population", ylab="Power", ylim=c(0, 100), type="l", col="red", main=paste("Power at OR ", OR))
lines(ordered.data[20:38, 1], ordered.data[20:38, 3], col="green")
lines(ordered.data[39:57, 1], ordered.data[39:57, 3], col="blue")
lines(ordered.data[58:76, 1], ordered.data[58:76, 3], col="orange")
lines(ordered.data[77:95, 1], ordered.data[77:95, 3], col="purple")
legend(0.83, 99, legend=unique(ordered.data[,2]), title="Controls", title.col="black", text.col=c("red", "green", "blue", "orange", "purple"), cex=0.75)
}
ORs <- seq(0.1, 2, 0.1)
p.frqs <- seq(0.05, 0.95, 0.05)
sample.sizes <- c(500, 1000, 1500, 2000, 3000)
for(OR in ORs){
power.plotter(OR)
}
power.plotter(1.5)
pwr.calc <- function(OR, p, ncontrol, ncase=1000, nsimul=100){
p.cases <- (OR*(p/(1-p)))/(1+(OR*(p/(1-p)))) #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul) #Initializes a blank vector for storing the chi-squared test p-values on the simulations of case and control data seen below.
for(i in 1:nsimul){ #Iterates through case and control data creation for nsimul number of times.
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val #Extracts the p-value from the chisq.test result and puts it into the chisq.pvals vector.
}
return((sum(chisq.pvals<=(0.05/500000)))/length(chisq.pvals))
}
power.plotter(1.5)
pwr.calc <- function(OR, p, ncontrol, ncase=1000, nsimul=100){
p.cases <- (OR*(p/(1-p)))/(1+(OR*(p/(1-p)))) #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul) #Initializes a blank vector for storing the chi-squared test p-values on the simulations of case and control data seen below.
for(i in 1:nsimul){ #Iterates through case and control data creation for nsimul number of times.
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val #Extracts the p-value from the chisq.test result and puts it into the chisq.pvals vector.
}
print(sum(chisq.pvals<=(0.05/500000)))
dim(chisq.pvals)
print(length(chisq.pvals))
return((sum(chisq.pvals<=(0.05/500000)))/length(chisq.pvals))
}
power.plotter(1.5)
pwr.calc <- function(OR, p, ncontrol, ncase=1000, nsimul=100){
p.cases <- (OR*(p/(1-p)))/(1+(OR*(p/(1-p)))) #Recall that the odds ratio is the probability of having the risk genotype given that you have the disease/probabilty of having the risk genotype given you don't have the disease. OR=p.cases/p.controls, hence p.cases=OR*p.controls
chisq.pvals <- vector(length=nsimul) #Initializes a blank vector for storing the chi-squared test p-values on the simulations of case and control data seen below.
for(i in 1:nsimul){ #Iterates through case and control data creation for nsimul number of times.
control.vec <- rbinom(n=(2*ncontrol), size=1, prob=p) #Creates a vector of controls.
case.vec <- rbinom(n=(2*ncase), size=1, prob=p.cases) #Creates a vector of cases.
testing.table <- matrix(c(sum(case.vec), sum(control.vec), sum(case.vec==0), sum(control.vec==0)), nrow=2, ncol=2) #Creates the 2x2 contingency table for chisq.test to be run on. Note that this has cases in the first row and controls in the second, and those with the risk genotype in the first column and those without in the second column. Not that it matters for the chisquared test, but you could reverse this by adding the byrow=TRUE flag to this call of "matrix()".
chisq.pvals[i] <- chisq.test(testing.table)$p.val #Extracts the p-value from the chisq.test result and puts it into the chisq.pvals vector.
}
print(sum(chisq.pvals<=(0.05/500000)))
print(length(chisq.pvals))
return((sum(chisq.pvals<=(0.05/500000)))/length(chisq.pvals))
}
power.plotter(1.5)
test3 <- individual.OR.table(1.5)
test3
fisher.test(c(.03, .1, .02, .03, .08))
OR.vec <- c(.03, .1, .02, .03, .08)
log.OR.to.pval <- function(OR.vector){
pchisq(OR.vector*(2*log(10)),df=1,lower.tail=FALSE)/2
}
log.OR.to.pval(mean(OR.vec))
pchisq(OR.vector*(2*log(10)),df=10,lower.tail=FALSE)/2
OR.vec <- c(.03, .1, .02, .03, .08)
log.OR.to.pval <- function(OR.vector){
pchisq(OR.vector*(2*log(10)),df=10,lower.tail=FALSE)/2
}
log.OR.to.pval(mean(OR.vec))
log.OR.to.pval(OR.vec)
y <- c(0.03, 0.1, 0.02, 0.03, 0.08)
sd <- c(0.04, 0.025, 0.04, 0.03, 0.025)
p <- 1 - pnorm(y, sd=sd)
p
p <- pnorm(y, sd=sd)
p
p <- pnorm(y, sd=sd, lower.tail=FALSE)
p
p <- 1 - pnorm(y, sd=sd) #Could also just take pnorm with lower.tail=FALSE.
p
1 - pchisq( -2*sum(log(p)), 2*length(p))
pchisq(-2*sum(ln(p)), 10)
pchisq(-2*sum(log(p)), 10, lower.tail=FALSE)
?log()
y <- c(0.03, 0.1, 0.02, 0.03, 0.08)
sd <- c(0.04, 0.025, 0.04, 0.03, 0.025)
# Fisher's method
p <- 1 - pnorm(y, sd=sd) #Could also just take pnorm with lower.tail=FALSE.
1 - pchisq( -2*sum(log(p)), 2*length(p)) #Again, could also just take pchisq with lower.tail=FALSE.
# Stouffer's method
Z <- qnorm(1-p)
n <- 1/sd^2
wz <- sqrt(n)
Z.summary <- sum(wz*Z) / sqrt(sum(wz^2))
p.summary <- 1 - pnorm(Z.summary,1)
Z.summary
p.summary
wz
n
?pnorm()
ORs <- c(0.03, 0.1, 0.02, 0.03, 0.08) #Set up a vector of the log(OR) values from the studies.
SEs <- c(0.04, 0.025, 0.04, 0.03, 0.025) #Set up a vector of the standard errors from the studies.
individual.pvals <- 1 - pnorm(ORs, sd=SEs) #Could also just take pnorm with lower.tail=FALSE.
individual.pvals
?pchisq()
individual.pvals <- 1 - pnorm(ORs, sd=SEs) #Gets individual p-values from each of the studies. Since log(OR) stats are normally distributed, we can find their p-values by inputting them into the normal distribution (with their standard deviations being the standard error for each respective study). Note that you could also just take pnorm with lower.tail=FALSE.
1 - pchisq( -2*sum(log(individual.pvals)), 2*length(individual.pvals)) #Now, with the individual p-vals, we just input them into the equation given for Fisher's method of combining p-values (slide 5 of GWAS meta analysis lecture), and since the test statistic this equation yields is chi-squared distributed with degrees of freedom=2k (k being the number of studies), we can find the overall p-value by looking at the chi-squared distribution with pchisq(). Again, one could also just take pchisq with lower.tail=FALSE, as opposed to subtracting it from 1.
?qnorm()
individual.pvals
Z <- qnorm(individual.pvals)
Z
Z <- qnorm(1-individual.pvals)
Z
Z.scores <- qnorm(1-individual.pvals) #We did 1-individual.pvals just to get the Z-scores to be positive, rather than negative. It mostly matters for the downstream calculations.
Z.scores <- qnorm(1-individual.pvals) #We did 1-individual.pvals just to get the Z-scores to be positive, rather than negative. It mostly matters for the downstream calculations.
weights <- sqrt(sample.sizes)
sample.sizes <- 1/SEs^2
ORs <- c(0.03, 0.1, 0.02, 0.03, 0.08) #Set up a vector of the log(OR) values from the studies.
SEs <- c(0.04, 0.025, 0.04, 0.03, 0.025) #Set up a vector of the standard errors from the studies.
# A) Fisher's method
individual.pvals <- 1 - pnorm(ORs, sd=SEs) #Gets individual p-values from each of the studies. Since log(OR) stats are normally distributed, we can find their p-values by inputting them into the normal distribution (with their standard deviations being the standard error for each respective study). Note that you could also just take pnorm with lower.tail=FALSE.
1 - pchisq( -2*sum(log(individual.pvals)), 2*length(individual.pvals)) #Now, with the individual p-vals, we just input them into the equation given for Fisher's method of combining p-values (slide 5 of GWAS meta analysis lecture), and since the test statistic this equation yields is chi-squared distributed with degrees of freedom=2k (k being the number of studies), we can find the overall p-value by looking at the chi-squared distribution with pchisq(). Again, one could also just take pchisq with lower.tail=FALSE, as opposed to subtracting it from 1.
#The output yields a combined p-value for the significance of the SNP across the 5 studies: ~2.9*10^-6
# B) Stouffer's method
#First, the hint tips us off that we need to find the sample sizes for each of these studies. Slide 10 of the same lecture notes that sampling variance of effect size (which is what our SEs are the square root of, since variance is standard deviation squared) is proportional to the inverse of sample size. Conversely, sample size is inversely proportion to SEs^2 (sampling variance of effect size, i.e. the log ORs). Hence:
sample.sizes <- 1/SEs^2
#We also know we need Z-scores, rather than p-values, in order to use Stouffer's method. We can convert our vector of individual p-values into a vector of Z scores (with the knowledge that Z-scores are normally distributed). In R, the qnorm() function can be used to convert probabilities such as p-values, akin to distribution quantiles, to Z-scores. We know to use this because it's given in slide 6 of Xin's meta-analysis lecture--to convert p-values to Z scores, we take 1-p of the inverse of the normal cumulative distribution function, phi. Hence:
Z.scores <- qnorm(1-individual.pvals) #We did 1-individual.pvals just to get the Z-scores to be positive, rather than negative. It mostly matters for the downstream calculations.
#Looking at the slide for Stouffer's method, we can see that the weight of each study, w, is proportional to the sample size of the study. Specifically, the weights are the square roots of the sample sizes. This, and the way to obtain the Z.scores above, are both nicely reviewed in Box 3 of the Evangelou and Ioannidis review on Chalk. Note that in that box, they divide the p-values by 2, assuming they are two-sided (ours are already one sided).
weights <- sqrt(sample.sizes)
weights*Z.scores
1 - pchisq( -2*sum(log(individual.pvals)), 2*length(individual.pvals)) #Now, with the individual p-vals, we just input them into the equation given for Fisher's method of combining p-values (slide 5 of GWAS meta analysis lecture), and since the test statistic this equation yields is chi-squared distributed with degrees of freedom=2k (k being the number of studies), we can find the overall p-value by looking at the chi-squared distribution with pchisq(). Again, one could also just take pchisq with lower.tail=FALSE, as opposed to subtracting it from 1.
meta.pval <- 1 - pchisq( -2*sum(log(individual.pvals)), 2*length(individual.pvals)) #Now, with the individual p-vals, we just input them into the equation given for Fisher's method of combining p-values (slide 5 of GWAS meta analysis lecture), and since the test statistic this equation yields is chi-squared distributed with degrees of freedom=2k (k being the number of studies), we can find the overall p-value by looking at the chi-squared distribution with pchisq(). Again, one could also just take pchisq with lower.tail=FALSE, as opposed to subtracting it from 1.
meta.pval
Z.summary
p.summary <- 1 - pnorm(Z.summary, 1) #Again, we could've also chosen to do lower.tail=FALSE instead of subtracting this value from 1.
p.summary
pnorm(Z.summary, 1, lower.tail=FALSE)
p.summary
summary.effect <- sum(sample.sizes * ORs) / sum(sample.sizes)
summary.effect
sd.summary.effect <- sqrt(1/sum(sample.sizes))
sd.summary.effect
1-pnorm(summary.effect/sd.summary.effect) #Again, could also have used lower.tail=FALSE instead of 1-this.
sample.sizes
ORs
Cochran.Q
Cochran.Q <- sum(((ORs-summary.effect)/SEs)^2) #Note that the standard errors of the studies are already the square roots of their individual variances, that y is, again, the response variable of interest here (the effect size), and that we utilize the summary effect we obtained from the fixed effects model.
Cochran.Q
1 - pchisq(Cochran.Q, length(y)-1)
τ.squared <- (Cochran.Q - (length(ORs) - 1)) / (sum(sample.sizes) - sum(sample.sizes^2)/sum(sample.sizes))
SEs
weights.RE <- 1/(SEs^2 + τ.squared)
summary.effect.RE <- sum(weights.RE * ORs)/sum(weights.RE)
summary.effect.RE
sd.summary.effect.RE <- sqrt(1/sum(weights.RE))
1-pnorm(summary.effect.FE/sd.summary.effect.FE) #Again, could also have used lower.tail=FALSE instead of 1-this.
summary.effect.FE <- sum(sample.sizes * ORs) / sum(sample.sizes)
summary.effect.FE #And we can see that the summary effect is ~0.063.
#Now, we want to obtain the significance of the summary effect--aka, its p-value. To find this, we need to get the Z-score based on its standard error. For fixed effect models, as per Box 3 from the Evangelou and Ioannidis review, the variance is equal to 1/sum(weights). Again, our weights here are the vector of sample.sizes, and since we are interested in the standard error/deviation, rather than the variance, we would take the square root of this. This should also be clear from slide 14 of Xin's lecture:
sd.summary.effect.FE <- sqrt(1/sum(sample.sizes))
sd.summary.effect.FE #Standard error of the summary effect is roughly 0.013.
#Now, we evaluate the significance of this summary effect. We obtain the Z-score by dividing the point estimate of the mean by its standard deviation (definition of a Z-score), and use pnorm() on that to obtain a p-value. I don't enter an sd into the pnorm() call, once again because it's a Z-score we have here, and those should be normally distributed with sd=1 (the default for the function) under the null hypothesis.
1-pnorm(summary.effect.FE/sd.summary.effect.FE) #Again, could also have used lower.tail=FALSE instead of 1-this.
sd.summary.effect.RE <- sqrt(1/sum(weights.RE))
1 - pnorm(mu.RE/sd.mu.RE)
1 - pnorm(summary.effect.RE/sd.summary.effect.RE)
weights.RE
weights.FE
sample.sizes
sample.size/sum(sample.sizes)
sample.sizes/sum(sample.sizes)
weights.RE
weights.RE/sum(weights.RE)
weights.RE/sum(weights.RE) #Vector of weight proportionalities under RE model.
sample.sizes/sum(sample.sizes) #Vector of weight proportionalities under FE model.
weights.RE/sum(weights.RE) #Vector of weight proportionalities under RE model.
summary.effect.FE
summary.effect.RE
sd.summary.effect.FE
sd.summary.effect.RE
pchisq(Cochran.Q, 4, lower.tail=FALSE)
pchisq(Cochran.Q, 4, lower.tail=FALSE)
Cochran.Q
(Cochran.Q-4)/Cochran.Q
T.squared <- 5
T.squared
T.squared <- (Cochran.Q-4)/(sum(sample.sizes)-(sum(sample.sizes^2)/sum(sample.sizes)))
T.squared
I.squared
I.squared <- ((Cochran.Q-4)/Cochran.Q)*100
I.squared #Roughtl 29%
pnorm()
pnorm(0)
pchisq( -2*sum(log(individual.pvals)), 2*length(individual.pvals))
Z.scores
1-qnorm(individual.pvals)
qnorm(individual.pvals)
individual.pvals
1-pnorm(ORs)
p.summary <- 1 - pnorm(Z.summary, 1) #Again, we could've also chosen to do lower.tail=FALSE instead of subtracting this value from 1.
p.summary #The output yields a combined p-value for the significance of the SNP across the 5 studies: ~9.6*10^-5...slightly less significant than the p-value obtained utilizing Fisher's method. Interesting!
1-pnorm(Z.summart)
1-pnorm(Z.summary)
1-pnorm(Z.summary, 1)
pnorm(Z.summary)
1-pnorm(Z.summary)
Z.summary
sd.summary.effect.FE #Standard error of the summary effect is roughly 0.013.
?pchisq()
Z.scores
qnorm(individual.pvals, lower.tail=FALSE)
qnorm(individual.pvals, lower.tail=TRUE)
Z.summary <- sum(weights*Z.scores) / sqrt(sum(weights^2)) #Note that, from looking at the equation, we see that the denominator is the square root of N, the total sample size of all the studies combined. Hence here we square the weights vector (since it's the sqrt of sample sizes), sum it all together to get the total sample size of all the studies combined, and then take the sqrt of that as per the equation. This gives us a meta Z-score:
Z.summary #Output is 4.728038.
Z.summary <- sum(weights*Z.scores) / sqrt(sum(sample.sizes)) #Note that, from looking at the equation, we see that the denominator is the square root of N, the total sample size of all the studies combined. Hence here we square the weights vector (since it's the sqrt of sample sizes), sum it all together to get the total sample size of all the studies combined, and then take the sqrt of that as per the equation. This gives us a meta Z-score:
Z.summary #Output is 4.728038.
?pnorm()
y <- c(0.03, 0.1, 0.02, 0.03, 0.08)
sd <- c(0.04, 0.025, 0.04, 0.03, 0.025)
# Fisher's method
p <- 1 - pnorm(y, sd=sd) #Could also just take pnorm with lower.tail=FALSE.
1 - pchisq( -2*sum(log(p)), 2*length(p)) #Again, could also just take pchisq with lower.tail=FALSE.
# Stouffer's method
Z <- qnorm(1-p)
n <- 1/sd^2
wz <- sqrt(n)
Z.summary <- sum(wz*Z) / sqrt(sum(wz^2))
p.summary <- 1 - pnorm(Z.summary,1)
p.summary
pnorm(Z.summary, lower.tail=FALSE)
w <- 1/sd^2
mu <- sum(w * y) / sum(w)
sd.mu <- sqrt(1/sum(w))
1-pnorm(mu / sd.mu)
p.summary <- 1 - pnorm(Z.summary) #Again, we could've also chosen to do lower.tail=FALSE instead of subtracting this value from 1.
p.summary #The output yields a combined p-value for the significance of the SNP across the 5 studies: ~9.6*10^-5...slightly less significant than the p-value obtained utilizing Fisher's method. Interesting!
pnorm(summary.effect.FE, summary.effect.FE, sd.summary.effect.FE)
1-pnorm(summary.effect.FE/sd.summary.effect.FE) #Again, could also have used lower.tail=FALSE instead of 1-this.
sqrt(T.squared)
Z.scores
meta.pval #The output yields a combined p-value for the significance of the SNP across the 5 studies: ~2.9*10^-6
p.summary
pnorm(1-Z.summary)
Z.summary
1-pnorm(Z.summary)
p.summary
pnorm(Z.summary, lower.tail=FALSE)
printwd()
getwd()
sys.getenv()
Sys.getenv()
dice
sample_attributes=read.delim("~/Desktop/eQTL/GTEx_Data_V6_Annotations_SampleAttributesDS.txt",sep="\t")
ngenes=1000
test=read.table("~/Desktop/eQTL/GTEx_Analysis_v6_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct.gz",header=TRUE,skip = 2)
library(reshape2)
x <- melt(test)
install.packages("Cairo")
install.packages("shiny")
install.packages("GOTHiC")
library(BiocInstaller)
biocLite()
biocLite("GOTHiC")
library("GOTHiC", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("GOTHiC")
setwd("~/Testing/GILAD/g-i-l-a-d")
setwd("/Users/ittaieres/Testing/GILAD/g-i-l-a-d ")
require(gplots)
install.packages("gplots")
require(gplots)
# Change file name at will:
word.list = as.matrix(read.table("bingo_words_y.txt"))
word.list <- gsub("\\n", "\n", word.list, fixed=T)
nw = dim(word.list)[1]
numbers.m <- matrix (1,nrow=5, ncol=5)
# To control the number of cards printed, change the number in the loop.
for (i in 1:11){
# Edit the following two lines if you want to change the central square, or if you simply want to sample 25 terms.
words.ind = sample(1:nw,24)
words = c(word.list[words.ind[1:12],],"FDR",word.list[words.ind[13:24],])
words.m = matrix(words,nrow=5,ncol=5)
#Insert heatmap here
pdf(paste("Bingo_Card_",i,".pdf", sep=""))
heatmap.2 (numbers.m, Rowv=NA, Colv=NA, col =colorRampPalette(c("white","white"))(5), scale="none", dendrogram="none", trace="none", density.info="none", key=F, cellnote=words.m, notecol="black", na.color="grey90",colsep=c(1,2,3,4,5), rowsep=c(1,2,3,4,5), sepcolor="black", sepwidth=c(0.02,0.02), labRow=F, labCol=F, main="G-I-L-A-D", lwid=c(0.2,1), lhei=c(0.2,1))
dev.off()
}