I'm trying to loop over multiple multilevel ordered logistic regressions with random intercepts on the country, dropping one observation at a time from the main dataset, while producing one massive, augmented tidy data frame with the results. Given that I am just having trouble with the looping, and I can produce everything one data frame at a time, I will first show what I am trying to do with two data frames. Then, I will (poorly) attempt to solve the problem with a loop, which is the part for which I would greatly appreciate your assistance.
First, let's load some libraries:
# load libraries
library(dplyr)
library(purrr)
library(tibble)
library(ordinal)
library(reprex)
library(magrittr)
Next, let's create some sample data:
# create original df in the right format
data0 = data.frame(country = rep(c("Algeria","Belgium","Canada","Denmark","England", "France"),times=10),
x1 = rep(0:1),times=30,
x2 = rnorm(n = 60, mean = 100, sd = 5),
y = rep(1:6),times=10)
data0$country = factor(data0$country)
data0$y = factor(data0$y)
data0 <- data0 %>% dplyr::select(country,x1,x2,y)
Building on this post, I'll create a custom tidier for clmm2 models:
tidy.clmm2 <-
function(fit){
results = as.data.frame(coefficients(summary(fit)))
colnames(results) = c("estimate","std.error","statistic","p.value")
results %>% tibble::rownames_to_column("term")
}
Estimate the first model, tidy it, and augment it with confidence intervals, odds ratios, and model observation number:
# model 1: get the dataset with one less observation, removing the top row
data1 <- data0 %>% dplyr::slice(-1)
# model 1: estimate
m1 <- clmm2(y ~ x1 + x2, random=country, data=data1, Hess = TRUE)
summary(m1)
# model 1: get odds ratio and renames columns for model 1
or1 <- as.data.frame(exp(coef(m1)))
or1 <- or1 %>% rename(estimate.odds =1)
or1$term <- row.names(or1)
# model 1: tidy the model, get the CIs, and store the observation/model number
tidy_m1 <- tidy.clmm2(m1)
tidy_m1$obs <- m1$nobs[[1]]
tidy_m1$conf.low <- tidy_m1$estimate - (1.96 * tidy_m1$std.error)
tidy_m1$conf.high <- tidy_m1$estimate + (1.96 * tidy_m1$std.error)
# model 1: merge over the odds ratios and make final data
tidy_m1 <- left_join(tidy_m1, or1, by=c("term"))
Do the same for model 2, and bind the final output with model 1:
# model 2: get the dataset with one less observation, removing the top row
data2 <- data1 %>% dplyr::slice(-1)
# model 2: estimate
m2 <- clmm2(y ~ x1 + x2, random=country, data=data2, Hess = TRUE)
summary(m2)
# model 2: get odds ratio and renames columns for model 1
or2 <- as.data.frame(exp(coef(m2)))
or2 <- or2 %>% rename(estimate.odds =1)
or2$term <- row.names(or2)
# model 2: tidy the model, get the CIs, and store the observation/model number
tidy_m2 <- tidy.clmm2(m2)
tidy_m2$obs <- m2$nobs[[1]]
tidy_m2$conf.low <- tidy_m2$estimate - (1.96 * tidy_m2$std.error)
tidy_m2$conf.high <- tidy_m2$estimate + (1.96 * tidy_m2$std.error)
# model 2: merge over the odds ratios and make final data
tidy_m2 <- left_join(tidy_m2, or2, by=c("term"))
## Bind everything together to get a final data frame
tidy_final <- dplyr::bind_rows(tidy_m1,tidy_m2)
Here is my (poor) attempt at a loop, building on this excellent post:
# create vector to store observation/data frame numbers
vector <- 1:59
# start of model
df_final <- purrr::map_dfr(1:59,
function(i) data.frame(model = i,
tidy.clmm2(glm(as.formula(paste0('y ~ ', i)),
random=country,
Hess = TRUE,
data = data[vector]))))
How do I fix the above in line with what I was doing one model at a time? Don't worry about the NaNs, as I have only provided a highly-stylized example to facilitate responses.
First I would improve your tidier function to do all of the extra tidying you were doing for each model (I'm assuming a tidier for this model class doesn't already exist in some package or other):
tidy.clmm2 <-
function(fit){
results = as.data.frame(coefficients(summary(fit)))
colnames(results) = c("estimate","std.error","statistic","p.value")
tidy_fit <- results %>% tibble::rownames_to_column("term")
or1 <- as.data.frame(exp(coef(fit)))
or1 <- or1 %>% rename(estimate.odds =1)
or1$term <- row.names(or1)
#tidy the model, get the CIs, and store the observation/model number
tidy_fit$obs <- fit$nobs[[1]]
tidy_fit$conf.low <- tidy_fit$estimate - (1.96 * tidy_fit$std.error)
tidy_fit$conf.high <- tidy_fit$estimate + (1.96 * tidy_fit$std.error)
#merge over the odds ratios and make final data
left_join(tidy_fit, or1, by=c("term"))
}
I don't know purrr, but you can then use the base function lapply (which I guess is similar) to estimate the model and the get the tidy parameters on successively smaller datasets, befor passing the results to bind_rows. Here I only do it for the first 10:
df_final <- lapply(1:10, function(i) {
mod <- clmm2(y ~ x1 + x2, random=country, data=data0[(i+1):60,], Hess = TRUE)
tidy.clmm2(mod)
}) %>% bind_rows(, .id="Dataset")
> df_final
Dataset term estimate std.error statistic p.value obs conf.low conf.high estimate.odds
1 1 1|2 -6.438544e+01 NaN NaN NaN 59 NaN NaN 1.090838e-28
2 1 2|3 -3.518858e+01 NaN NaN NaN 59 NaN NaN 5.221475e-16
3 1 3|4 2.961620e-01 NaN NaN NaN 59 NaN NaN 1.344688e+00
4 1 4|5 2.971124e+01 2.207366e+01 1.346004e+00 1.783011e-01 59 -13.55313129 72.97561932 8.006253e+12
5 1 5|6 7.273474e+01 NaN NaN NaN 59 NaN NaN 3.875214e+31
6 1 x1 2.054609e+01 NaN NaN NaN 59 NaN NaN 8.376306e+08
7 1 x2 -8.899366e-02 3.978817e-03 -2.236686e+01 8.275733e-111 59 -0.09679214 -0.08119518 9.148514e-01
8 2 1|2 -7.093800e+01 NaN NaN NaN 58 NaN NaN 1.556034e-31
9 2 2|3 -2.937359e+01 7.293763e+00 -4.027221e+00 5.644000e-05 58 -43.66936969 -15.07781910 1.750693e-13
10 2 3|4 -4.496101e+00 5.083061e+00 -8.845263e-01 3.764122e-01 58 -14.45889969 5.46669816 1.115240e-02
11 2 4|5 2.637413e+01 1.028467e+00 2.564412e+01 4.917528e-145 58 24.35833179 28.38992198 2.845364e+11
12 2 5|6 8.163547e+01 5.497370e-02 1.484991e+03 0.000000e+00 58 81.52771782 81.74321472 2.843364e+35
13 2 x1 1.821630e+01 5.424801e+00 3.357967e+00 7.851802e-04 58 7.58369175 28.84891081 8.151530e+07
14 2 x2 -7.476304e-02 NaN NaN NaN 58 NaN NaN 9.279633e-01
15 3 1|2 -1.105259e+02 NaN NaN NaN 57 NaN NaN 9.981611e-49
16 3 2|3 -4.660975e+01 NaN NaN NaN 57 NaN NaN 5.723279e-21
17 3 3|4 -1.713909e+00 5.880686e+00 -2.914471e-01 7.707094e-01 57 -13.24005308 9.81223521 1.801602e-01
18 3 4|5 4.382707e+01 NaN NaN NaN 57 NaN NaN 1.081073e+19
19 3 5|6 1.265335e+02 3.425069e-03 3.694335e+04 0.000000e+00 57 126.52681215 126.54023843 8.970400e+54
20 3 x1 3.988252e+01 5.999608e-04 6.647520e+04 0.000000e+00 57 39.88133918 39.88369103 2.092937e+17
21 3 x2 -2.175413e-01 2.615869e-04 -8.316215e+02 0.000000e+00 57 -0.21805405 -0.21702863 8.044943e-01
22 4 1|2 -6.646351e+01 NaN NaN NaN 56 NaN NaN 1.365411e-29
23 4 2|3 -3.064726e+01 1.722949e+01 -1.778768e+00 7.527783e-02 56 -64.41706724 3.12253810 4.898489e-14
24 4 3|4 -3.327152e+00 NaN NaN NaN 56 NaN NaN 3.589518e-02
25 4 4|5 3.188512e+01 NaN NaN NaN 56 NaN NaN 7.039323e+13
26 4 5|6 7.214246e+01 NaN NaN NaN 56 NaN NaN 2.143249e+31
27 4 x1 2.524262e+01 1.420819e+01 1.776624e+00 7.563011e-02 56 -2.60544100 53.09068131 9.177632e+10
28 4 x2 -1.139253e-01 5.874164e-04 -1.939430e+02 0.000000e+00 56 -0.11507663 -0.11277396 8.923246e-01
29 5 1|2 -5.520472e+01 1.024174e-03 -5.390171e+04 0.000000e+00 55 -55.20673057 -55.20271580 1.058994e-24
30 5 2|3 -3.324923e+01 NaN NaN NaN 55 NaN NaN 3.631150e-15
31 5 3|4 7.955460e-01 5.817109e+00 1.367597e-01 8.912208e-01 55 -10.60598823 12.19708026 2.215650e+00
32 5 4|5 2.816002e+01 5.061437e+00 5.563641e+00 2.642027e-08 55 18.23960201 38.08043320 1.697228e+12
33 5 5|6 6.669283e+01 NaN NaN NaN 55 NaN NaN 9.211492e+28
34 5 x1 3.351064e+01 1.024452e-03 3.271080e+04 0.000000e+00 55 33.50863540 33.51265125 3.576741e+14
35 5 x2 -1.850497e-01 1.173298e-03 -1.577175e+02 0.000000e+00 55 -0.18734938 -0.18275005 8.310630e-01
36 6 1|2 -3.286871e+01 NaN NaN NaN 54 NaN NaN 5.312531e-15
37 6 2|3 -1.394418e+01 NaN NaN NaN 54 NaN NaN 8.792619e-07
38 6 3|4 1.402239e+01 NaN NaN NaN 54 NaN NaN 1.229831e+06
39 6 4|5 4.198923e+01 NaN NaN NaN 54 NaN NaN 1.720640e+18
40 6 5|6 6.091421e+01 NaN NaN NaN 54 NaN NaN 2.849095e+26
41 6 x1 2.791695e+01 9.151861e+00 3.050412e+00 2.285273e-03 54 9.97930271 45.85459762 1.330998e+12
42 6 x2 6.368881e-04 NaN NaN NaN 54 NaN NaN 1.000637e+00
43 7 1|2 -8.540551e+01 NaN NaN NaN 53 NaN NaN 8.106962e-38
44 7 2|3 -4.981897e+01 NaN NaN NaN 53 NaN NaN 2.311503e-22
45 7 3|4 1.323464e-02 NaN NaN NaN 53 NaN NaN 1.013323e+00
46 7 4|5 4.199212e+01 3.930682e+00 1.068317e+01 1.220384e-26 53 34.28798714 49.69626009 1.725630e+18
47 7 5|6 9.757102e+01 NaN NaN NaN 53 NaN NaN 2.368942e+42
48 7 x1 2.742401e+01 NaN NaN NaN 53 NaN NaN 8.130075e+11
49 7 x2 -1.149872e-01 NaN NaN NaN 53 NaN NaN 8.913776e-01
50 8 1|2 -3.428518e+01 2.418647e+01 -1.417536e+00 1.563264e-01 52 -81.69065068 13.12029594 1.288655e-15
51 8 2|3 -1.468742e+01 2.158542e+01 -6.804326e-01 4.962306e-01 52 -56.99484187 27.61999711 4.181514e-07
52 8 3|4 -2.123523e+00 2.585943e+01 -8.211793e-02 9.345529e-01 52 -52.80800668 48.56096089 1.196095e-01
53 8 4|5 1.428943e+01 2.595126e+01 5.506258e-01 5.818902e-01 52 -36.57503632 65.15390302 1.606283e+06
54 8 5|6 3.811502e+01 1.071360e+00 3.557631e+01 3.257489e-277 52 36.01515820 40.21488826 3.573915e+16
55 8 x1 8.597913e+00 2.399747e+01 3.582841e-01 7.201307e-01 52 -38.43712837 55.63295352 5.420333e+03
56 8 x2 -3.759098e-02 NaN NaN NaN 52 NaN NaN 9.631068e-01
57 9 1|2 -1.381079e+02 NaN NaN NaN 51 NaN NaN 1.048377e-60
58 9 2|3 -5.910938e+01 NaN NaN NaN 51 NaN NaN 2.133649e-26
59 9 3|4 -9.084277e+00 NaN NaN NaN 51 NaN NaN 1.134354e-04
60 9 4|5 5.875811e+01 NaN NaN NaN 51 NaN NaN 3.298542e+25
61 9 5|6 1.564752e+02 NaN NaN NaN 51 NaN NaN 9.043358e+67
62 9 x1 4.528686e+01 NaN NaN NaN 51 NaN NaN 4.654054e+19
63 9 x2 -2.132578e-01 2.314244e-05 -9.215010e+03 0.000000e+00 51 -0.21330318 -0.21321246 8.079478e-01
64 10 1|2 -5.537366e+01 2.398432e+01 -2.308744e+00 2.095778e-02 50 -102.38292260 -8.36439034 8.943892e-25
65 10 2|3 -2.567911e+01 1.644416e+01 -1.561594e+00 1.183836e-01 50 -57.90967366 6.55145288 7.042130e-12
66 10 3|4 -3.561160e+00 1.421535e+01 -2.505152e-01 8.021890e-01 50 -31.42323652 24.30091746 2.840587e-02
67 10 4|5 2.622061e+01 1.739894e+01 1.507024e+00 1.318045e-01 50 -7.88130243 60.32253224 2.440441e+11
68 10 5|6 6.147091e+01 NaN NaN NaN 50 NaN NaN 4.971402e+26
69 10 x1 1.720323e+01 NaN NaN NaN 50 NaN NaN 2.959845e+07
70 10 x2 -6.799849e-02 NaN NaN NaN 50 NaN NaN 9.342619e-01