← Back to branch summary

~ubuntu-branches/ubuntu/utopic/r-cran-genabel/utopic

« back to all changes in this revision

Viewing changes to R/lilog.R

  • Committer: Package Import Robot
  • Author(s): Andreas Tille, Charles Plessy, Andreas Tille
  • Date: 2014-08-07 17:30:04 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20140807173004-24c20tmrw61yl5dz
Tags: 1.8-0-1
[ Charles Plessy ]
* debian/control: removed myself from Uploaders.

[ Andreas Tille ]
* New upstream version
* Moved debian/upstream to debian/upstream/metadata
* cme fix dpkg-control
* More detailed copyright

Show diffs side-by-side

added added

removed removed

Lines of Context:
R/lilog.R
 
1
#' LiLog (Linear to Logistic Beta estimator)
 
2
#'
 
3
#' Estimating Logistic Betas from a linear regression.
 
4
#'
 
5
#' The function transforms the betas from the linear regression of the
 
6
#' binomial trait into its logistic counterpart.
 
7
#' It is based on the following relationship: if we denote with
 
8
#' \eqn{\alpha} the coefficients of the logistic regression and with
 
9
#' \eqn{\beta} those coming from the linear regression we have that:
 
10
#'
 
11
#' \deqn{z=\alpha_0 + \alpha_{SNP} + \alpha_{COV1} \ldots +
 
12
#'  \alpha_{COVi} = \ln \left( \frac{p}{1-p} \right),}{%
 
13
#'  z = alpha0 + alphaSNP + alphaCOV1 + ... + alphaCOVi = ln(p/(1-p)),}
 
14
#' with
 
15
#' \deqn{p = \beta_0 + \beta_{SNP}}{%
 
16
#'       p = beta0 + betaSNP}
 
17
#' Covariates are disregarded and the relationship
 
18
#' \deqn{\alpha_0 + \alpha_{SNP} = \mathrm{logit}(\beta_0 + \beta_{SNP})}{%
 
19
#' alpha0 + \alphaSNP = logit(beta0 + betaSNP)}
 
20
#' holds.
 
21
#'
 
22
#' \eqn{\beta_0} is not the coefficient coming from the linear
 
23
#' regression but it is estimated as:
 
24
#' \deqn{\mathrm{prev} - 2  q  (1-q)  \beta - q^2  \beta}{%
 
25
#'       prev - 2 * q * (1-q) * beta - q^2 * beta}
 
26
#'
 
27
#' and \eqn{\alpha_0} is estimated as
 
28
#' \eqn{\mathrm{logit}(\beta_0)}{logit(beta0)}.
 
29
#' So \eqn{\alpha_{SNP}}, which is the parameter of interest, will be:
 
30
#' \deqn{\alpha_{SNP} = \mathrm{logit}(\beta_0 + \beta_{SNP}) -
 
31
#'            \mathrm{logit}(\beta_0).}{%
 
32
#'       alphaSNP = logit(beta0 + betaSNP) - logit(beta0).}
 
33
#'
 
34
#' \eqn{SE(\alpha_{SNP})}{SE(\alphaSNP)} can be estimated considering
 
35
#' that the following relationship must be true:
 
36
#' \deqn{\frac{\alpha_{SNP}}{SE(\alpha_{SNP})} =
 
37
#'       \frac{\beta_{SNP}}{SE(\beta_{SNP})},}{%
 
38
#'       alphaSNP/SEalphaSNP = betaSNP/SEbetaSNP,} which is equivalent to
 
39
#' \deqn{SE(\alpha_{SNP}) = \alpha_{SNP}
 
40
#'         \frac{SE(\beta_{SNP})}{\beta_{SNP}}.}{%
 
41
#'       SE(alphaSNP) = alphaSNP * SE(betaSNP)/betaSNP.}
 
42
#'
 
43
#' @param prev: Prevalence of the disease as observed from the data
 
44
#' @param beta: \eqn{\beta_{SNP}}{beta_SNP} as estimated from a linear
 
45
#'  regression of the binomial trait
 
46
#' @param SEbeta: Standard Error of \eqn{\beta}{beta} coming from the
 
47
#'  linear regression
 
48
#' @param q: Frequency of the coded allele.
 
49
#'
 
50
#' @return Dataframe containing \eqn{\beta}{beta} and
 
51
#' \eqn{SE(\beta)}{SE(beta)}.
 
52
#'
 
53
#' @author Nicola Pirastu, Lennart Karssen, Yurii Aulchenko
 
54
#'
 
55
#' @references
 
56
#'
 
57
#' \emph{Genome-wide feasible and meta-analysis oriented methods for
 
58
#' analysis of binary traits using mixed models},
 
59
#' Nicola Pirastu, Lennart C. Karssen, Pio D'adamo, Paolo Gasparini,
 
60
#' Yurii Aulchenko (Submitted).
 
61
#'
 
62
#' @seealso
 
63
#' For transforming results from \code{\link{formetascore}} see
 
64
#' \code{\link{formetascore2lilog}}
 
65
#' while for tranforming the output from ProbABEL use \code{\link{palinear2LiLog}}
 
66
#'
 
67
#' @examples
 
68
#'
 
69
#' # simulate a snp with q = 0.3
 
70
#' snp <- rbinom(n=1000, size=2, prob=0.3)
 
71
#' # estimate the probability according to a logistic model with beta 0.4
 
72
#' probabil <- 1/(1 + exp( -(snp * 0.4 - 1) ))
 
73
#' # simulate the binary trait
 
74
#' bin.trait <- rbinom(n=1000, size=1, prob=probabil)
 
75
#' # run linear regression
 
76
#' res.lin <- glm(bin.trait~snp)
 
77
#' # run logistic regression
 
78
#' res.log <- glm(bin.trait~snp, family=binomial)
 
79
#'
 
80
#' # recover information from linear regression
 
81
#' snp.info <- summary(res.lin)$coefficients["snp",]
 
82
#'
 
83
#' # transform linear regression coefficients into logistic regression scale
 
84
#' LiLog(prev=mean(bin.trait), beta=snp.info[1], SEbeta=snp.info[2],
 
85
#'       q=mean(snp)/2)
 
86
#'
 
87
#' summary(res.log)$coefficients["snp", 1:2]
 
88
#'
 
89
#' # Example with qtscore
 
90
#' data(srdta)
 
91
#'
 
92
#' snp <- as.numeric(srdta[, 100])
 
93
#' probabil <- 1/(1 + exp( -(snp * 0.4 - 1) ))
 
94
#' bin.trait <- rbinom(n=nids(srdta), size=1, prob=probabil)
 
95
#' srdta@phdata$binary <- bin.trait
 
96
#' res <- qtscore(binary, data=srdta)
 
97
#' sum <- summary(srdta)
 
98
#'
 
99
#' lilog.res <- LiLog(prev=mean(bin.trait), beta=res@results$effB,
 
100
#'                    SEbeta=res@results$se_effB, q=sum$Q.2)
 
101
#' summary(glm(bin.trait~snp,
 
102
#'         family=binomial)$coefficients)["snp"1:2,] # logistic regression
 
103
#' lilog.res[100, ] # lilog result
 
104
#'
 
105
#' @keywords
 
106
#' LiLog
 
107
#'
 
108
 
 
109
LiLog <- function(prev, beta, SEbeta, q){
 
110
  Icept <- (prev - 2 * q * beta)
 
111
  alfa0 <- log(Icept / (1-Icept))
 
112
  T1    <- (Icept + beta)
 
113
  T1[which(T1<0 | T1>1)] <- NA
 
114
 
 
115
  log.beta   <- binomial()$linkfun(T1) - alfa0
 
116
  se.log     <- (log.beta * SEbeta) / beta
 
117
  res        <- data.frame(log.beta, se.log)
 
118
  names(res) <- c("beta", "sebeta")
 
119
  return(res)
 
120
}