-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathmixed-shrinkage.Rmd
40 lines (32 loc) · 1.03 KB
/
mixed-shrinkage.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
---
title: 'Shrinkage'
---
```{r, include=FALSE}
library(tidyverse)
library(pander)
library(lmerTest)
```
```{r}
set.seed(1)
within.sd <- 1
between.sd <- 1
df <- expand.grid(person=1:10, trial=1:10, group=factor(1:3)) %>%
mutate(person = factor(as.numeric(factor(paste(person, group))))) %>%
group_by(person) %>%
mutate(u = rnorm(1, 0, between.sd)) %>%
ungroup() %>%
mutate(y = rnorm(n(), 0, within.sd) + u + 1*as.numeric(group))
m <- lmer(y ~ group + (1|person), data=df)
df <- left_join(df,
ranef(m)$person %>% rownames_to_column(var="person") %>% rename(uhat=`(Intercept)`))
peeps.df <- df %>% group_by(person) %>% summarise_all(funs(first))
xb.person <- predict(m, re.form=NA, newdata=peeps.df)
ranef(m)$person %>%
rownames_to_column(var="person") %>%
left_join(peeps.df, .) %>%
mutate(person=person,
xb= xb.person, shrunk = `(Intercept)`+xb) %>%
ggplot(aes(x=person, y=shrunk, xend=person, yend=xb, color=group)) +
geom_segment(arrow=arrow(length = unit(0.03, "npc"))) +
facet_wrap(~group)
```