|
1 |
| -filter(group=="block") %>% |
2 |
| -mutate(grp = str_remove(term, "var_")) %>% |
3 |
| -select(group, grp, estimate) %>% |
4 |
| -rename(variance = estimate, |
5 |
| -effect = group) |
6 |
| -, |
7 |
| -diag_lme$modelStruct$varStruct %>% |
8 |
| -coef(unconstrained = FALSE, allCoef = TRUE) %>% |
9 |
| -enframe(name = "grp", value = "varStruct") %>% |
10 |
| -mutate(sigma = diag_lme$sigma) %>% |
11 |
| -mutate(StandardError = sigma * varStruct) %>% |
12 |
| -mutate(variance = StandardError ^ 2) %>% |
13 |
| -mutate(effect = "Residual") %>% |
14 |
| -select(effect, grp, variance) |
15 |
| -) |
16 |
| -bind_rows( |
17 |
| -diag_glmmTMB %>% |
18 |
| -tidy(effects = "ran_pars", scales = "vcov") %>% |
19 |
| -filter(group == "block") %>% |
20 |
| -mutate(grp = str_remove(term, "var__")) %>% |
21 |
| -select(group, grp, estimate) %>% |
22 |
| -rename(variance = estimate, |
23 |
| -effect = group) |
24 |
| -, |
25 |
| -diag_glmmTMB %>% |
26 |
| -mixedup::extract_cor_structure(which_cor="diag") %>% |
27 |
| -pivot_longer( cols = 2:3, values_to ="variance", names_to="grp") %>% |
28 |
| -rename(effect = group) |
29 |
| -) |
30 |
| -bind_rows( |
31 |
| -diag_lme %>% |
32 |
| -tidy(effects = "ran_pars", scales = "vcov") %>% |
33 |
| -filter(group=="block") %>% |
34 |
| -mutate(grp = str_remove(term, "var_")) %>% |
35 |
| -select(group, grp, estimate) %>% |
36 |
| -rename(variance = estimate, |
37 |
| -effect = group) |
38 |
| -, |
39 |
| -diag_lme$modelStruct$varStruct %>% |
40 |
| -coef(unconstrained = FALSE, allCoef = TRUE) %>% |
41 |
| -enframe(name = "grp", value = "varStruct") %>% |
42 |
| -mutate(sigma = diag_lme$sigma) %>% |
43 |
| -mutate(StandardError = sigma * varStruct) %>% |
44 |
| -mutate(variance = StandardError ^ 2) %>% |
45 |
| -mutate(effect = "Residual") %>% |
46 |
| -select(effect, grp, variance) |
47 |
| -) |
48 |
| -bind_rows( |
49 |
| -diag_glmmTMB %>% |
50 |
| -tidy(effects = "ran_pars", scales = "vcov") %>% |
51 |
| -filter(group == "block") %>% |
52 |
| -mutate(grp = str_remove(term, "var__")) %>% |
53 |
| -select(group, grp, estimate) %>% |
54 |
| -rename(variance = estimate, |
55 |
| -effect = group) |
56 |
| -, |
57 |
| -diag_glmmTMB %>% |
58 |
| -mixedup::extract_cor_structure(which_cor="diag") %>% |
59 |
| -pivot_longer( cols = 2:3, values_to ="variance", names_to="grp") %>% |
60 |
| -rename(effect = group) |
61 |
| -) |
62 |
| -diag_glmmTMB %>% |
63 |
| -mixedup::extract_cor_structure(which_cor="diag") |
64 |
| -bind_rows( |
65 |
| -diag_glmmTMB %>% |
66 |
| -tidy(effects = "ran_pars", scales = "vcov") %>% |
67 |
| -filter(group == "block") %>% |
68 |
| -mutate(grp = str_remove(term, "var__")) %>% |
69 |
| -select(group, grp, estimate) %>% |
70 |
| -rename(variance = estimate, |
71 |
| -effect = group) |
72 |
| -, |
73 |
| -diag_glmmTMB %>% |
74 |
| -mixedup::extract_cor_structure(which_cor="diag", scales = "vcov") %>% |
75 |
| -pivot_longer( cols = 2:3, values_to ="variance", names_to="grp") %>% |
76 |
| -rename(effect = group) |
77 |
| -) |
78 |
| -bind_rows( |
79 |
| -diag_glmmTMB %>% |
80 |
| -tidy(effects = "ran_pars", scales = "vcov") %>% |
81 |
| -filter(group == "block") %>% |
82 |
| -mutate(grp = str_remove(term, "var__")) %>% |
83 |
| -select(group, grp, estimate) %>% |
84 |
| -rename(variance = estimate, |
85 |
| -effect = group) |
86 |
| -, |
87 |
| -diag_glmmTMB %>% |
88 |
| -mixedup::extract_cor_structure(which_cor="diag") %>% |
89 |
| -mutate(estimate = estimate ^2) %>% |
90 |
| -pivot_longer( cols = 2:3, values_to ="variance", names_to="grp") %>% |
91 | 1 | rename(effect = group)
|
92 | 2 | )
|
93 | 3 | diag_glmmTMB %>%
|
@@ -510,3 +420,93 @@ knitr::purl(input = purl_files %>% slice(i) %>% pull(Rmd),
|
510 | 420 | output = purl_files %>% slice(i) %>% pull(R),
|
511 | 421 | documentation = 0)
|
512 | 422 | }
|
| 423 | +pacman::p_load(conflicted, |
| 424 | +tidyverse, |
| 425 | +nlme, glmmTMB, |
| 426 | +broom.mixed, |
| 427 | +emo, flair) |
| 428 | +# package function conflicts |
| 429 | +conflict_prefer("filter", "dplyr") |
| 430 | +pacman::p_load(kableExtra) |
| 431 | +options(knitr.kable.NA = '') |
| 432 | +# this is never seen but needed for code below |
| 433 | +dat <- iris %>% |
| 434 | +mutate(unit = 1:n() %>% as.factor()) %>% |
| 435 | +rename(y=Sepal.Length, |
| 436 | +fixedeffects=Species, |
| 437 | +randomeffects=Sepal.Width) %>% |
| 438 | +as_tibble() |
| 439 | +# Chunk 1 |
| 440 | +pacman::p_load(conflicted, |
| 441 | +tidyverse, |
| 442 | +nlme, glmmTMB, |
| 443 | +broom.mixed, |
| 444 | +emo, flair) |
| 445 | +# package function conflicts |
| 446 | +conflict_prefer("filter", "dplyr") |
| 447 | +# Chunk 2 |
| 448 | +pacman::p_load(kableExtra) |
| 449 | +options(knitr.kable.NA = '') |
| 450 | +# Chunk 3: hideme |
| 451 | +# this is never seen but needed for code below |
| 452 | +dat <- iris %>% |
| 453 | +mutate(unit = 1:n() %>% as.factor()) %>% |
| 454 | +rename(y=Sepal.Length, |
| 455 | +fixedeffects=Species, |
| 456 | +randomeffects=Sepal.Width) %>% |
| 457 | +as_tibble() |
| 458 | +# Chunk 4: StandErrMod |
| 459 | +# |
| 460 | +# |
| 461 | +StandErrMod <- glmmTMB( |
| 462 | +y ~ |
| 463 | +fixedeffects + |
| 464 | +(1 | randomeffects), |
| 465 | +REML = TRUE, |
| 466 | +data = dat |
| 467 | +) |
| 468 | +# Chunk 5 |
| 469 | +decorate("StandErrMod") |
| 470 | +# Chunk 6: PseudErrMod |
| 471 | +dat <- dat %>% |
| 472 | +mutate(unit = as.factor(1:n())) |
| 473 | +PseudErrMod <- glmmTMB( |
| 474 | +y ~ |
| 475 | +fixedeffects + |
| 476 | +(1 | randomeffects) + |
| 477 | +(1 | unit), # Pseudo Err |
| 478 | +dispformula = ~ 0, # ErrVar = 0 |
| 479 | +REML = TRUE, |
| 480 | +data = dat |
| 481 | +) |
| 482 | +decorate("PseudErrMod") %>% flair_lines(c(1,2,8,9)) |
| 483 | +purl_files <- list.files(pattern = ".Rmd") %>% |
| 484 | +tibble(Rmd = .) %>% |
| 485 | +filter(Rmd %not_in% c("_hilang_setup.Rmd", |
| 486 | +"0contactinfo.Rmd", |
| 487 | +"index.Rmd", |
| 488 | +"model_selection.Rmd", |
| 489 | +"sources.Rmd", |
| 490 | +"variance_structures.Rmd", |
| 491 | +"weighted_two_stage.Rmd")) %>% |
| 492 | +mutate(R = paste0("Rpurl/",str_sub(Rmd, 1, -3))) |
| 493 | +setwd("D:/Coding/MMFAIR") |
| 494 | +library("rmarkdown") |
| 495 | +`%not_in%` = Negate(`%in%`) |
| 496 | +### create purled R files ### |
| 497 | +purl_files <- list.files(pattern = ".Rmd") %>% |
| 498 | +tibble(Rmd = .) %>% |
| 499 | +filter(Rmd %not_in% c("_hilang_setup.Rmd", |
| 500 | +"0contactinfo.Rmd", |
| 501 | +"index.Rmd", |
| 502 | +"model_selection.Rmd", |
| 503 | +"sources.Rmd", |
| 504 | +"variance_structures.Rmd", |
| 505 | +"weighted_two_stage.Rmd")) %>% |
| 506 | +mutate(R = paste0("Rpurl/",str_sub(Rmd, 1, -3))) |
| 507 | +for(i in 1:nrow(purl_files)){ |
| 508 | +knitr::purl(input = purl_files %>% slice(i) %>% pull(Rmd), |
| 509 | +output = purl_files %>% slice(i) %>% pull(R), |
| 510 | +documentation = 0) |
| 511 | +} |
| 512 | +purl_files |
0 commit comments