---
title: "Loop"
author: "Susi/Felix"
date: "5th Dec 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

```

```{r}
getwd()
```

```{r}
library(readr)
library(dplyr)
library(metafor)
library(devtools)
library(patchwork)
library(purrr)
library(tidyverse)
library(tibble)
```

#save in a folder called "export"

```{r}
source("D:/Susi/garvan/Github/mice_sex_diff/mice_sex_diff_Nov/R/meta_analysis_grouping.R")
source("D:/Susi/garvan/Github/mice_sex_diff/mice_sex_diff_Nov/R/data_load_clean_grouping.R")
source("D:/Susi/garvan/Github/mice_sex_diff/mice_sex_diff_Nov/R/calc_pop_stats_grouping.R")
```

# Prepare data 
```{r ignore}
# Load raw data - save cleaned dataset as RDS for reuse
#data_raw <- load_raw("data/dr7.0_all_control_data.csv") %>%
#    clean_raw_data()
#dir.create("export", F, F)
#saveRDS(data_raw, "export/data_clean.rds")

data1 <- readRDS("export/data_clean.rds")
```


CLEAN DATA: select traits with at least 2 centers
```{r ignore}
dat1 <-
  data1 %>%
  group_by(parameter_name) %>% # am keeping this here - not based on parameter groups
  summarize(center_per_trait = length(unique(production_center, na.rm = TRUE)))
#dat1$center_per_trait

```

```{r ignore}
dat2 <- merge(data1, dat1) #sollte auch im grossen Datensatz so funktionieren, also ohne die uebereinstimmende Variable (trait) anzugeben
dat_moreThan1center <-
  dat2  %>%
  filter(center_per_trait >= 2)
```


```{r ignore}
data2 <- dat_moreThan1center
min(data2$center_per_trait) #=2;ok!
```


# Define population variable

```{r ignore}
data3 <- data2 %>%
mutate(population = sprintf("%s-%s", production_center, strain_name))
```
 
 add grouping stuff
```{r}  
 group <- read.csv("ParemeterGrouping.csv") 
 data$parameterGroup <- group$parameter[match(data$parameter_name, group$parameter_name)] # needs to be in certain order, doesn't work if group dataset is first!
 # data4 <- merge(data, group, by="parameter_name") #alternative, works also.

```

# Assign each unique parameter_name (=trait,use trait variable) a unique number ('id')

```{r}
# newly created grouping stuff

names(data)[1] <- "parameter_name"  # this is the original parameter_name  - 232 traits
names(data)[16] <- "parameter_group"    # my new variable, where redundant traits are combined
data <- transform(data, id = match(parameter_name, unique(parameter_name)))

```

```{r}
head(data)

```

```{r}
n <- length(unique(data$id)) # 168 now
n
```

# Create matrix to store results for all traits

```{r}
results.alltraits.grouping <- as.data.frame(cbind(c(1:n), matrix(rep(0, n*13), ncol = 13))) #number of individual results per trait = 10
names(results.alltraits.grouping) <- c("id", "lnCVR", "lnCVR_lower", "lnCVR_upper", "lnCVR_se", "lnVR", "lnVR_lower", "lnVR_upper", "lnVR_se", "lnRR", "lnRR_lower", "lnRR_upper" ,"lnRR_se" , "sampleSize")
```


`


# LOOP
I am removing the variable "weight", which has some missing data, got some error message before that data were excluded below, as there are "NA"s in the rows.
```{r ignore}
#data <- data[,c(1:11,13:16)]
```


```{r}
#write.csv(unique(data$parameter_name, "parametersToCheck.csv"))
```


```{r}

for(t in 1:n) {
  
  tryCatch({
    
    data_par_age <- data_subset_parameterid_individual_by_age(data, t, age_min = 0, age_center = 100)
    
    population_stats <- calculate_population_stats(data_par_age)
    
    results <- create_meta_analysis_effect_sizes(population_stats)
    
#lnCVR,  log repsonse-ratio of the coefficient of variance    
    cvr <- metafor::rma.mv(yi = effect_size_CVR, V = sample_variance_CVR, random = list(~1| strain_name, ~1|production_center,    #~1|parameter_name,
                                                                                        ~1|err), control=list(optimizer="optim", optmethod="Nelder-Mead"), data = results)
    results.alltraits.grouping[t, 2] <- cvr$b
    results.alltraits.grouping[t, 3] <- cvr$ci.lb
    results.alltraits.grouping[t, 4] <- cvr$ci.ub
    results.alltraits.grouping[t, 5] <- cvr$se
    
    cvr
    
    #lnVR, comparison of standard deviations   
    
cv <- metafor::rma.mv(yi = effect_size_VR, V = sample_variance_VR, random = list(~1| strain_name, ~1|production_center, #~1|parameter_name,
                                                                                 ~1|err), control=list(optimizer="optim", optmethod="Nelder-Mead"),  verbose=F, data = results)
    results.alltraits.grouping[t, 6] <- cv$b
    results.alltraits.grouping[t, 7] <- cv$ci.lb
    results.alltraits.grouping[t, 8] <- cv$ci.ub
    results.alltraits.grouping[t, 9] <- cv$se
    
    # for means, lnRR

means <- metafor::rma.mv(yi = effect_size_RR, V = sample_variance_RR, random = list(~1| strain_name, ~1|production_center,
                                                                                    #~1|parameter_name,
                                                                                    ~1|err),
    control=list(optimizer="optim", optmethod="Nelder-Mead"), data = results)
    results.alltraits.grouping[t, 10] <- means$b
    results.alltraits.grouping[t, 11] <- means$ci.lb
    results.alltraits.grouping[t, 12] <- means$ci.ub
    results.alltraits.grouping[t, 13] <- means$se
     
      # for the associated sample sizes, i.e. number of centres contributing ? or combination centes/ strains (doesn't lookm like it)
        results.alltraits.grouping[t, 14] <- means$k
   
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}
 
```

for those traits where convergence wasn't achieved, manual calculation, with 1000 iterations

Trait 125
```{r} 

    data_par_age <- data_subset_parameterid_individual_by_age(data, 125, age_min = 0, age_center = 100)
    
    population_stats <- calculate_population_stats(data_par_age)
    
    results <- create_meta_analysis_effect_sizes(population_stats)
    
#lnCVR,  log repsonse-ratio of the coefficient of variance    
    cvr <- metafor::rma.mv(yi = effect_size_CVR, V = sample_variance_CVR, random = list(~1| strain_name, ~1|production_center,    #~1|parameter_name,
                                                                                        ~1|err), control=list(optimizer="optim", optmethod="Nelder-Mead",  maxit= 1000), data = results)
    results.alltraits.grouping[125, 2] <- cvr$b
    results.alltraits.grouping[125, 3] <- cvr$ci.lb
    results.alltraits.grouping[125, 4] <- cvr$ci.ub
    results.alltraits.grouping[125, 5] <- cvr$se
    
    cvr
    
    #lnVR, comparison of standard deviations   
    
cv <- metafor::rma.mv(yi = effect_size_VR, V = sample_variance_VR, random = list(~1| strain_name, ~1|production_center, #~1|parameter_name,
                                                                                 ~1|err), control=list(optimizer="optim", optmethod="Nelder-Mead", maxit= 1000),  verbose=F, data = results)
    results.alltraits.grouping[125, 6] <- cv$b
    results.alltraits.grouping[125, 7] <- cv$ci.lb
    results.alltraits.grouping[125, 8] <- cv$ci.ub
    results.alltraits.grouping[125, 9] <- cv$se
    
    # for means, lnRR

means <- metafor::rma.mv(yi = effect_size_RR, V = sample_variance_RR, random = list(~1| strain_name, ~1|production_center,
                                                                                    #~1|parameter_name,
                                                                                    ~1|err),
    control=list(optimizer="optim", optmethod="Nelder-Mead", maxit= 1000), data = results)
    results.alltraits.grouping[125, 10] <- means$b
    results.alltraits.grouping[125, 11] <- means$ci.lb
    results.alltraits.grouping[125, 12] <- means$ci.ub
    results.alltraits.grouping[125, 13] <- means$se
     
      # for the associated sample sizes, i.e. number of centres contributing ? or combination centes/ strains (doesn't lookm like it)
        results.alltraits.grouping[125, 14] <- means$k
   


```
  
  
  to adjust for traits 8,18,84,125,132, 144 (doesn't work), 155, 231 (doesn't work)
  
```{r}
data_par_age <- data_subset_parameterid_individual_by_age(data, 8, age_min = 0, age_center = 100)
    
    population_stats <- calculate_population_stats(data_par_age)
    
    results <- create_meta_analysis_effect_sizes(population_stats)
    
#lnCVR,  log repsonse-ratio of the coefficient of variance    
    cvr <- metafor::rma.mv(yi = effect_size_CVR, V = sample_variance_CVR, random = list(~1| strain_name, ~1|production_center,    #~1|parameter_name,
                                                                                        ~1|err), control=list(optimizer="optim", optmethod="Nelder-Mead",  maxit= 1000), data = results)
    results.alltraits.grouping[8, 2] <- cvr$b
    results.alltraits.grouping[8, 3] <- cvr$ci.lb
    results.alltraits.grouping[8, 4] <- cvr$ci.ub
    results.alltraits.grouping[8, 5] <- cvr$se
    
    cvr
    
    #lnVR, comparison of standard deviations   
    
cv <- metafor::rma.mv(yi = effect_size_VR, V = sample_variance_VR, random = list(~1| strain_name, ~1|production_center, #~1|parameter_name,
                                                                                 ~1|err), control=list(optimizer="optim", optmethod="Nelder-Mead", maxit= 1000),  verbose=F, data = results)
    results.alltraits.grouping[8, 6] <- cv$b
    results.alltraits.grouping[8, 7] <- cv$ci.lb
    results.alltraits.grouping[8, 8] <- cv$ci.ub
    results.alltraits.grouping[8, 9] <- cv$se
    
    # for means, lnRR

means <- metafor::rma.mv(yi = effect_size_RR, V = sample_variance_RR, random = list(~1| strain_name, ~1|production_center,
                                                                                    #~1|parameter_name,
                                                                                    ~1|err),
    control=list(optimizer="optim", optmethod="Nelder-Mead", maxit= 1000), data = results)
    results.alltraits.grouping[8, 10] <- means$b
    results.alltraits.grouping[8, 11] <- means$ci.lb
    results.alltraits.grouping[8, 12] <- means$ci.ub
    results.alltraits.grouping[8, 13] <- means$se
     
      # for the associated sample sizes, i.e. number of centres contributing ? or combination centes/ strains (doesn't lookm like it)
        results.alltraits.grouping[8, 14] <- means$k
```

 I still get the NA warnings / removals, but have extensively checked the data set - there are no more NAs present.
 
 I also don't understand the k<= 1, as traits that were collected by one center only had been excluded in the previous steps.
 Also, "Single-level factor(s) found in 'random' argument" is strange, as some trais that actually worked and were included have been done in one strain only. so I don't understand what's egtting thrown out here.
 # due to "parameter_name" being included as random factor; which  doesn't really work because for the loop the analyses are within each parameter name? 
 also, number of iterations has beeen increased to 500, so that those that disn't converge at 100 interations may be included.
 (however, this significantly affects running time! takes very long...)
 
 Next step is to clean the data, i.e. remove any traits that don't make sense (i.e. number of ribs / digits , numbers of events etc)
 
 
```{r}
procedures <- read.csv("procedures.csv")
```

```{r}  
 results.alltraits.grouping$parameter_group <- data$parameter_group[match(results.alltraits.grouping$id, data$id)]
 results.alltraits.grouping$procedure <- data$procedure_name[match(results.alltraits.grouping$id, data$id)]
```


```{r}   
results.alltraits.grouping$GroupingTerm <-  procedures$GroupingTerm[match(results.alltraits.grouping$procedure, procedures$procedure)]



```

check - remove traits #  144,158,160,161,162,163,165,166,167,168, 221,222,231
```{r}
meta1 <- results.alltraits.grouping
```

```{r}
meta1b <- meta1[ !(meta1$id %in% c(144,158,160,161,162,163,165,166,167,168, 221,222,231)), ]

meta1b$parameter_name <- data$parameter_name[match(meta1b$id, data$id)]

```

```{r}
write.csv(meta1b, "meta1.csv")
```




## COLLAPSING DOWN CORRELATED  PARAMETER_NAMES


```{r}
meta1 <- read.csv("meta.csv") # cleaned up data/ cleaned up parameter groups. 
unique(meta1$procedure) #19. Can still get reduced.
unique(meta1$GroupingTerm) #10. sounds reasonable
unique(meta1$parameter_group) # 152 levels. To be used as random factor in meta-meta!
##  test for meta-meta-analysis / collapse confounded parameters!
```


```{r}
meta1 <- read.csv("meta.csv") # cleaned up data/ cleaned up parameter groups. 
unique(meta1$procedure) #19. Can still get reduced.
unique(meta1$GroupingTerm) #10. sounds reasonable
unique(meta1$parameter_group) # 152 levels. To be used as random factor in meta-meta!
```

```{r meta-meta}
# use meta1 , newly loaded (as has been cleaned / checked in excel)

# need to write function to collapse the correlated parametrs down (based on parameter_groups).  For the parameter_groups that are represented more than 1, run meta analysis, the other traits stay as are.  rma(yi = estimate i.e. ln RR / ln CVR / ln VR, se = lnRR_se / lnCVR_se / lnVR_se, method = "FE" )


meta1[,"MetaMeanV"] <-meta1$lnRR_se^2
meta1[,"MetaCVRV"] <-meta1$lnCVR_se^2
meta1[,"MetaVRV"] <-meta1$lnVR_se^2
```


workflow:
count parameter_names per parameter_group, save in new list "par_group_size":
```{r}
meta1 %>%  count(parameter_group) # 152
```
dataframe: most groups have only 1 , it;s only a couple where there's multiple parameter names.


```{r}
meta1b <-
  meta1 %>%
  group_by(parameter_group) %>% 
  summarize(par_group_size = length(unique(parameter_name, na.rm = TRUE)))
#this gives me a summary of number of parameter names in each parameter group, now I neeed to merge it back together


meta1$par_group_size <- meta1b$par_group_size[match(meta1$parameter_group, meta1b$parameter_group)]

```

create subsets: needs to use those that have more than 1 count (par_group_size)


```{r}
meta1_sub<-subset(meta1,par_group_size >1) # 86 observations
```

NEST
```{r}
n_count <- meta1_sub %>% 
  group_by(parameter_group) %>% 
  mutate(raw_N = sum(sampleSize)) %>%  # get the sample sizes as well...??
  
  nest()
```


FIXED EFFECTS MODEL

#temporary: SE is currently *count, but really should be using the covariance matrix - CAN BE DONE IN  METAFOR??
```{r}
model_count <- n_count %>% 
  mutate(model_lnRR = map(data, ~ rma(yi = lnRR, sei= lnRR_se*par_group_size, method = "FE", data = .)),
        model_lnCVR = map(data, ~ rma(yi = lnCVR, sei= lnCVR_se*par_group_size, method = "FE", data = .)),
        model_lnVR = map(data, ~ rma(yi = lnVR, sei= lnVR_se*par_group_size, method = "FE", data = .)))


 # OLD, to check of ses have changed
sub_model3 <- n_count %>%
  mutate(model_lnRR= map(data, ~ rma(yi = lnRR, sei = lnRR_se, method = "FE", data= .)),
  model_lnCVR= map(data, ~ rma(yi = lnCVR, sei = lnCVR_se, method = "FE", data= .)),
  model_lnVR= map(data, ~ rma(yi = lnVR, sei = lnVR_se, method = "FE", data= .)))


```



 unknot this somehow
 
```{r}

count_fun <- function(mod_sub)
  coef(summary(mod_sub))[1, c(1,5,6, 2)]

sub_results <- model_count %>% 
  transmute(parameter_group, estimatelnCVR = map(model_lnCVR, count_fun),
            parameter_group, estimatelnVR = map(model_lnVR, count_fun),
            parameter_group, estimatelnRR = map(model_lnRR, count_fun)) %>%
   unnest()

names(sub_results) <- c("parameter_group","lnCVR","lnCVR_lower","lnCVR_upper","lnCVR_se","lnVR",
				"lnVR_lower","lnVR_upper","lnVR_se", "lnRR","lnRR_lower","lnRR_upper","lnRR_se")

########################### seems to be working. 

sub3_results <- sub_model3 %>% 
  transmute(parameter_group, estimatelnCVR = map(model_lnCVR, count_fun),
            parameter_group, estimatelnVR = map(model_lnVR, count_fun),
            parameter_group, estimatelnRR = map(model_lnRR, count_fun)) %>%
   unnest()

#print(sub_model$model)

names(sub3_results) <- c("parameter_group","lnCVR","lnCVR_lower","lnCVR_upper","lnCVR_se","lnVR",
				"lnVR_lower","lnVR_upper","lnVR_se", "lnRR","lnRR_lower","lnRR_upper","lnRR_se")

# ok, what is the sample size here??? just the counts of the combined data? But the sample sizes of the single parameters is much higher then???? n_count contains the raw N as well, but don't kno how to extract it


```


merge the two (the new and the initial) datasets back together

```{r}
#5. Merge with set with count = 1
meta_all <- meta1 %>% filter(par_group_size == 1) %>% as_tibble

str(meta_all)
str(sub_results)

col.diff<- which(is.na(match(names(meta_all),names(sub_results))))
```

# combinding
```{r}
combined<-bind_rows(sub_results, meta_all)
glimpse(combined)

```

```{r}
combo <- combined
combo$counts <- meta1$par_group_size[match( combo$parameter_group, meta1$parameter_group)]

combo$procedure2 <-meta1$procedure[match( combo$parameter_group, meta1$parameter_group)]
combo$GroupingTerm2 <-meta1$GroupingTerm[match( combo$parameter_group, meta1$parameter_group)]


```

```{r}
str(combo)
head(combo)
combo <-combo[, c(1:13, 15, 16, 24:26)]
combo <-combo[, c(1, 14:18, 2:13)]
#rename
colnames(combo)[which(names(combo) == c("procedure2", "GroupingTerm2"))] <- c("procedure", "GroupingTerm")
```
# save combo, this is to be used for graphs.

```{r}
#write.csv(combo, "DataMetaCombined.csv"))
```
### PLOTS

```{r}
meta_all <- combo
str(meta_all)
head(meta_all)
```
```{r}
table(meta_all$GroupingTerm)
```


######################################################################
#Create 6 plots of lnCVR, lnCV, lnRR, separate for males and females #

# 1. Males biased (lnCVR, lnCV, lnRR > 0) #

```{r}
meta_all0 <- meta_all[, c("lnCVR", "lnVR", "lnRR", "GroupingTerm")]  
```

```{r}
meta_all1 <- gather(meta_all0, trait, value, lnCVR:lnRR)
meta_all1 <- with(meta_all1, meta_all1[order(trait, GroupingTerm),])

meta_all_malebias0 <- meta_all1 %>%
				group_by_at(vars(trait, GroupingTerm)) %>%
				summarise(malebias0= sum(value > 0), femalebias0= sum(value<= 0), total= sum(value<1000), 
					malepercent= malebias0*100/total, femalepercent= femalebias0*100/total)
print.data.frame(meta_all_malebias0)

malebias0_plot <- ggplot(meta_all_malebias0, aes(GroupingTerm, malepercent)) + 
    	geom_bar(aes(fill = trait), stat= "identity", position= "dodge") +   
	ylab("Percentage of parameters \nwith male bias")  + 
	ylim(0, 100) +
	ggtitle("Male biased sex difference (trait > 0)") + 
	theme_bw(base_size = 18) +
	theme(axis.title.x = element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
		axis.text.x = element_text(angle = -30), panel.border=element_blank(), axis.line = element_line(colour = "black") ) 
malebias0_plot 


```



#--Save plot to pdf:
```{r}
ggsave(malebias0_plot , file="malebias0_plot .pdf", width= 10, height= 5)
```



# 2. Females biased (lnCVR, lnCV, lnRR > 0) #

```{r}
meta_all_femalebias0 <- meta_all1 %>%
				group_by_at(vars(trait, GroupingTerm)) %>%
				summarise(malebias0= sum(value > 0), femalebias0= sum(value<= 0), total= sum(value<1000), 
					malepercent= malebias0*100/total, femalepercent= femalebias0*100/total)
print.data.frame(meta_all_femalebias0)

femalebias0_plot <- ggplot(meta_all_femalebias0, aes(GroupingTerm, femalepercent)) + 
    	geom_bar(aes(fill = trait), stat= "identity", position = "dodge") +
	ylab("Percentage of parameters \nwith female bias")  + 
	ylim(0, 100) +
	ggtitle("Female biased sex difference (trait < 0)") + 
	theme_bw(base_size = 18) +
	theme(axis.title.x = element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
		axis.text.x = element_text(angle = -30), panel.border=element_blank(), axis.line = element_line(colour = "black") ) 
femalebias0_plot 

```

#--Save plot to pdf:
```{r}
ggsave(femalebias0_plot , file="femalebias0_plot .pdf", width= 10, height= 5)

```

###########################################
# 3. Male biased ( delta lnCVR, lnCV, lnRR > 10%) #
 (as above)
meta_all0 <- meta_all[, c("lnCVR", "lnVR", "lnRR", "GroupingTerm")]  #"sampleSize", 

meta_all1 <- gather(meta_all0, trait, value, lnCVR:lnRR)
meta_all1 <- with(meta_all1, meta_all1[order(trait, GroupingTerm),])

#choosing male biased ln-ratios of a larger difference in untransformed m/f ratios than 10%
```{r}
meta_all_malebias10perc <- meta_all1 %>%
				group_by_at(vars(trait, GroupingTerm)) %>%
				summarise(malebias10perc= sum(value > log(11/10)), femalebias10perc= sum(value < -log(11/10)), total= sum(value<1000), 
					malepercent= malebias10perc*100 / total, femalepercent= femalebias10perc*100 / total)
print.data.frame(meta_all_malebias10perc)

malebias10perc_plot <- ggplot(meta_all_malebias10perc, aes(GroupingTerm, malepercent)) + 
    	geom_bar(aes(fill = trait), stat= "identity", position = "dodge") +
	ylab("Percentage of parameters \nwith male bias")  + 
	ylim(0, 100) +
	ggtitle("Male biased sex difference (> 10%)") + 
	theme_bw(base_size = 18) +
	theme(axis.title.x = element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
		axis.text.x = element_text(angle = -30), panel.border=element_blank(), axis.line = element_line(colour = "black") ) 
malebias10perc_plot #missing bars represent 0%


#--Save plot to pdf:
ggsave(malebias10perc_plot , file="malebias10perc_plot .pdf", width= 10, height= 5)

```



###########################################
# 4. Female biased ( delta lnCVR, lnCV, lnRR > 10%) #

#choosing male biased ln-ratios of a larger difference in untransformed m/f ratios than 10%
```{r}
meta_all_malebias10perc <- meta_all1 %>%
				group_by_at(vars(trait, GroupingTerm)) %>%
				summarise(malebias10perc= sum(value > log(11/10)), femalebias10perc= sum(value < -log(11/10)), total= sum(value<1000), 
					malepercent= malebias10perc*100 / total, femalepercent= femalebias10perc*100 / total)
print.data.frame(meta_all_malebias10perc)

femalebias10perc_plot <- ggplot(meta_all_malebias10perc, aes(GroupingTerm, femalepercent)) + 
    	geom_bar(aes(fill = trait), stat= "identity", position = "dodge") +
	ylab("Percentage of parameters \nwith female bias")  + 
	ylim(0, 100) +
	ggtitle("Female biased sex difference (> 10%)") + 
	theme_bw(base_size = 18) +
	theme(axis.title.x = element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
		axis.text.x = element_text(angle = -30), panel.border=element_blank(), axis.line = element_line(colour = "black") ) 
femalebias10perc_plot #missing bars represent 0%

#--Save plot to pdf:
ggsave(femalebias10perc_plot , file="femalebias10perc_plot.pdf", width= 10, height= 5)
```




###########################################
# 5. Male biased ( delta lnCVR, lnCV, lnRR different from zero) #

#create column with 1= different from zero, 0= zero included in CI

```{r}
meta_all_sig <- meta_all %>% 
    mutate(lnCVRsig = ifelse(lnCVR_lower*lnCVR_upper >0, 1, 0), lnVRsig = ifelse(lnVR_lower*lnVR_upper >0, 1, 0), 
	     lnRRsig = ifelse(lnRR_lower*lnRR_upper > 0, 1,0))

meta_all_sig2 <- meta_all_sig[, c("lnCVR", "lnVR", "lnRR", "lnCVRsig", "lnVRsig", "lnRRsig", "GroupingTerm")]  #"sampleSize", 

meta_all_sig3 <- gather(meta_all_sig2, trait, value, lnCVR:lnRR)
meta_all_sig4 <- with(meta_all_sig3, meta_all_sig3[order(trait, GroupingTerm),])
meta_all_sig4$sig <- "placeholder"

meta_all_sig4$sig <- ifelse(meta_all_sig4$trait == "lnCVR", meta_all_sig4$lnCVRsig,
				ifelse(meta_all_sig4$trait == "lnVR", meta_all_sig4$lnVRsig, meta_all_sig4$lnRRsig))

```

#choosing sex biased ln-ratios significantly larger than 0
```{r}
meta_all_malebiasSig0 <- meta_all_sig4 %>%
				group_by_at(vars(trait, GroupingTerm)) %>%
				filter(sig== 1) %>%
				summarise(male_sig= sum(value > 0), female_sig= sum(value < 0)) 
meta_all_malebiasSig1 <- ungroup(meta_all_malebiasSig0) %>%
				add_row(trait = "lnCVR", GroupingTerm = "Hearing", male_sig = 0, female_sig= 0, .before = 4) %>% #add "Hearing" for lnCVR (not filtered as only zeros)
				mutate(total = meta_all_malebias0$total, malepercent= male_sig*100 / total, femalepercent= female_sig*100 / total)

print.data.frame(meta_all_malebiasSig1) #check

malebiasSig1_plot <- ggplot(meta_all_malebiasSig1, aes(GroupingTerm, malepercent)) + 
    	geom_bar(aes(fill = trait), stat= "identity", position = "dodge") +
	ylab("Percentage of parameters \nwith male bias")  + 
	ylim(0, 100) + 
	ggtitle("Male biased sex difference (significant)") + 
	theme_bw(base_size = 18) +
	theme(axis.title.x = element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
		axis.text.x = element_text(angle = -30), panel.border=element_blank(), axis.line = element_line(colour = "black") ) 
malebiasSig1_plot #missing bars represent 0%

#--Save plot to pdf:
ggsave(malebiasSig1_plot , file="malebiasSig1_plot.pdf", width= 10, height= 5)
```


###########################################
# 6. Female biased ( delta lnCVR, lnCV, lnRR > 10%) #

```{r}
femalebiasSig1_plot <- ggplot(meta_all_malebias10perc, aes(GroupingTerm, femalepercent)) + 
    	geom_bar(aes(fill = trait), stat= "identity", position = "dodge") +
	ylab("Percentage of parameters \nwith female bias")  +
	ylim(0, 100) + 
	ggtitle("Female biased sex difference (significant)") + 
	theme_bw(base_size = 18) +
	theme(axis.title.x = element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
		axis.text.x = element_text(angle = -30), panel.border=element_blank(), axis.line = element_line(colour = "black") ) 
femalebiasSig1_plot #missing bars represent 0%

#--Save plot to pdf:
ggsave(femalebiasSig1_plot , file="femalebiasSig1_plot.pdf", width= 10, height= 5)


```