daDoctoR/R/rep_glm.R

92 lines
2.2 KiB
R
Raw Normal View History

2018-10-03 10:32:10 +02:00
#' A repeated logistic regression function
2018-10-02 21:07:43 +02:00
#'
#' For bivariate analyses.
#' @param y Effect meassure.
#' @param v1 Main variable in model
#' @keywords logistic regression
#' @export
#' @examples
#' rep_glm()
rep_glm<-function(y,v1,string,ci=FALSE,data,v2=NULL,v3=NULL){
## x is data.frame of predictors, y is vector of an aoutcome as a factor
## output is returned as coefficient, or if or=TRUE as OR with 95 % CI.
## The confint() function is rather slow, causing the whole function to hang when including many predictors and calculating the ORs with CI.
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
require(broom)
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
d<-data
x<-select(d,one_of(c(string)))
m1<-length(coef(glm(y~v1,family = binomial())))
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
if (!is.factor(y)){stop("Some kind of error message would be nice, but y should be a factor!")}
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
if (ci==TRUE){
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
df<-data.frame(matrix(ncol = 4))
names(df)<-c("pred","or_ci","pv","t")
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
for(i in 1:ncol(x)){
m<-glm(y~v1+x[,i],family = binomial())
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
l<-suppressMessages(round(exp(confint(m))[-c(1:m1),1],2))
u<-suppressMessages(round(exp(confint(m))[-c(1:m1),2],2))
or<-round(exp(coef(m))[-c(1:m1)],2)
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
or_ci<-paste0(or," (",l," to ",u,")")
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
pv<-round(tidy(m)$p.value[-c(1:m1)],3)
pv<-ifelse(pv<0.001,"<0.001",pv)
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
t <- ifelse(pv<=0.1|pv=="<0.001","include","drop")
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
pv <- ifelse(pv<=0.05|pv=="<0.001",paste0("*",pv),
ifelse(pv>0.05&pv<=0.1,paste0(".",pv),pv))
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
v<-x[,i]
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
if (is.factor(v)){
pred<-paste(names(x)[i],levels(v)[-1],sep = "_")
}
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
else {pred<-names(x)[i]}
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
df<-rbind(df,cbind(pred,or_ci,pv,t))
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
}}
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
if (ci==FALSE){
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
df<-data.frame(matrix(ncol = 4))
names(df)<-c("pred","b","pv","t")
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
for(i in 1:ncol(x)){
m<-glm(y~v1+x[,i],family = binomial())
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
b<-round(coef(m)[-c(1:m1)],3)
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
pv<-round(tidy(m)$p.value[-c(1:m1)],3)
pv<-ifelse(pv<0.001,"<0.001",pv)
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
t <- ifelse(pv<=0.1|pv=="<0.001","include","drop")
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
pv <- ifelse(pv<=0.05|pv=="<0.001",paste0("*",pv),
ifelse(pv>0.05&pv<=0.1,paste0(".",pv),pv))
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
v<-x[,i]
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
if (is.factor(v)){
pred<-paste(names(x)[i],levels(v)[-1],sep = "_")
}
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
else {pred<-names(x)[i]}
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
df<-rbind(df,cbind(pred,b,pv,t))
2018-10-03 10:32:10 +02:00
2018-10-02 21:07:43 +02:00
}}
2018-10-03 10:32:10 +02:00
result<-df
2018-10-02 21:07:43 +02:00
return(df)
}