## author: Hae Kyung Im
## profile likelihood for estimating effect of gene expression on drug effect
## IC50 ~ gene.exp with trio covariance structure
##
## profile restricted likelihood
## input Y (e.g. IC50), M (3m by q matrix)
## m = number of trios
## q number of explanatory variables including constant 1
## rho is the correlation between parent and child
pRL = function(rho, Y, M, n, m, q)
{
## rho2 = rho^2
## iRho3 = matrix( c(1-rho2, -rho, rho2, -rho, 1, -rho, rho2, -rho, 1-rho2),
## 3,3, byrow=T)
## iRho3 = iRho3 / (1-2*rho2)
## iRho = kronecker(diag(m), iRho3)
iRho = iRhofun(rho, m)
logdetRho = m * log(1 - 2*rho^2)
Wrho = t(M) %*% iRho %*% M
iWrho = solve(Wrho)
Qrho = diag(n) - M %*% iWrho %*% t(M) %*% iRho
res = - (n-q)/2 * log( t(Y) %*% iRho %*% Qrho %*% Y ) - 1/2 * logdetRho +
- 1/2 * log( det(Wrho) )
return(res)
}
## restricted likelihood. Use this for calculating standard errors of estimates rhohat and sig2hat
RL = function(rho, sig2, Y, M, n, m, q)
{
iRho = iRhofun(rho, m)
logdetRho = m * log(1 - 2*rho^2)
Wrho = t(M) %*% iRho %*% M
iWrho = solve(Wrho)
Qrho = diag(n) - M %*% iWrho %*% t(M) %*% iRho
res = - (n-q)/2 * log(sig2) - 1/2 * logdetRho - 1/2 * log( det(Wrho) ) +
t(Y) %*% iRho %*% Qrho %*% Y/ (2*sig2)
return(res)
}
## ener function to minimize ener = - loglik
ener = function(rho,Y,M,n,m,q)
{
- pRL(rho,Y,M,n,m,q)
}
## sig2hat maximizes RL for given rho
sig2fun = function(rho, Y, M, n, m, q)
{
## rho2 = rho^2
## iRho3 = matrix( c(1-rho2, -rho, rho2, -rho, 1, -rho, rho2, -rho, 1-rho2),
## 3,3, byrow=T)
## iRho3 = iRho3 / (1-2*rho2)
## iRho = kronecker(diag(m), iRho3)
iRho = iRhofun(rho,m)
Wrho = t(M) %*% iRho %*% M
iWrho = solve(Wrho)
Qrho = diag(n) - M %*% iWrho %*% t(M) %*% iRho
res = t(Y) %*% iRho %*% Qrho %*% Y / (n-q)
return(res)
}
## betahat maximizes likelihood given sig2 and rho
betafun = function(rho, Y, M, n, m, q)
{
iRho = iRhofun(rho, m)
Wrho = t(M) %*% iRho %*% M
iWrho = solve(Wrho)
res = iWrho %*% t(M) %*% iRho %*% Y
}
iRhofun = function(rho,m)
{
rho2 = rho^2
iRho3 = matrix( c(1-rho2, -rho, rho2, -rho, 1, -rho, rho2, -rho, 1-rho2),
3,3, byrow=T)
iRho3 = iRho3 / (1-2*rho2)
res = kronecker(diag(m), iRho3)
return(res)
}
Tuesday, February 9, 2010
simulation of trios and estimation of rho
## simulation of trios and estimation of rho
## test
rm(list=ls())
options(digits=3,digits.sec=4)
source('/Volumes/iDisk/Documents/R_functions/myfunctions2.r')
pre = '/Users/haky/Clients/'
input.dir = pre %&% 'Dolan.Eileen/miscovariance/'
output.dir = input.dir
source.dir = pre %&% 'Dolan.Eileen/miscovariance/'
source(source.dir %&% "functions.r")
setwd(input.dir)
## simulate trios
## Y = mu + X + epsi
## epsi ~ N(0,rho)
newsim = T
if(newsim)
{
n = 300
m = n/3
rho = .6
sig2 = 10
sig = sqrt(sig2)
Rho3 = matrix( c(1, rho, 0, rho, 1, rho, 0, rho, 1),3,3, byrow=T)
Rho = kronecker(diag(m), Rho3)
X = rnorm(n)
Z = rnorm(n)
C = chol(Rho)
epsi = t(C) %*% Z
Y = X + sig * epsi
M = matrix( c(rep(1,n), X), n,NCOL(X)+1, byrow=F )
triono = rep(1:m,each=3)
sortno = rep(1:3, m)
data = data.frame(drug=Y, gene=M[,2], triono, sortno)
write.table(data,file=input.dir %&% "simdata.csv",sep=",",
col.names=T,quote=F, row.names=F)
} else {
data = read.table(file=input.dir %&% "simdata.csv",header=T,sep=",")
}
Y = data$drug
n = NROW(Y)
X = data$gene
M = matrix( c(rep(1,n),X), n,NCOL(X)+1, byrow=F )
m = n/3
q = NCOL(M)
system.time( optimize(ener,c(-1/sqrt(2),1/sqrt(2)),Y, M, n,m,q))
opres = optimize(ener,c(-1/sqrt(2),1/sqrt(2)),Y, M, n,m,q)
rhohat = opres$min
sig2hat = sig2fun(rhohat, Y, M, n, m, q)
betahat = betafun(rhohat, Y, M, n, m, q)
print( c(rhohat, sqrt(sig2hat), betahat ) )
## compare with simple linear regression
fit = lm(Y~M-1)
print(summary(fit)$coef)
print(summary(fit)$sigma)
## install nlme
##install.packages("nlme")
##library(nlme)
##lmefit = lme(drug~gene, data=data, corr = corExp(1,~sortno|triono), random=~1|triono)
## nlme package doesn't allow for toeplitz correlation structure
Tuesday, February 2, 2010
q-values
q-values gives natural way to choose cutoff values. If a small increase in FDR will increase number of SNPs drastically then it makes sense to include more.
Thursday, December 3, 2009
linear regression
Linear regression [1] "" [1] "-----resp1 and resp2" p12019 p13881 pall qall rs4946514 0.00566 0.0183 4.9e-05 0.000945 p12019 p13881 pall qall rs4946514 0.00566 0.0183 0.000049 0.000945 rs1934632 0.33406 0.0163 0.000258 0.001892 rs10863854 0.33704 0.0163 0.000294 0.001892 rs1934626 0.33262 0.0603 0.004931 0.023770 rs8021893 0.36013 0.0748 0.007678 0.027420 rs3813437 0.06367 0.1076 0.009607 0.027420 rs1317255 0.58202 0.0787 0.009955 0.027420 rs1934633 0.58538 0.0787 0.011499 0.027712 rs11590447 0.21793 0.1076 0.013921 0.028690 rs11120986 0.22251 0.1076 0.014881 0.028690 [1] "" [1] "-----twbc1 and twbc2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs11994137 0.00754 0.11935 0.00843 0.152 rs16882778 0.00850 0.11935 0.00936 0.152 rs9429328 0.00438 0.38810 0.01376 0.152 rs10434944 0.55585 0.00431 0.01739 0.152 rs2463460 0.00408 0.86010 0.02005 0.152 rs7197684 0.01017 0.65459 0.03840 0.185 rs10798738 0.32148 0.01924 0.04178 0.185 rs10271646 0.61609 0.01313 0.04429 0.185 rs10241173 0.89424 0.01313 0.04852 0.185 rs11120989 0.07202 0.09951 0.04873 0.185 [1] "" [1] "-----thgb1 and thgb2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs10798738 0.1388 0.000738 0.00248 0.0825 rs10271646 0.0126 0.203671 0.02006 0.2647 rs10434944 0.8602 0.008577 0.02668 0.2647 rs6973410 0.0134 0.483187 0.03936 0.2647 rs6691275 0.2119 0.030941 0.04708 0.2647 rs7134205 0.5373 0.015264 0.05073 0.2647 rs2250910 0.6368 0.016720 0.06148 0.2647 rs11669028 0.1245 0.119283 0.08096 0.2647 rs11587056 0.1852 0.071636 0.08099 0.2647 rs10241173 0.0661 0.203671 0.08192 0.2647 [1] "" [1] "-----tplt1 and tplt2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs6870861 0.70463 0.00518 0.00310 0.0521 rs2551038 0.51618 0.00518 0.00426 0.0521 rs7134205 0.06731 0.01811 0.00447 0.0521 rs6872533 0.85295 0.00518 0.00492 0.0521 rs10434944 0.69144 0.00851 0.00555 0.0521 rs1649942 0.44185 0.01161 0.00998 0.0780 rs16882778 0.00830 0.30535 0.02724 0.1328 rs1398281 0.04628 0.10223 0.02731 0.1328 rs8021893 0.65432 0.02684 0.02817 0.1328 rs11994137 0.00878 0.30535 0.02830 0.1328 [1] "" [1] "-----tanc1 and tanc2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs2463460 0.00233 0.7516 0.0131 0.521 rs9429328 0.02347 0.1225 0.0225 0.521 rs10434944 0.23557 0.0377 0.0452 0.521 rs2551038 0.11110 0.0894 0.0597 0.521 rs7134205 0.29028 0.0429 0.0613 0.521 rs2250910 0.21582 0.0622 0.0715 0.521 rs10798738 0.36856 0.0597 0.0890 0.521 rs10241173 0.27446 0.0664 0.0893 0.521 rs7098297 0.17447 0.1106 0.0984 0.521 rs17467610 0.03265 0.6213 0.1026 0.521
proportional odds logistic regression
proportional odds logistic regression [1] "" [1] "-----resp1 and resp2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs1934632 0.8157 0.0337 0.000230 0.00432 rs10863854 0.8124 0.0337 0.000249 0.00432 rs4946514 0.1510 0.0260 0.003979 0.03996 rs8021893 0.0985 0.1251 0.006291 0.03996 rs3813437 0.1198 0.1101 0.007889 0.03996 rs11590447 0.1478 0.1101 0.008099 0.03996 rs11120986 0.1472 0.1101 0.008384 0.03996 rs11587056 0.1498 0.1101 0.009217 0.03996 rs7519687 0.4270 0.1101 0.012985 0.05004 rs1934626 0.8071 0.1251 0.014763 0.05030 [1] "" [1] "-----twbc1 and twbc2" p12019 p13881 pall qall rs9429328 0.0150 0.0369 0.0117 0.0407 p12019 p13881 pall qall rs11994137 0.00589 0.10460 0.00245 0.0147 rs16882778 0.00622 0.10460 0.00265 0.0147 rs9429328 0.01503 0.03689 0.01174 0.0407 rs10798738 0.32392 0.01827 0.01472 0.0407 rs7197684 0.00549 0.49408 0.01940 0.0429 rs11120989 0.03439 0.18783 0.02826 0.0502 rs7134205 0.33881 0.00284 0.03494 0.0502 rs2463460 0.00535 0.99993 0.03689 0.0502 rs3751803 0.07388 0.10460 0.04087 0.0502 rs7867305 0.23608 0.05348 0.04579 0.0506 [1] "" [1] "-----thgb1 and thgb2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs10798738 0.26391 0.000279 5.01e-05 0.00195 rs633280 0.75176 0.023558 1.42e-02 0.18566 rs2250910 0.70816 0.026157 1.67e-02 0.18566 rs10271646 0.01116 0.170033 2.05e-02 0.18566 rs1398281 0.00747 0.559424 2.74e-02 0.18566 rs7867305 0.05931 0.057104 2.86e-02 0.18566 rs6973410 0.00940 0.559424 4.09e-02 0.22773 rs11669028 0.52399 0.020350 5.77e-02 0.28073 rs6691275 0.40797 0.058301 7.02e-02 0.30353 rs7519687 0.61042 0.058301 8.61e-02 0.32282 [1] "" [1] "-----tplt1 and tplt2" p12019 p13881 pall qall rs7134205 0.0310 0.0384 0.00888 0.115 p12019 p13881 pall qall rs7134205 0.03104 0.03837 0.00888 0.115 rs2551038 0.46418 0.00563 0.01196 0.115 rs6870861 0.78226 0.00563 0.01232 0.115 rs6872533 0.65253 0.00563 0.01421 0.115 rs1426897 0.11996 0.04167 0.02674 0.148 rs1398281 0.14305 0.03696 0.02815 0.148 rs7519687 0.00742 0.79520 0.03208 0.148 rs16882778 0.03559 0.15603 0.04312 0.166 rs11994137 0.03871 0.15603 0.04616 0.166 rs10863854 0.38567 0.03913 0.05968 0.193 [1] "" [1] "-----tanc1 and tanc2" [1] p12019 p13881 pall qall <0 rows> (or 0-length row.names) p12019 p13881 pall qall rs11669028 0.03672 0.0828 0.0252 0.235 rs2463460 0.00737 0.9530 0.0337 0.235 rs10798738 0.38750 0.0183 0.0357 0.235 rs7134205 0.37458 0.0347 0.0502 0.235 rs2551038 0.07542 0.1092 0.0510 0.235 rs7867305 0.72487 0.0253 0.0650 0.235 rs11120989 0.07055 0.1385 0.0656 0.235 rs9429328 0.57531 0.0355 0.0699 0.235 rs2250910 0.25919 0.0549 0.0754 0.235 rs11994137 0.13532 0.1046 0.0911 0.235
Subscribe to:
Posts (Atom)