hysteresis_modeling.Rmd

---
title: "Hysteresis modeling"
output: html_notebook
---



```{r}
library(ggplot2)
library(dplyr)
library(hysteresis)

```



```{r}
load("akk_df.rda")
  
  
  
```



```{r, fig.height=4, fig.width=8}

pertubation_press_plot =
akk_df %>%
  select(Time,DSSpress, press_cumul_day) %>%
  unique %>%
  
  ggplot(aes(y=Time+1, x=press_cumul_day %>% as.character %>% as.numeric)) +
  #geom_point() +
  geom_path(color=c(rep(c("red","pink"),2),"red",rep("blue",4)), arrow = arrow(length=unit(0.20,"cm"), type = "closed")) + 
  ggrepel::geom_text_repel(aes(y=Time+1, x=press_cumul_day %>% as.character %>% as.numeric, label=Time, col=DSSpress), alpha = 1) +
  geom_point(aes(y=Time+1, x=press_cumul_day %>% as.character %>% as.numeric, col=DSSpress), alpha = 0) +
  coord_flip() + xlab("press\nperturbation") +
  geom_vline(xintercept = 0, linetype =2) + ylab("Time\n(days)") +
  scale_color_manual("", values=c("blue","red"), labels=c("Washout period\n(Water)","Chronic challenge\n(DSS)")) +
  guides(colour = guide_legend("Treatment", override.aes = list(size = 3, alpha = 1)))

pertubation_press_plot + danr::theme_dan()  +
  theme(legend.position = c(0.7,0.75))

```

```{r}

akk_df


```




```{r}


akk_df %>%
  arrange(press_cumul_day) %>%
  filter(hysteresis_grp != "resistant\nprofile") %>%
  #filter(hysteresis_grp == "no resilient\nprofile") %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  ggplot() + 
  geom_point(aes(y=log10(genus+10^-6),x=press_cumul_day %>% as.character %>% as.numeric, col=DSSpress)) +
  geom_path(aes(y=log10(genus+10^-6),x=press_cumul_day %>% as.character %>% as.numeric, group=Rat_ID)) +
  facet_wrap(~hysteresis_grp)



```


```{r}

akk_df %>%
  filter(press_cumul_day != -6) %>%
  filter(hysteresis_grp == "resistant\nprofile") %>%
  #filter(hysteresis_grp != "no resilient\nprofile") %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  select(genus,press_cumul_day,Rat_ID,hysteresis_grp) %>%
  mutate(genus_norm = -log10(genus+10^-6) %>% scale(), press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() %>% scale() ) %>%
  #ggplot() + geom_point(aes(x=press_cumul_day,y=genus_norm))
  
  
  #filter(Rat_ID %in% c("R13","R17")) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  with(.,floop(.$press_cumul_day, .$genus_norm,  method="harmonic2",n=3,m=3)) -> models_resistant


akk_df %>%
  filter(press_cumul_day != -6) %>%
  filter(hysteresis_grp != "resistant\nprofile") %>%
  filter(hysteresis_grp != "no resilient\nprofile") %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  select(genus,press_cumul_day,Rat_ID,hysteresis_grp) %>%
  mutate(genus_norm = -log10(genus+10^-6) %>% scale(), press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() %>% scale() ) %>%
  #ggplot() + geom_point(aes(x=press_cumul_day,y=genus_norm))
  
  
  #filter(Rat_ID %in% c("R13","R17")) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  with(.,floop(.$press_cumul_day, .$genus_norm,  method="harmonic2", n=3,m=3)) -> models_resilient



akk_df %>%
  filter(press_cumul_day != -6) %>%
  filter(hysteresis_grp != "resistant\nprofile") %>%
  filter(hysteresis_grp == "no resilient\nprofile") %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  select(genus,press_cumul_day,Rat_ID,hysteresis_grp) %>%
  mutate(genus_norm = -log10(genus+10^-6) %>% scale(), press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() %>% scale() ) %>%
  #ggplot() + geom_point(aes(x=press_cumul_day,y=genus_norm))
  
  
  #filter(Rat_ID %in% c("R13","R17")) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  with(.,floop(.$press_cumul_day, .$genus_norm,  method="harmonic2", n=3,m=3)) -> models_no_resilient


plot(models_resistant)
plot(models_resilient)
plot(models_no_resilient)


data.frame(x=models_no_resilient$x, y=models_no_resilient$y) %>%
    ggplot() + geom_point(aes(x=y,y=-x))



mloop(b.x=, n.points = 24)


models$Estimates[,"R13"]


summary(models_resistant)
summary(models_resilient)
summary(models_no_resilient)


```



```{r}


akk_df %>%
  #filter(press_cumul_day != -6) %>%
  filter(Study_group %in% c("G5","G6","G7")) %>%
  #filter(hysteresis_grp == "resistant\nprofile") %>%
  #filter(hysteresis_grp != "no resilient\nprofile") %>%
  filter(!(Rat_ID == "R67" & Time == "63")) %>%
  filter(!(Rat_ID == "R66" & Time == "68")) %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  select(genus,press_cumul_day,Rat_ID,hysteresis_grp, Time) %>%
  mutate(genus_norm = log10(genus+10^-6) , press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() ) %>%
  #ggplot() + geom_point(aes(x=press_cumul_day,y=genus_norm))
  
  
  #filter(Rat_ID %in% c("R13","R17")) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  #with(.,floop(.$press_cumul_day, .$genus_norm,  method="harmonic2", times=Time+1, subjects = .$Rat_ID,n=1,m=1)) -> models_all
#with(.,floop(.$press_cumul_day, .$genus_norm,  method= "geometric", subjects = .$Rat_ID)) -> models_all
with(.,floop(.$press_cumul_day, .$genus_norm,  method="harmonic2", times=Time+1, subjects = .$Rat_ID, n=3,m=1)) -> models_all
  
models_all_summary = summary(models_all, seed=444)

#plot(models_all_summary)

akk_df %>%
  select(Rat_ID, hysteresis_grp) %>%
  unique()


models_all_summary$values %>%
  filter(Parameter=="retention") %>%
  merge(.,
        akk_df %>%
  select(Rat_ID, hysteresis_grp) %>%
  unique(),
  by.x="Subject", by.y="Rat_ID"
  ) %>%
  mutate(B.q0.025 = ifelse(Boot.Estimate<0, -B.q0.025, B.q0.025),
         B.q0.975 = ifelse(Boot.Estimate<0, -B.q0.975, B.q0.975)) %>%
  mutate(Boot.Estimate = ifelse(Boot.Estimate<0, -Boot.Estimate, Boot.Estimate)) %>%
  ggplot() + geom_errorbar(aes(x=forcats::fct_reorder(Subject,Boot.Estimate), ymin = B.q0.025, ymax = B.q0.975, col=hysteresis_grp), width = 0.6) + 
  geom_point(aes(x=forcats::fct_reorder(Subject,Boot.Estimate),y=Boot.Estimate, col=hysteresis_grp)) +
  #facet_wrap(~hysteresis_grp, scales = "free_x") +
  geom_hline(yintercept = 0, col="red") +
  ylab("Hysteresis retention parameter\n(n=1000 Bootstrap)") + 
  xlab("Rats") + cowplot::theme_cowplot() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size=8))
  
models_all_summary$values %>%
  filter(Parameter=="b.y") %>%
  merge(.,
        akk_df %>%
  select(Rat_ID, hysteresis_grp) %>%
  unique(),
  by.x="Subject", by.y="Rat_ID"
  ) %>%
  mutate(B.q0.025 = ifelse(Boot.Estimate<0, -B.q0.025, B.q0.025),
         B.q0.975 = ifelse(Boot.Estimate<0, -B.q0.975, B.q0.975)) %>%
  mutate(Boot.Estimate = ifelse(Boot.Estimate<0, -Boot.Estimate, Boot.Estimate)) %>%
  ggplot() + geom_errorbar(aes(x=forcats::fct_reorder(Subject,Boot.Estimate), ymin = B.q0.025, ymax = B.q0.975, col=hysteresis_grp), width = 0.6) + 
  geom_point(aes(x=forcats::fct_reorder(Subject,Boot.Estimate),y=Boot.Estimate, col=hysteresis_grp)) +
  #facet_wrap(~hysteresis_grp, scales = "free_x") +
  geom_hline(yintercept = 0, col="red") +
  ylab("Hysteresis saturation parameter\n(n=1000 Bootstrap)") + 
  xlab("Rats") + cowplot::theme_cowplot() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size=8))



models_all_summary$values %>%
  filter(Parameter=="cy") %>%
  merge(.,
        akk_df %>%
  select(Rat_ID, hysteresis_grp) %>%
  unique(),
  by.x="Subject", by.y="Rat_ID"
  ) %>%
  #mutate(B.q0.025 = ifelse(Boot.Estimate<0, -B.q0.025, B.q0.025),
  #       B.q0.975 = ifelse(Boot.Estimate<0, -B.q0.975, B.q0.975)) %>%
  #mutate(Boot.Estimate = ifelse(Boot.Estimate<0, -Boot.Estimate, Boot.Estimate)) %>%
  ggplot() + geom_errorbar(aes(x=forcats::fct_reorder(Subject,Boot.Estimate), ymin = B.q0.025, ymax = B.q0.975, col=hysteresis_grp), width = 0.6) + 
  geom_point(aes(x=forcats::fct_reorder(Subject,Boot.Estimate),y=Boot.Estimate, col=hysteresis_grp)) +
  #facet_wrap(~hysteresis_grp, scales = "free_x") +
  geom_hline(yintercept = 0, col="red") +
  ylab("Hysteresis center parameter\n(n=1000 Bootstrap)") + 
  xlab("Rats") + cowplot::theme_cowplot() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size=8))


models_all_summary$values %>%
  filter(Parameter=="retention")


models_all_summary$values %>%
  filter(Parameter=="cy")



#plot(models_all)
```

```{r}

akk_df %>%
  filter(press_cumul_day != -6) %>%
  filter(Study_group %in% c("G5","G6","G7")) %>%
  #filter(hysteresis_grp == "resistant\nprofile") %>%
  #filter(hysteresis_grp != "no resilient\nprofile") %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  select(genus,press_cumul_day,Rat_ID,hysteresis_grp) %>%
  mutate(genus_norm = -log10(genus+10^-6) %>% scale(), press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() %>% scale() ) %>%
  #ggplot() + geom_point(aes(x=press_cumul_day,y=genus_norm))
  
  
  filter(Rat_ID %in% c("R62")) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  with(.,floop(.$press_cumul_day, .$genus_norm,  method="harmonic2" ,n=3,m=3)) -> model_rat62


akk_df %>%
  #filter(press_cumul_day != -6) %>%
  filter(Study_group %in% c("G5","G6","G7")) %>%
  #filter(hysteresis_grp == "resistant\nprofile") %>%
  #filter(hysteresis_grp != "no resilient\nprofile") %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  #select(genus,press_cumul_day,Rat_ID,hysteresis_grp) %>%
  mutate(genus_norm = -log10(genus+10^-6) %>% scale(), press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() %>% scale() ) %>%
  #ggplot() + geom_point(aes(x=press_cumul_day,y=genus_norm))
  
  
  filter(Rat_ID %in% c("R62")) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  ggplot() + geom_text(aes(x=press_cumul_day, y=-genus_norm, label=Time, col=DSSpress))






```

## tipping point figure
```{r}

tip_akk = 0.011480422 #from bimodality test

akk_df %>%
  #filter(DSSpulse=="Water") %>%
  filter(Study_group %in% c("G5","G6","G7")) %>%
  ggplot() + 
  geom_line(aes(x=forcats::fct_reorder(Rat_ID,genus,mean) , y=genus+10^-6, group = forcats::fct_reorder(Rat_ID,genus,mean), col=hysteresis_grp)) +
  geom_point(aes(x=forcats::fct_reorder(Rat_ID,genus,mean) , y=genus+10^-6)) +
  geom_hline(yintercept = tip_akk, linetype=2) +
  scale_y_log10() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size=8))


```


## akkermansia hysteresis and tipping point



```{r}


tip_akk = 0.011480422 #from bimodality test


merge(
  
akk_df %>%
  filter(!(Rat_ID == "R67" & Time == "63")) %>%
  filter(!(Rat_ID == "R66" & Time == "68")) %>%
  #filter(DSSpulse=="Water") %>%
  filter(Study_group %in% c("G5","G6","G7")) %>%
  select(Rat_ID,genus,hysteresis_grp) %>%
  group_by(Rat_ID) %>%
  summarise(genus_min = min(genus), genus_max=max(genus)) %>%
  mutate(tipping_point = tip_akk) %>%
  mutate(tipping_point_crosses = ifelse(genus_min < tipping_point, "yes", "no")),




models_all_summary$values %>%
  filter(Parameter=="retention") %>%
  merge(.,
        akk_df %>%
          filter(!(Rat_ID == "R67" & Time == "63")) %>%
  filter(!(Rat_ID == "R66" & Time == "68")) %>%
  select(Rat_ID, hysteresis_grp) %>%
  unique(),
  by.x="Subject", by.y="Rat_ID"
  ) %>%
  mutate(B.q0.025 = ifelse(Boot.Estimate<0, -B.q0.025, B.q0.025),
         B.q0.975 = ifelse(Boot.Estimate<0, -B.q0.975, B.q0.975)) %>%
  mutate(Boot.Estimate = ifelse(Boot.Estimate<0, -Boot.Estimate, Boot.Estimate)) %>%
  select(Subject,B.q0.025,B.q0.975,Boot.Estimate,hysteresis_grp) %>%
  mutate(retention = ifelse((B.q0.025 > 0 &  B.q0.975 < 0) | (B.q0.025 < 0 &  B.q0.975 > 0)   , "no","yes")), by.x="Rat_ID",by.y="Subject") -> akk_tipping_retention

akk_tipping_retention %>%
  mutate(hysteresis_grp = case_when(
    
tipping_point_crosses == "no" ~ "resistant",
tipping_point_crosses == "yes" & retention == "yes" ~ "no resilient",
tipping_point_crosses == "yes" & retention == "no" ~ "resilient"
    
    
    
  )) -> akk_tipping_retention



```


tipping point and retention plots
```{r, fig.height=4, fig.width=4}

akk_tipping_retention %>%
  mutate(hysteresis_grp = hysteresis_grp %>% forcats::as_factor(levels=c("resistant","resilien","no resilient"))) %>%
  ggplot() + 
  geom_errorbar(aes(x=forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric) , ymin=genus_min+10^-6, ymax=genus_max+10^-6, group = forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric), col=hysteresis_grp), width=0) + 
  geom_point(aes(x=forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric) , y=genus_min+10^-6,
                 group = forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric), col=hysteresis_grp)) + 
  
  geom_point(aes(x=forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric) , y=genus_max+10^-6,
                 group = forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric), col=hysteresis_grp)) +
  #facet_wrap(~tipping_point_crosses, scales="free_x") + 
  scale_y_log10() +
  geom_hline(yintercept = tip_akk, linetype=2) +
  xlab("Rat id") +
  ylab("Akkermansia\nrel. abundance") +
  danr::theme_dan() +
  theme(axis.text.y = element_text(angle = 0, vjust = 0.5, size=6)) +
  coord_flip() +
  guides(col=FALSE) -> plot_tipping_point_crosses




```


```{r, fig.height=4, fig.width=5.5}

akk_tipping_retention %>%
  mutate(hysteresis_grp = hysteresis_grp %>% forcats::as_factor(levels=c("resistant","resilien","no resilient"))) %>%
  ggplot() + 
  geom_errorbar(aes(x=forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric) , ymin=B.q0.025, ymax=B.q0.975, group = forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric), col=hysteresis_grp)) + 
  geom_point(aes(x=forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric) , y=Boot.Estimate, group = forcats::fct_reorder(Rat_ID,hysteresis_grp %>% as.factor %>% as.numeric), col=hysteresis_grp)) +
  #facet_wrap(~retention, scales="free_x") + 
  #scale_y_log10() +
  geom_hline(yintercept = 0, linetype=2) +
  xlab("Rat id") +
  ylab("Hysteresis\nretention parameter") +
  danr::theme_dan() +
  theme(axis.text.y = element_text(angle = 0, vjust = 0.5, size=6)) +
  coord_flip() -> plot_hysteresis_retention





```


```{r, fig.height=3, fig.width=12}

cowplot::plot_grid(
  pertubation_press_plot + danr::theme_dan(),
  plot_tipping_point_crosses,
  plot_hysteresis_retention + scale_color_discrete("detected profile"), 
  rel_widths = c(6,4,8), nrow = 1, labels="AUTO")




```

## hysteresis models



```{r, fig.height=4, fig.width=12}


models_all_summary$values %>%
  merge( akk_tipping_retention %>%
  select(Rat_ID, hysteresis_grp), by.x="Subject", by.y="Rat_ID") %>%
  select(hysteresis_grp,Parameter,Boot.Estimate) %>%
  mutate(Boot.Estimate = ifelse(Parameter=="retention", abs(Boot.Estimate), Boot.Estimate)) %>%
  group_by(hysteresis_grp,Parameter) %>%
  summarise(m=mean(Boot.Estimate), sd=sd(Boot.Estimate))  %>% filter(hysteresis_grp=="no resilient")


rbind(
  
mloop(retention=0,n.points=100,period=99, m=1,n=3, b.y=0, b.x=15.4778485, cx = 10.9715306, cy=-0.9264833, phase.angle = -7.3910444) %>% with(., data.frame(x=.$x,y=.$y, hysteresis_grp="resistant")),


mloop(retention=0,n.points=100,period=99, m=1,n=3, b.y=-0.6499827, b.x=16.5271817, cx = 11.8892321, cy=-1.7996530, phase.angle = -7.3336891) %>% with(., data.frame(x=.$x,y=.$y, hysteresis_grp="resilient")),

mloop(retention=1.3664291,n.points=100,period=99, m=1,n=3, b.y=-0.5516496, b.x=15.6476289, cx = 11.1731736, cy=-1.9507013, phase.angle = -7.4588586) %>% with(., data.frame(x=.$x,y=.$y, hysteresis_grp="no resilient"))

) -> model_hysteresis_df

model_hysteresis_df %>%
  ggplot() + geom_path(aes(x=x,y=y)) + facet_wrap(~hysteresis_grp) +
  danr::theme_dan() + 
  xlab("Press perturbation\n(days)") + ylab("Akkermansia\nrel. abundance")


```



```{r}

akk_df %>%
  filter(Study_group %in% c("G5","G6","G7")) %>%
  filter(!(Rat_ID == "R67" & Time == "63")) %>%
  filter(!(Rat_ID == "R66" & Time == "68")) %>%
  select(-A1,-A2,-A3,-Study_group,-diet) %>%
  select(genus,press_cumul_day,Rat_ID,hysteresis_grp, Time) %>%
  mutate(genus_norm = log10(genus+10^-6) , press_cumul_day = press_cumul_day %>% as.character %>% as.numeric() ) %>%
  mutate(Rat_ID = Rat_ID %>% as.character) %>%
  select(-hysteresis_grp) %>%
  merge( akk_tipping_retention %>%
  select(Rat_ID, hysteresis_grp),by="Rat_ID") %>%
  mutate(hysteresis_grp = hysteresis_grp %>% forcats::as_factor(levels=c("resistant","resilien","no resilient"))) %>%
  group_by(hysteresis_grp,press_cumul_day,Time) %>%
  mutate(genus_log=log10(genus+10^-6)) %>%
  summarise(m=mean(genus_log), sd=sd(genus_log), n=n()) %>%
  mutate(sem=sd/sqrt(n)) %>%
  mutate(CI_inf=m-1.96*sem, CI_sup=m+1.96*sem) -> observed_hystesis_data_df

 
observed_hystesis_data_df %>%

 ggplot() + 
  geom_point(aes(x=press_cumul_day,y=m, col=Time>30)) +
  ggrepel::geom_text_repel(aes(x=press_cumul_day,y=m, label=Time, col=Time>30)) +
  geom_errorbar(aes(x=press_cumul_day,ymin=CI_inf, ymax=CI_sup, col=Time>30), width=0.4) +
  facet_wrap(~hysteresis_grp) + danr::theme_dan() +
  ylab("Akkermansia\nrel. abundance") + xlab("Press perturbation") +
  guides(col=FALSE) +
  annotate("segment", x = 0, xend = 4, y = -4, yend = -4, size=4, colour = "grey") + 
  annotate("text", x = 10, y = -4,  label = "model")





```


```{r}


observed_hystesis_data_df %>%

 ggplot() + 
  geom_point(aes(x=press_cumul_day,y=m, col=Time>31)) +
  ggrepel::geom_text_repel(aes(x=press_cumul_day,y=m, label=Time, col=Time>31)) +
  geom_errorbar(aes(x=press_cumul_day,ymin=CI_inf, ymax=CI_sup, col=Time>31), width=0.4) +
  facet_wrap(~hysteresis_grp) + danr::theme_dan() +
  xlab("Akkermansia\nrel. abundance") + ylab("Press perturbation") +
  guides(col=FALSE) +
  geom_path(data=model_hysteresis_df, aes(x=x,y=y), col="grey", size=3, alpha=0.5)  +
  danr::theme_dan() + 
  xlab("Press perturbation\n(days)") + ylab("Akkermansia\nrel. abundance (Log10 scale)") -> plot_hysteresis_model



```


### final plot
```{r, fig.height=8, fig.width=12}

dat_text = data.frame(
  label = c("Hysteresis model fit", "", ""),
  hysteresis_grp   = c("resistant", "resilient", "no resilient"),
  x = c(-9,NA,NA),
  y = c(-4,NA,NA),
  xend = c(-7,NA,NA),
  yend = c(-4,NA,NA)
)



cowplot::plot_grid(
cowplot::plot_grid(
  pertubation_press_plot + danr::theme_dan()  +
  theme(legend.position = c(0.65,0.75)),
  plot_tipping_point_crosses + ylab("Akkermansia rel. abund.\n and Tipping point") + scale_color_brewer("detected profile", type = "qual"),
  plot_hysteresis_retention + scale_color_brewer("detected profile", type = "qual") + theme(legend.position = c(0.45, 0.25)), 
  rel_widths = c(6,4,4), nrow = 1, labels="AUTO"),
plot_hysteresis_model +  
  scale_color_manual("", values=c("red","blue"), labels=c("Chronic challenge (DSS)","Washout period (Water)")) +
  theme(legend.position = c(0, 0.25)) +
  #annotate("segment", x = -9, xend = -7, y = -4, yend = -4, size=3, colour = "grey") + 
  geom_segment(data=dat_text, mapping = aes(x = x, xend = xend, y = y, yend = yend), size=3, colour = "grey") +
  geom_text(data=dat_text, mapping = aes(x = -6, y = -4, label = label, hjust = 0)) +
  #annotate("text", x = -1, y = -4,  label = "model") +
  guides(colour = guide_legend("Treatment", override.aes = list(size = 1, alpha = 1))), nrow=2, labels=c(NA,"D"))




ggsave("Figure3_hysteresis_modeling.pdf", h=8, w=12)

```



```{r}


data(EllipseData)
EllipseData
models <- fel(EllipseData$X,EllipseData$Y,method="harmonic2",subjects=EllipseData$subjects,subset=EllipseData$repeated==1)
models

summary(models)


loopf <- mloop(sd.x=0.07,sd.y=0.05,n=3,m=3, retention=.5)
loopf.model <- floop(loopf$x,loopf$y,n=3,m=3)


for (i in c(-1,1,0,-4,4)){
  obj<-mloop(retention=i,n.points=100,period=99, m=3,n=3)
  plot(obj$x,obj$y,xlim=c(-1,1),ylim=c(-5,5),type="l",main=paste("retention=",i,sep=""),xlab="x",ylab="y")
  segments(0,0,0,i)
  text(0.3,0.5,"Retention")
}

```