A friend asked me whether I could create a loop that would run multiple regression models. She wanted to evaluate the association between 100 dependent variables (outcomes) and 100 independent variables (exposures), which means 10,000 regression models. Regression analyses with many dependent (outcome) and independent (exposure) variables are common in genetics.
So the models will look something like this (dx is the dependent variable, ix is the independent variable, and v are the other covariates):
dx1 = ix1 + v1 + v2 + v3 dx1 = ix2 + v1 + v2 + v3 dx1 = ix3 + v1 + v2 + v3 ... dx1 = ix100 + v1 + v2 + v3 dx2 = ix1 + v1 + v2 + v3 ... dx100 = ix100 + v1 + v2 + v3
The output should be a data frame with five columns: the dependent variable, the independent variable, the beta estimate, the standard error, and the p-value.
It will look something like this (the numbers are just for illustration purposes):
d i beta se pvalue 1 dx1 ix1 0.1 0.002 0.950 2 dx2 ix2 0.2 0.002 0.826 3 dx3 ix3 0.3 0.005 0.123
OK, now let’s begin. The dataset I received had all the variables in columns and the observations in rows (the data are not real, just random numbers for illustration purposes):
id dx1 dx2 ... dx100 ix1 ... 1x100 v1 v2 v3 10 324 124 ... 214 32 ... 32 ax b4 c3 11 431 982 ... 114 12 ... 77 ce b2 c5 12 545 123 ... 104 34 ... 11 ar c2 a5 ....
Position of the variables
First, create vectors that hold the positions of the dependent and independent variables in your dataset, along with empty vectors to store the results.
# outcome out_start=2 out_end= 101 out_nvar=out_end-out_start+1 out_variable=rep(NA, out_nvar) out_beta=rep(NA, out_nvar) out_se = rep(NA, out_nvar) out_pvalue=rep(NA, out_nvar) # exposure exp_start=102 exp_end=201 exp_nvar=exp_end-exp_start+1 exp_variable=rep(NA, exp_nvar) exp_beta=rep(NA, exp_nvar) exp_se = rep(NA, exp_nvar) exp_pvalue=rep(NA, exp_nvar) number=1
For Loop
I used a linear mixed-effects model, and therefore I loaded the lme4 library. The loop should also work with other regression analyses (e.g., linear regression) if you modify it according to your regression model. If you don’t know which part to modify, leave a comment below and I will try to help.
Like other loops, this one calls the variables of interest one by one and, for each of them, extracts and stores the beta, the standard error, and the p-value. Remember, this code is specific to linear mixed-effects models.
library(lme4)
for (i in out_start:out_end){
outcome = colnames(dat)[i]
for (j in exp_start:exp_end){
exposure = colnames(dat)[j]
model <- lmer(get(outcome) ~ get(exposure) + v1 + (1|v2) + (1|v3),
na.action = na.exclude,
data=dat)
Vcov <- vcov(model, useScale = FALSE)
beta <- fixef(model)
se <- sqrt(diag(Vcov))
zval <- beta / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
out_beta[number] = as.numeric(beta[2])
out_se[number] = as.numeric(se[2])
out_pvalue[number] = as.numeric(pval[2])
out_variable[number] = outcome
number = number + 1
exp_beta[number] = as.numeric(beta[2])
exp_se[number] = as.numeric(se[2])
exp_pvalue[number] = as.numeric(pval[2])
exp_variable[number] = exposure
number = number + 1
}
}
After the loop finishes, we create two data frames with the results:
outcome = data.frame(out_variable, out_beta, out_se, out_pvalue) exposure = data.frame(exp_variable, exp_beta, exp_se, exp_pvalue)
Management of the data frame
We now have two different data frames with our results, and we need to combine them into one. With the help of the tidyverse package, this is a simple task. Basically, we rename the variables so that both data frames use the same names, and then we merge them together.
library(tidyverse)
outcome = outcome %>%
rename(
variable = out_variable,
beta = out_beta,
se = out_se,
pvalue = out_pvalue
)
exposure = exposure %>%
rename(
variable = exp_variable,
beta = exp_beta,
se = exp_se,
pvalue = exp_pvalue
)
all = rbind(outcome, exposure)
all = na.omit(all)
head(all)
variable beta se pvalue
1 dx1 0.1 0.002 0.950
3 dx2 0.2 0.002 0.826
........
2 ix1 0.1 0.002 0.950
4 ix2 0.2 0.002 0.826
........
However, this is still not the data frame we are looking for. We need a data frame that has both the dependent and the independent variable in the same row. Therefore, we do the final transformation as follows:
data = all %>%
mutate(
type = substr(variable, 1, 2)
) %>%
spread(type, variable) %>%
rename(
d = dx,
i = ix
) %>%
mutate (
beta = round(beta, 5),
se = round(se, 5),
pvalue = round(pvalue, 5)
) %>%
select(d, i, beta, se, pvalue)
head(data)
d i beta se pvalue
1 dx1 ix1 0.1 0.002 0.950
2 dx2 ix2 0.2 0.002 0.826
3 dx3 ix3 0.3 0.005 0.123
The final data frame lists, for each model, the outcome, the exposure, and the corresponding estimates. I hope you find this post useful for your research and data analysis!
Thank you for this. I am trying to create a loop that can loop in 5 minute data observations and run regressions for each day between 2008 and 2011. Could you kindly spare a minute to spare me out?
Dear Dhana,
Hi! Thanks for your code and it did help me a lot! But I have questions when I tried to run the code. I am trying to run multiple uni variate analysis with one same outcome but thousands of variables for each model. I added ‘NA’=as.character(NA) to the position of variables and the code gives me error saying ” variable lengths differ (found for ‘get(exposure)’)”. I don’t know what’s going on. Could you help me figure it out the problem? I really appreciate it! Thanks!
The code is paste below, and my rows (13940 rows) are the variables whereas columns are the IDs (50 columns)
dat3.nna <- dat3.nna %>% mutate(id = row_number()) ##will change original ID to 1:n
# outcome
out_start=2
out_end= 13941
out_nvar=out_end-out_start+1
‘NA’=as.character(NA)
out_variable=rep(NA, out_nvar)
out_beta=rep(NA, out_nvar)
out_se = rep(NA, out_nvar)
out_pvalue=rep(NA, out_nvar)
# exposure
exp_start=13942
exp_end=13944
exp_nvar=exp_end-exp_start+1
exp_variable=rep(NA, exp_nvar)
exp_beta=rep(NA, exp_nvar)
exp_se = rep(NA, out_nvar)
exp_pvalue=rep(NA, exp_nvar)
number=1
library(lme4)
for (i in out_start:out_end){
outcome = colnames(dat3.nna)[i]
for (j in exp_start:exp_end){
exposure = colnames(dat3.nna)[j]
model <- lm(get(outcome) ~ get(exposure), na.action = na.exclude, data=dat3.nna) Vcov <- vcov(model, useScale = FALSE) beta <- fixef(model) se <- sqrt(diag(Vcov)) zval <- beta / se pval <- 2 * pnorm(abs(zval), lower.tail = FALSE) out_beta[number] = as.numeric(beta[2]) out_se[number] = as.numeric(se[2]) out_pvalue[number] = as.numeric(pval[2]) out_variable[number] = outcome number = number + 1 exp_beta[number] = as.numeric(beta[2]) exp_se[number] = as.numeric(se[2]) exp_pvalue[number] = as.numeric(pval[2]) exp_variable[number] = exposure number = number + 1 } }
Many thanks for the code, it turned out to be really useful! It is very handy and works perfectly for lm and glm models (and other models, with appropriate modifications).
I have a question regarding models with a more complex structure than the one in the code above. For example, I am trying to include an interaction term between the variable that is being replaced (“exposure” in your code above) and a co-factor (a fixed categorical variable with two levels). I have no problem at all running the model, which in fact works fine. However, I have not accomplished yet to extract the beta/se/pvalue for the interacting term. This surpass my level of skill on R, and I would like to know if you have any idea on how to do it.
I really hope you read and reply to my comment.
Thanks in advance.
Hi,
Thanks for the code. I am modifying it for glm and it works fine for me up until I hit the ‘tidyverse’ package part. Specifically, where does the “out_nobs” come from? I get a “Error: `out_nobs` contains unknown variables” when I hit this section.
Thanks for the code. How can one modify this script for lme?
Thanks for the fantastic code. I am trying to modify this code for the felm (panel regresion). could you please help me?
I would suggest first to use the same regression model that I have used in this example; then try to explore how the summary of flem works; finally replace the code in the script. Hope this helps. Good luck.
I, too, would love help modifying this code for felm – panel data with fixed effects. Any insight? Thanks!
Thanks for the code! It is very helpful but I got the error message saying “Error: bad ‘data’, but variables found in environment of formula: try specifying ‘formula’ as a formula rather than a string in the original model” when I was trying to run it. Do you know what is going on and how to fix that ?
Thank you.
Ann
Test your model without using the function. After you make sure works well, select few variables to test function.
Thanks for the code. It’s very helpful. I’m not great at coding, so am trying to understand why “number = number + 1” is repeated in the for loop.
+1 tells the script to look to the subsequent variable.
I understand, but you **repeat** number = number+1 twice, that’s why I was confused. Wouldn’t the below code work as well?
# load libraries
require(“lme4”)
require(“tidyverse”)
# simulate dataset
set.seed(1000)
n <- 100
df1 <- data.frame(
x1 = rbinom(n, 1, 0.25),
x2 = rbinom(n, 1, 0.25),
x3 = rbinom(n, 1, 0.25),
x4 = rbinom(n, 1, 0.25),
x5 = rbinom(n, 1, 0.25),
y1 = rnorm(n, 0, 1),
y2 = rnorm(n, 0, 1),
y3 = rnorm(n, 0, 1),
y4 = rnorm(n, 0, 1),
y5 = rnorm(n, 0, 1),
male = rbinom(n, 1, 0.5),
z1 = factor(sample(1:10, 100, replace = TRUE))
)
dat <- df1
# outcome index
out_start = 6
out_end= 10
out_nvar = out_end-out_start + 1
# exposure index
exp_start = 1
exp_end = 5
exp_nvar = exp_end – exp_start + 1
# make empty vectors to store results
out_variable = rep(NA, out_nvar*exp_nvar)
exp_variable = rep(NA, out_nvar*exp_nvar)
exp_beta = rep(NA, out_nvar*exp_nvar)
exp_se = rep(NA, out_nvar*exp_nvar)
exp_pvalue = rep(NA, out_nvar*exp_nvar)
# start index
number = 1
# for loop using lme4 for mixed effects model where y ~ x + male + (1|z1)
for (i in out_start:out_end){
outcome = colnames(dat)[i]
for (j in exp_start:exp_end){
exposure = colnames(dat)[j]
model <- lmer(get(outcome) ~ get(exposure) + male + (1|z1),
na.action = na.exclude,
data=dat)
Vcov <- vcov(model, useScale = FALSE)
beta <- fixef(model)
se <- sqrt(diag(Vcov))
zval <- beta / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
exp_beta[number] = as.numeric(beta[2])
exp_se[number] = as.numeric(se[2])
exp_pvalue[number] = as.numeric(pval[2])
out_variable[number] = outcome
exp_variable[number] = exposure
number = number + 1
}
}
# output results
all <- data.frame(y = out_variable, x = exp_variable, beta = exp_beta, se = exp_se, pvalue = exp_pvalue)
all
Because we have repeated exposure and outcomes. I did not test your code.
This is extremely helpful. Thank you for this.
I am new to R and have a similar dataset. I however have only one dependent variable (dx) and would like to know how will the code change if I need to do linear regression.
I would appreciate your help.
Sorry for late reply! You have to slightly modify the code above.
First remove this part of code:
for (i in out_start:out_end){
outcome = colnames(dat)[i]
Second replace get(outcome) with your dependent variable.
Third remove this code:
out_beta[number] = as.numeric(beta[2])
out_se[number] = as.numeric(se[2])
out_pvalue[number] = as.numeric(pval[2])
out_variable[number] = outcome
number = number + 1
Finally remove this closing } brackets from the script.