Skip to content

Commit

Permalink
Add cowres results and finalise simres
Browse files Browse the repository at this point in the history
  • Loading branch information
haziqj committed Aug 11, 2022
1 parent 78d7344 commit f94b8be
Show file tree
Hide file tree
Showing 8 changed files with 236 additions and 3 deletions.
19 changes: 19 additions & 0 deletions 01-cow.R
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ cattle$id <- as.factor(cattle$id) # convert to factors
# levels(cattle$group) <- c("Treatment A", "Treatment B")
# str(cattle)
# write_csv(cattle, file = "cattle.csv")
N <- nrow(cattle)

## ---- plot_cow_data --------
ggplot(cattle, aes(time, weight, group = id, col = group)) +
Expand Down Expand Up @@ -67,12 +68,30 @@ est_all <- function(kernel = "fbm,0.5", est.hurst = FALSE) {

# Hurst = 0.5
res1 <- est_all(kernel = "fbm,0.5", est.hurst = FALSE)
res1 <- res1 %>%
mutate(
k = c(2, 3, 3, 4, 4),
AIC = 2 * k - 2 * loglik,
BIC = k * log(N) -2 * loglik
)

# Hurst = 0.3
res2 <- est_all(kernel = "fbm,0.3", est.hurst = FALSE)
res2 <- res2 %>%
mutate(
k = c(2, 3, 3, 4, 4),
AIC = 2 * k - 2 * loglik,
BIC = k * log(N) -2 * loglik
)

# Estimate hurst
res3 <- est_all(kernel = "fbm", est.hurst = TRUE)
res3 <- res3 %>%
mutate(
k = c(3, 4, 4, 5, 5),
AIC = 2 * k - 2 * loglik,
BIC = k * log(N) -2 * loglik
)

save(res1, res2, res3, file = "cowres.RData")

Expand Down
19 changes: 19 additions & 0 deletions README.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -95,11 +95,30 @@ res_tab1[, -1] %>%
<<plot_cow_data>>
```

The model fitted was an I-prior model with ANOVA kernel (Pearson & fBm).
The results are tabulated below.

### Fixed hurst = 0.5

```{r}
load("cowres.RData")
mutate(res1, formula = paste0("`", formula, "`")) %>%
kbl(., format = "pipe", digits = 2)
```

### Fixed hurst = 0.3

```{r}
mutate(res2, formula = paste0("`", formula, "`")) %>%
kbl(., format = "pipe", digits = 2)
```

### Estimated hurst value

```{r}
mutate(res3, formula = paste0("`", formula, "`")) %>%
kbl(., format = "pipe", digits = 2)
```



Expand Down
33 changes: 33 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -93,3 +93,36 @@ The geometric mean
## Functional response model

![](figure/cow_plot-1.png)<!-- -->

The model fitted was an I-prior model with ANOVA kernel (Pearson & fBm).
The results are tabulated below.

### Fixed hurst = 0.5

| model | formula | loglik | error | lambda | psi | hurst | k | AIC | BIC |
|------:|:----------------------|---------:|------:|:---------------------|:--------|------:|----:|--------:|--------:|
| 1 | `time` | -2789.23 | 16.25 | 0.837 | 0.00375 | 0.5 | 2 | 5582.46 | 5591.45 |
| 2 | `group * time` | -2789.20 | 16.24 | 0.019,-0.836 | 0.00375 | 0.5 | 3 | 5584.40 | 5597.88 |
| 3 | `id * time` | -2295.16 | 2.89 | -0.203,-0.088 | 0.07384 | 0.5 | 3 | 4596.33 | 4609.81 |
| 4 | `(group + id) * time` | -2270.85 | 2.62 | -1.019,-0.187,-0.085 | 0.08711 | 0.5 | 4 | 4549.70 | 4567.67 |
| 5 | `group * id * time` | -2249.00 | 3.09 | -1.057,4.918,0.047 | 0.06538 | 0.5 | 4 | 4506.00 | 4523.97 |

### Fixed hurst = 0.3

| model | formula | loglik | error | lambda | psi | hurst | k | AIC | BIC |
|------:|:----------------------|---------:|------:|:--------------------|:--------|------:|----:|--------:|--------:|
| 1 | `time` | -2792.78 | 16.22 | 4.465 | 0.00375 | 0.3 | 2 | 5589.56 | 5598.54 |
| 2 | `group * time` | -2792.73 | 16.20 | 0.03,-4.462 | 0.00376 | 0.3 | 3 | 5591.47 | 5604.94 |
| 3 | `id * time` | -2266.39 | 1.62 | -0.163,-0.381 | 0.13445 | 0.3 | 3 | 4538.79 | 4552.26 |
| 4 | `(group + id) * time` | -2242.30 | 1.44 | 0.708,-0.152,-0.36 | 0.15896 | 0.3 | 4 | 4492.60 | 4510.57 |
| 5 | `group * id * time` | -2238.78 | 2.16 | -1.184,-1.265,0.248 | 0.09450 | 0.3 | 4 | 4485.56 | 4503.53 |

### Estimated hurst value

| model | formula | loglik | error | lambda | psi | hurst | k | AIC | BIC |
|------:|:----------------------|---------:|------:|:------------------|:--------|------:|----:|--------:|--------:|
| 1 | `time` | -2788.77 | 16.28 | 0.347 | 0.00374 | 0.62 | 3 | 5583.53 | 5597.01 |
| 2 | `group * time` | -2788.75 | 16.27 | 0.013,-0.35 | 0.00375 | 0.61 | 4 | 5585.49 | 5603.46 |
| 3 | `id * time` | -2253.21 | 0.16 | -0.065,0.83 | 1.24112 | 0.17 | 4 | 4514.43 | 4532.40 |
| 4 | `(group + id) * time` | -2231.13 | 0.13 | 0.102,0.058,0.745 | 1.59394 | 0.18 | 5 | 4472.27 | 4494.73 |
| 5 | `group * id * time` | -2232.78 | 0.18 | 0.11,0.058,0.751 | 1.19639 | 0.18 | 5 | 4475.55 | 4498.01 |
Binary file modified cowres.RData
Binary file not shown.
158 changes: 158 additions & 0 deletions simres_corr.Rmd
Original file line number Diff line number Diff line change
@@ -0,0 +1,158 @@
---
title: "Simulation Results for Linear Models (correlated errors)"
# author: "Haziq"
date: "`r Sys.Date()`"
output: pdf_document
---

```{r setup, include = FALSE}
library(tidyverse)
library(kableExtra)
model_codes <- readxl::read_excel("model_codes.xlsx") %>%
mutate(code = paste0(x1, x2, x3, x1x2, x1x3, x2x3, x1x2x3, sep = ""))
beta_vals <- readxl::read_excel("model_codes.xlsx", sheet = "beta") %>%
mutate(code = paste0(
ifelse(x1 > 0, 1, 0),
ifelse(x2 > 0, 1, 0),
ifelse(x3 > 0, 1, 0),
ifelse(x1x2 > 0, 1, 0),
ifelse(x1x3 > 0, 1, 0),
ifelse(x2x3 > 0, 1, 0),
ifelse(x1x2x3 > 0, 1, 0),
sep = "") %>% as.character()
)
res_cleanup <- function(x) {
# mutate(x, no = ifelse(prop < 0.05, NA, no)) %>%
# drop_na()
x
# x %>%
# drop_na() %>%
# mutate(prop2 = prop,
# prop = prop / sum(prop))
}
res_final_tab <- function(x, true.mod = "1111111") {
the.names <- names(x[c("iprior", "lasso", "spikeslab", "gprior")])
do.call(bind_rows, mapply(cbind, x[c("iprior", "lasso", "spikeslab", "gprior")],
"method" = the.names, SIMPLIFY = F)) %>%
mutate(method = factor(method, levels = the.names)) %>%
filter(mod == true.mod) %>%
group_by(method) %>%
select(mod, prop)
}
```


# Summary

## Results for error sd = 3 and corr = 0.5


```{r, message = FALSE, echo = FALSE}
load("simres.RData")
res_tab <- NULL
for (i in seq_len(length(myres))) {
res_tab <- bind_rows(res_tab,
res_final_tab(myres[[i]], beta_vals[i, 8, drop = TRUE])
)
}
res_tab %>%
pivot_wider(names_from = method, values_from = prop) %>%
replace(is.na(.), 0) %>%
kbl(booktabs = TRUE, position = "h", digits = 2) %>%
kable_styling(position = "center", latex_options = "HOLD_position")
```


\newpage

### Model `r beta_vals[1, 8, drop = TRUE]`

- True value: `r myres[[1]]$beta.true`
- `r paste0("nsim = ", myres[[1]]$nsim, ", n = ", myres[[1]]$n, ", corr = ", myres[[1]]$corr, ", err.sd = ", myres[[1]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[1]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage


### Model `r beta_vals[2, 8, drop = TRUE]`

- True value: `r myres[[2]]$beta.true`
- `r paste0("nsim = ", myres[[2]]$nsim, ", n = ", myres[[2]]$n, ", corr = ", myres[[2]]$corr, ", err.sd = ", myres[[2]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[2]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage


### Model `r beta_vals[3, 8, drop = TRUE]`

- True value: `r myres[[3]]$beta.true`
- `r paste0("nsim = ", myres[[3]]$nsim, ", n = ", myres[[3]]$n, ", corr = ", myres[[3]]$corr, ", err.sd = ", myres[[3]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[3]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage

### Model `r beta_vals[4, 8, drop = TRUE]`

- True value: `r myres[[4]]$beta.true`
- `r paste0("nsim = ", myres[[4]]$nsim, ", n = ", myres[[4]]$n, ", corr = ", myres[[4]]$corr, ", err.sd = ", myres[[4]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[4]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage

### Model `r beta_vals[5, 8, drop = TRUE]`

- True value: `r myres[[5]]$beta.true`
- `r paste0("nsim = ", myres[[5]]$nsim, ", n = ", myres[[5]]$n, ", corr = ", myres[[5]]$corr, ", err.sd = ", myres[[5]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[5]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage


### Model `r beta_vals[6, 8, drop = TRUE]`

- True value: `r myres[[6]]$beta.true`
- `r paste0("nsim = ", myres[[6]]$nsim, ", n = ", myres[[6]]$n, ", corr = ", myres[[6]]$corr, ", err.sd = ", myres[[6]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[6]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage


### Model `r beta_vals[7, 8, drop = TRUE]`

- True value: `r myres[[7]]$beta.true`
- `r paste0("nsim = ", myres[[7]]$nsim, ", n = ", myres[[7]]$n, ", corr = ", myres[[7]]$corr, ", err.sd = ", myres[[7]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[7]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage

### Model `r beta_vals[8, 8, drop = TRUE]`

- True value: `r myres[[8]]$beta.true`
- `r paste0("nsim = ", myres[[8]]$nsim, ", n = ", myres[[8]]$n, ", corr = ", myres[[8]]$corr, ", err.sd = ", myres[[8]]$err.sd)`

```{r echo = FALSE}
lapply(myres[[8]][c("iprior", "lasso", "spikeslab", "gprior")], res_cleanup)
```
\newpage




Binary file modified simres_corr.pdf
Binary file not shown.
10 changes: 7 additions & 3 deletions simres_uncorr.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,13 @@ beta_vals <- readxl::read_excel("model_codes.xlsx", sheet = "beta") %>%
res_cleanup <- function(x) {
# mutate(x, no = ifelse(prop < 0.05, NA, no)) %>%
# drop_na()
x
mutate(x, no = ifelse(prop < 0.05, NA, no)) %>%
drop_na()
# x
# x %>%
# drop_na() %>%
# mutate(prop2 = prop,
# prop = prop / sum(prop))
}
res_final_tab <- function(x, true.mod = "1111111") {
Expand Down
Binary file modified simres_uncorr.pdf
Binary file not shown.

0 comments on commit f94b8be

Please sign in to comment.