---
title: "EPIDEMIOLOGY OF CANINE MAMMARY TUMOURS ON THE CANARY ARCHIPELAGO IN SPAIN"
subtitle: 
header-includes:
   - \usepackage{lmodern}
   - \usepackage[spanish]{babel}
output:
  word_document:
    toc: yes
    toc_depth: 3
  html_document:
    highlight: pygments
    theme: cerulean
    toc: yes
    toc_float: yes
    toc_depth: 5
  pdf_document:
    toc: yes
    toc_depth: 5
    latex_engine: xelatex
fontsize: 12pt 
geometry: margin=1in
mainfont: Times New Roman     
---

\fontsize{13}{16} \selectfont

<style type="text/css">
 .tab{ text-indent: 40px; }
 body{ font-size: 170%; }
 .question{background-color: white; font-size: 110%; font-weight: normal;}
</style>



```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)
library(readxl)
library(tidyverse)
library(janitor)
library(flextable)
library(rstatix)
library(DescTools)
library(ggsci)
library(forestplot)

# Funciones
# ---------

# Formatea P-valores
formatPval <- function(p) ifelse(p<0.0001,"<0.0001",sprintf("%.4f",p))
meanSd <- function(x) sprintf("%.2f ± %.2f",mean(x,na.rm=TRUE),sd(x,na.rm=TRUE))

# Compara la primera y última proporciones de una tabla de Tendencia
# La primera columna es el año, la segunda son los casos de interés y la tercera el
# total de casos:
comparaProp <- function(tablaTendencia){
  dif <- tablaTendencia[c(1,nrow(tablaTendencia)),] 
  names(dif) <- c("Year","x","n")
  test <- with(dif,prop.test(x,n))
  difProp <- 100*diff(test$estimate)
  intervDifProp <- -100*test$conf.int[2:1]
  p <- formatPval(test$p.value)
  pDifCI <- sprintf("%.2f%% (%.2f%%, %.2f%%)",difProp,intervDifProp[1],intervDifProp[2])
  probs <- sprintf("%.2f%%",100*test$estimate)
  result <- data.frame(t(c(probs,pDifCI,p)))
  names(result) <- c(dif$Year,"Difference (IC95%)","p")
  result
}


```


## Reading data.

```{r}
zoocan <- read.csv("zoocan.csv")
zoocan <- zoocan %>% rename(Island=Isle)

all_tumoursFD <- read.csv("all_tumoursFD.csv") 

all_tumoursFD<- all_tumoursFD %>% 
  mutate(tumourType=ifelse(Location=="Mammary gland","CMT","Non-CMT"))

FD <- all_tumoursFD %>% unique()

mammary_tumours <- all_tumoursFD %>% 
  filter(Location=="Mammary gland")

mammary_tumours %>% tabyl(Histo_Cyto)
mammary_tumours %>% tabyl(Benign_Malignant)

```


## CMT by year 

```{r}
# Number of tumours by year, and number of Female Dogs with any tumour by year
all_tumoursFD_year <- all_tumoursFD %>% 
  group_by(ReportRef,ResultDate) %>% 
  summarize(nTumours=n()) %>% 
  group_by(ResultDate) %>% 
  summarize(all_tumoursFD=sum(nTumours),FD_with_any_tumour=n())

# Number of CMT by year, and number of Female Dogs with CMT by year
CMT_FD_year <- all_tumoursFD %>% 
  filter(tumourType=="CMT") %>% 
  group_by(ReportRef,ResultDate) %>% 
  summarize(nTumours=n()) %>% 
  group_by(ResultDate) %>% 
  summarize(CMT_year=sum(nTumours),FD_CMT=n())

# Percentage of CMT and percentage of FD with CMT
pctTum <- full_join(all_tumoursFD_year,CMT_FD_year) %>% 
  mutate(`%FD_with_CMT`=round(100*FD_CMT/FD_with_any_tumour,2),
         `%Tumours_CMT`=round(100*CMT_year/all_tumoursFD,2)) 

pctTum %>% 
  flextable() %>% 
  autofit()

```



```{r}
## Neutering
# Neutering r
Neutering_FD_all_tumoursFD <- all_tumoursFD %>% 
  select(ResultDate, ReportRef, Neuter_status) %>% 
  unique()%>%
  tabyl(ResultDate, Neuter_status) %>% 
  adorn_percentages() %>% 
  select(ResultDate,"Neutering rate all FD"=Neutered) %>%
  adorn_pct_formatting(digits = 2)


Neutering_FD_CMT <-all_tumoursFD %>% 
  filter(Location=="Mammary gland") %>% 
  select(ResultDate, ReportRef, Neuter_status) %>% 
  unique()%>%tabyl(ResultDate, Neuter_status) %>% 
  adorn_percentages() %>% 
  select(ResultDate,"Neutering rate FD_CMT"=Neutered) %>% 
  adorn_pct_formatting(digits = 2)

full_join(Neutering_FD_all_tumoursFD, Neutering_FD_CMT) %>% 
  flextable() %>% 
  autofit()

```

\  

## Table 1: Proportion of pathology reports according to the number of CMT.

```{r}
 all_tumoursFD %>% 
  rename(Year=ResultDate) %>% 
  filter(Location=="Mammary gland")%>% 
  group_by(ReportRef)%>% 
  summarize(nall_tumours=n())%>% 
  mutate(nall_tumours=ifelse(nall_tumours>=5,"≥5",nall_tumours),
         nall_tumours=factor(nall_tumours,levels=c(1:4,"≥5"))) %>% 
  tabyl(nall_tumours)%>% 
  adorn_pct_formatting(digits=2) %>% 
  select(-n) %>%
  pivot_wider(everything(),names_from = nall_tumours,values_from = percent) %>% 
  mutate(`Number of reported CMT nodules`="Proportion of reports") %>% 
  select(6,everything()) %>% 
  flextable() %>% 
  autofit()
```


## Table 2: Relative proportions of the different kind of malignant CMT.

```{r}
bm <- mammary_tumours %>% 
  filter(Benign_Malignant=="Malignant") %>% 
  tabyl(ResultDate, Primary_tumour) %>% 
  adorn_totals("col")

pt <- NULL
for (k in 2:(ncol(bm)-1)){
  pTrend <- bm %>% select(all_of(k),Total) %>% 
    prop_trend_test() %>% 
    pull(p) %>% 
    formatPval
  pt <- c(pt,pTrend)
}

pTrend <- tibble(Tumor=names(bm)[2:(ncol(bm)-1)],pTrend=pt) 

tb <- mammary_tumours %>% 
  filter(Benign_Malignant=="Malignant") %>% 
  tabyl(`Primary_tumour`) %>% 
  adorn_totals("row") %>% 
  rename(Tumor=`Primary_tumour`)

pcis <- NULL
for (k in 1:(nrow(tb)-1)){
  bt <- binom.test(tb$n[k],tb$n[nrow(tb)])
  ci <- bt$conf.int
  pci <- sprintf("%.2f%% (%.2f%%; %.2f%%)",100*bt$estimate,100*ci[1],100*ci[2])
  pcis <- c(pcis,pci)
}

tb$pci=c(pcis,"")
tb <- full_join(pTrend,tb) %>% 
  arrange(desc(n)) %>% 
  slice(-1)

tb %>% select(`Histological type`=Tumor,pci,pTrend) %>% 
  rename(`Proportion (CI95%)`=pci) %>% 
  rename(`Trend test p-value`=pTrend) %>% 
  slice(-nrow(tb)) %>% 
  flextable() %>% 
  autofit() %>% 
  set_caption("Relative proportions of the different kind of malignant CMT")


```


```{r}
bm <- mammary_tumours %>% 
  filter(Benign_Malignant=="Malignant") %>% 
  tabyl(ResultDate,`Primary_tumour`) %>% 
  adorn_percentages() %>% 
  select(ResultDate,`Carcinoma arising in mixed benign tumour`,`Sarcoma NOS`)

mdlC <- lm(`Carcinoma arising in mixed benign tumour`~ResultDate, data=bm)
mdlS <- lm(`Sarcoma NOS`~ ResultDate,data=bm)
bciC <- 100*cbind(coef(mdlC),confint(mdlC))[2,]
bciS <- 100*cbind(coef(mdlS),confint(mdlS))[2,]
tasaInc <- c(sprintf("%.2f%% (%.2f%%, %.2f%%)", bciC[1], bciC[2],bciC[3]),
          sprintf("%.2f%% (%.2f%%, %.2f%%)", bciS[1], bciS[2],bciS[3]))
tibble(Tumor=c("Carcinoma arising in mixed benign tumour","Sarcoma NOS"),
       incRate=tasaInc) %>% 
  flextable() %>% 
  autofit() %>% 
  set_caption("Tasa de incremento anual de la proporción de all_tumours")

```

## Table 3: Relative proportions of the different kind of benign CMT. 

```{r}
bm <- mammary_tumours %>% 
  filter(Benign_Malignant=="Benign") %>% 
  tabyl(ResultDate,`Primary_tumour`) %>% 
  adorn_totals("col")

pt <- NULL
for (k in 2:(ncol(bm)-1)){
  pTrend <- bm %>% select(all_of(k),Total) %>% 
    prop_trend_test() %>% 
    pull(p) %>% 
    formatPval
  pt <- c(pt,pTrend)
}

pTrend <- tibble(Tumor=names(bm)[2:(ncol(bm)-1)],pTrend=pt) 

tb <- mammary_tumours %>% 
  filter(Benign_Malignant=="Benign") %>% 
  tabyl(`Primary_tumour`) %>% 
  adorn_totals("row") %>% 
  rename(Tumor=`Primary_tumour`)

pcis <- NULL
for (k in 1:(nrow(tb)-1)){
  bt <- binom.test(tb$n[k],tb$n[nrow(tb)])
  ci <- bt$conf.int
  pci <- sprintf("%.2f%% (%.2f%%; %.2f%%)",100*bt$estimate,100*ci[1],100*ci[2])
  pcis <- c(pcis,pci)
}

tb$pci=c(pcis,"")
tb <- full_join(pTrend,tb)

tb %>% select(Tumor,pci,pTrend) %>% 
  rename(`p (CI95%)`=pci) %>% 
  slice(-nrow(tb)) %>% 
  flextable() %>% 
  autofit() %>% 
  set_caption("Relative proportions of the different kind of benign CMT")

bm <- mammary_tumours %>% 
  filter(Benign_Malignant=="Benign") %>% 
  tabyl(ResultDate,`Primary_tumour`) %>% 
  adorn_totals("col")

pt <- NULL
for (k in 2:(ncol(bm)-1)){
  pTrend <- bm %>% select(all_of(k),Total) %>% 
    prop_trend_test() %>% 
    pull(p) %>% 
    formatPval
  pt <- c(pt,pTrend)
}

pTrend <- tibble(Tumor=names(bm)[2:(ncol(bm)-1)],pTrend=pt) 

tb <- mammary_tumours %>% 
  filter(Benign_Malignant=="Benign") %>% 
  tabyl(`Primary_tumour`) %>% 
  adorn_totals("row") %>% 
  rename(Tumor=`Primary_tumour`)

pcis <- NULL
for (k in 1:(nrow(tb)-1)){
  bt <- binom.test(tb$n[k],tb$n[nrow(tb)])
  ci <- bt$conf.int
  pci <- sprintf("%.2f%% (%.2f%%; %.2f%%)",100*bt$estimate,100*ci[1],100*ci[2])
  pcis <- c(pcis,pci)
}

tb$pci=c(pcis,"")
tb <- full_join(pTrend,tb) %>% 
arrange(desc(n)) %>% 
slice(-1)

tb %>% select(`Histological type`=Tumor,pci,pTrend) %>% 
  rename(`Proportion (CI95%)`=pci) %>% 
  rename(`Trend test p-value`=pTrend) %>% 
  slice(-nrow(tb)) %>% 
  flextable() %>% 
  autofit() %>% 
  set_caption("Relative proportions of the different kind of malignant CMT")

## El único que cambia con el tiempo es el intraductal papillary, pero el cambio es que los únicos 9 casos que había se registraron entre 2003 y 2007 y después no se ha vuelto a detectar ninguno.

```

## Table 4: Number of cases and controls by the year of birth.

```{r}
# Selection of cases and controls.
# ==============================
controls <- zoocan %>% 
  filter(Year_of_birth>=2003 & Year_of_birth<=2013) %>% 
  select(Breed, Year_of_birth,Island) 
controls <- controls %>% 
  mutate(id=paste("CTRL",1:nrow(controls),sep="-"),
         group="Control")

cases <- mammary_tumours %>% 
  filter(Year_of_birth>=2003&Year_of_birth<=2013)%>% 
  mutate(id=as.character(ReportRef)) %>% 
  select(id,Breed,Year_of_birth,Island) %>% 
  unique() %>% 
  mutate(group="Case")

study <- full_join(controls,cases) %>% 
  mutate(group=factor(group, levels=c("Control","Case")),
         Breed=ifelse(Breed=="No_info",NA,Breed))

study %>% 
  group_by(Year_of_birth,group) %>% 
  summarize(n=n()) %>% 
  pivot_wider(everything(),names_from = group,values_from = n) %>%
 adorn_totals(where="row") %>% 
  flextable()

```

##Fig. 1: Breed distribution in cases and controls.

```{r fig.width=8,fig.height=9}
breedFreqs <- study %>% 
  mutate(Breed=ifelse(Breed=="Other","Other breeds",Breed),
         group=factor(group,levels=c("Case","Control"),labels=c("Cases","Controls"))) %>% 
  tabyl(Breed,group,show_na = FALSE) %>% 
  adorn_percentages("col") %>% 
  pivot_longer(cols =-1, names_to = "Group", values_to = "pct") %>% 
  mutate(pct=100*pct) 

brFreqs <- breedFreqs %>% filter(pct>=1.5) %>% 
  group_by(Breed) %>% 
  summarise(pct=sum(pct)) %>% 
  pull(Breed)

breedFreqsT <- breedFreqs %>% 
  mutate(Breed=ifelse(Breed %in% brFreqs,Breed,"Other breeds")) %>% 
  group_by(Breed,Group) %>% 
  summarize(pct=sum(pct,na.rm=TRUE)) %>% 
  ungroup()

brLevels <- breedFreqsT %>% 
  group_by(Breed) %>% 
  summarise(pct=sum(pct)) %>% 
  arrange(pct) %>% 
  pull(Breed)

breedFreqsT %>% 
  mutate(Breed=factor(Breed,levels=brLevels)) %>% 
  ggplot(aes(x=Breed, y=pct, fill=Group)) +
  geom_col(position = "dodge") +
  geom_text(aes(label=sprintf("%.2f%%",pct)), 
            position=position_dodge(width=0.9),
            hjust=-0.1, col="black", size=3.5) +
  coord_flip() +
  ylim(0,45) +
  ylab("Percentage in group") +
  theme_bw() +
  scale_color_aaas() +
  ##ggtitle("Fig. 1: Breed distribution in cases and controls.")
  theme(panel.border = element_blank(), axis.line = element_line())
ggsave("Fig1.tiff",dpi=300,width=6,height=8)
```

```{r}
breedFreqsT %>% 
  mutate(pct=sprintf("%.2f%%",pct)) %>% 
  pivot_wider(everything(),names_from = Group, values_from = pct) %>% 
  mutate(Breed=factor(Breed,levels=brLevels)) %>% 
  arrange(desc(Breed)) %>% 
  flextable %>% 
  autofit()

```


## Fig. 2v3: Relative proportion of CMT and neutering rate evolution in FD. 

```{r}
pctNeuter <- all_tumoursFD %>% 
  tabyl(ResultDate,Neuter_status) %>% 
  adorn_percentages() %>% 
  rename(Year=ResultDate)

pctTum <- pctTum %>% rename(Year=ResultDate,Prop_CMT=`%Tumours_CMT`)%>% 
  mutate(Prop_CMT=Prop_CMT/100)%>% 
  select(Year,Prop_CMT)

full_join(pctNeuter, pctTum)%>% 
  select(`Year of diagnosis`=Year,Prop_CMT, Neutered)%>% 
  pivot_longer(-`Year of diagnosis`,names_to = "Variable",values_to = "Pct") %>% 
  mutate(Variable=factor(Variable, levels=c("Prop_CMT","Neutered"),
                         labels=c("Relative proportion of CMT","Neutered rate of FD suffering from any tumour"))) %>%
  ggplot(aes(x=`Year of diagnosis`,y=100*Pct,color=Variable)) +
  geom_point()+
  geom_line()+
  geom_smooth(alpha=0.3) +
  labs(color="")+
  theme_bw() +
  scale_color_aaas() +
  theme(panel.border = element_blank(), axis.line = element_line()) +
  ylab("%") +
  theme(legend.position="bottom") +
  guides(color=guide_legend(nrow=2,byrow=TRUE))  
  ##ggtitle("Fig. 2: Relative proportion of CMT and neutering rate evolution in FD")
  ggsave("Fig2.tiff",dpi=300,width=7,height=6)


```



## Evolution of the relative proportion of CMT in FD.
```{r}
ntumoursH <- all_tumoursFD %>% 
  mutate(Location=ifelse(Location=="Mammary gland","Mammary gland","Other Location"),
         Location=ifelse(is.na(Location),"Other Location",Location))

nMama <- ntumoursH %>% 
  tabyl(ResultDate,Location) %>% 
  adorn_totals("col")%>% 
  rename(Year=ResultDate)

pctMama <- nMama %>% 
  select(-Total) %>% 
  adorn_percentages() 

pTrend <- nMama %>% select(`Mammary gland`,Total) %>% 
  prop_trend_test() %>% 
  pull(p) %>% 
  formatPval

tbPctMama <- pctMama %>% 
  select(Year,`Mammary gland`) %>% 
  adorn_pct_formatting(2) %>% 
  rename(pct_Mamm_Tum=`Mammary gland`)

comparaProp(nMama[,c(1,2,4)]) %>% 
  flextable() %>% 
  autofit()
```

# Correlation between CMT proportions and neutering rate in tumours diagnosed from 2003 to 2020

```{r}
all_tumoursFD <- all_tumoursFD %>% 
  mutate(location=ifelse(Location=="Mammary gland","Mammary gland","Other location"),
         location=ifelse(is.na(location),"Other location",location))

tt <- all_tumoursFD %>% select(ReportRef,ResultDate,Neuter_status,location) %>% unique()
ttm <- tt %>% filter(location=="Mammary gland") %>% rename(location1=location)
tto <- tt %>% filter(location!="Mammary gland") %>% rename(location2=location)
tta <- full_join(ttm,tto)
tta <- tta %>% 
  mutate(location1=ifelse(is.na(location1),"",location1),
         location2=ifelse(is.na(location2),"",location2),
         location=paste(location1,location2),
         mammaryTumour=ifelse(location1=="Mammary gland","CMT","non-CMT"))

# Hembras que tienen al menos un tumor de mama.
cmt <- tta %>% filter(mammaryTumour=="CMT") %>% select(ResultDate,Neuter_status)

pTrendNEUT <- tta  %>% tabyl(ResultDate,Neuter_status)%>% 
  adorn_totals("col")%>% 
  select(c(Neutered,Total))%>% 
  prop_trend_test() %>% 
  pull(p) %>% 
  formatPval

pTrendCMT <- tta %>% tabyl(ResultDate,mammaryTumour) %>% 
  adorn_totals("col") %>% 
  select(c(CMT,Total)) %>% 
  prop_trend_test() %>% 
  pull(p) %>% 
  formatPval

nNeut1 <- cmt  %>% tabyl(ResultDate,Neuter_status) %>% 
  adorn_totals("col")

pctNeut1 <- tta  %>% tabyl(ResultDate,Neuter_status) %>% 
  adorn_percentages() 

neutCMT <- pctNeut1 %>% 
  rename(Year=ResultDate) %>% 
  full_join(pctMama) %>% 
  select(Neutered,`Mammary gland`) 

cor.test(neutCMT$Neutered,neutCMT$`Mammary gland`)
```



## Fig. 3: Relative proportion of malignant CMT diagnosed by year.

```{r}

bm <- mammary_tumours %>% tabyl(ResultDate,Benign_Malignant) %>% 
  adorn_totals("col")

pTrend <- bm %>% select(Malignant,Total) %>% 
  prop_trend_test() %>% 
  pull(p) %>% 
  formatPval

bm %>% 
  select(-Total) %>% 
  adorn_percentages()%>% 
  ggplot(aes(x=ResultDate,y=Malignant))+
  geom_point() +
  geom_line() +
  geom_smooth(method="lm",lty=2, alpha=0.3) +
  annotate("text",x=2015,y=0.82,label=paste("Cochran-Armitage trend test p =",pTrend)) + 
  ##ggtitle("Fig. 3 Relative proportion of malignant CMT diagnosed by year.")+
  scale_color_aaas() +
  theme_bw() +
  theme(panel.border = element_blank(), axis.line = element_line()) +
  labs(y="Proportion of malignant CMT.", x="Year of diagnosis") 
ggsave("Fig3.tiff",dpi=300,width=7,height=5)

tb <- mammary_tumours %>% tabyl(Benign_Malignant,show_na=FALSE) %>% 
  adorn_totals("row")


```


## Proporción de tumores de mama malignos

```{r}
bt <- binom.test(tb$n[2],tb$n[3])
ci <- bt$conf.int
pci <- sprintf("%.2f%% (%.2f%%; %.2f%%)",100*bt$estimate,100*ci[1],100*ci[2])
tibble(`p (CI95%)`=pci) %>% 
  flextable() %>% 
  autofit() %>% 
  set_caption("Proporción de tumores de mama malignos")
```



## Fig. 4: Proportion of FD with multiple CMT.

```{r}
nMultiple <- all_tumoursFD %>% 
  rename(Year=ResultDate) %>% 
  filter(Location=="Mammary gland") %>% 
  group_by(Year,ReportRef) %>% 
  summarize(n=n()) %>% 
  mutate(type=ifelse(n==1,"Single","Multiple")) %>% 
  tabyl(Year,type) %>% 
  adorn_totals("col")

pctMultiple <- nMultiple %>% 
  select(-Total) %>% 
  adorn_percentages() 

pTrend <- nMultiple %>% select(Multiple,Total) %>% 
  prop_trend_test() %>% 
  pull(p) %>% 
  formatPval

pctMultiple %>% 
  adorn_percentages("row") %>% 
  ggplot(aes(x=Year,y=100*Multiple)) +
  geom_point()+
  geom_line()+
  geom_smooth(alpha=0.3) +
  #geom_smooth(method="lm",lty=2,se=FALSE) +
  annotate("text",x=2016,y=0.11,label=paste("Cochran-Armitage trend test p =",pTrend))+
  theme_bw() +
  scale_color_aaas() +
  theme(panel.border = element_blank(), axis.line = element_line()) +
  labs(y="Percentage of FD suffering from greater than 1 CMT at diagnosis",x="Year of diagnosis") 
  ##ggtitle("Fig. 4: Proportion of FD with multiple CMT.")+
ggsave("Fig4.tiff",dpi=300,width=7,height=5)

```


## Fig. 5: Age distribution according to the type of CMT.

```{r}
mammary_tumours %>% 
  group_by("Histological type"=Etiqueta3_morfologia) %>% 
  summarise(n=n(), `Mean ± sd`=meanSd(Age)) %>% 
  flextable() %>% 
  merge_v(j=1) %>% 
  autofit() %>% 
  hline(i=3) %>% 
  fix_border_issues()

```


```{r}
mammary_tumours %>% 
  ggplot(aes(x=Etiqueta3_morfologia,y=Age, fill=Etiqueta3_morfologia)) +
  geom_boxplot() +  xlab("Histological type of tumour") + guides(fill="none")+
  #ggtitle("Age distribution by type of tumour")+
  scale_fill_aaas(alpha=0.8) +
  theme_bw() +
  theme(panel.border = element_blank(), axis.line = element_line()) +
  labs(y="Age") 
  ##ggtitle("Fig. 5: Age distribution according to the type of CMT.") +
ggsave("Fig5.tiff",dpi=300,width=7, height=5)

```

## Table 5: Ages of diagnoses of the different CMT in both neutered and entire FD.
 
```{r}
mammary_tumours %>% 
  group_by("Neuter status"=Neuter_status,"Histological type of tumour"=Etiqueta3_morfologia) %>% 
  summarise(n=n(), `Mean ± sd`=meanSd(Age)) %>% 
  flextable() %>% 
  merge_v(j=1) %>% 
  autofit() %>% 
  hline(i=3) %>% 
  fix_border_issues()
```

## Table 6: Differences on ages at which the different CMT are diagnosed.

```{r}
library(broom)
mammary_tumours %>% 
  group_by(Etiqueta3_morfologia) %>% 
  summarize(tidy(t.test(Age~Neuter_status,conf.level=0.95^(1/3)))) %>% 
  mutate(`CI95%`=sprintf("(%.2f,%.2f)",conf.low,conf.high),
         p=formatPval(p.adjust(p.value)))%>% 
  select("Histological type"=Etiqueta3_morfologia,Difference=estimate, `CI95%`,
         `Entire mean age`=estimate1,`Neutered mean age`=estimate2,`P-value`=p) %>% 
  flextable() %>% 
  autofit()

```


```{r}
mdl <- lm(Age~Etiqueta3_morfologia*Neuter_status,data=mammary_tumours)
anova(mdl)

mdl <- aov(Age~Etiqueta3_morfologia+Neuter_status,data=mammary_tumours)
tukTum <- TukeyHSD(mdl,"Etiqueta3_morfologia")
tt1 <- rownames_to_column(data.frame(tukTum$Etiqueta3_morfologia),"Comparison")
tukNeut <- TukeyHSD(mdl,"Neuter_status")
tt2 <- rownames_to_column(data.frame(tukNeut$Neuter_status),"Comparison")
rbind(tt1,tt2) %>% 
  mutate(`Difference (CI95%)`=sprintf("%.2f (%.2f, %.2f)",diff,lwr,upr),
         p=formatPval(p.adj)) %>% 
  select(Comparison, `Difference (CI95%)`,"P-value"=p) %>% 
  flextable() %>% 
  autofit()
```


## Table 7: OR among different breeds when compared with crossbreed dogs, adjusted by island.

```{r}
study <- study %>% 
 mutate(Island=relevel(factor(Island),ref="Gran Canaria"))
# Filtro las razas con al menos 5 casos
ncc <- study %>% 
  tabyl(Breed,group) %>% 
  filter(Case>4 &!is.na(Breed)) %>% 
  select(Breed) %>% 
  pull() %>% as.character()

study <- study %>% 
  filter(Breed%in%ncc) %>% 
  mutate(Breed=relevel(factor(Breed),ref="Crossbreed"))
```

```{r}
# Modelo logístico Incluyendo la isla
modelo <- glm(group~Island+Breed,data=study,family="binomial")
```

```{r}
# Tabla de Odds-Ratios, CI, y número de casos y controles
# =======================================================

# Casos y controles por raza
tbr <- study %>% tabyl(Breed,group,show_na=FALSE) %>% 
  adorn_percentages("col") %>% 
  adorn_pct_formatting(affix_sign=FALSE) %>% 
  adorn_ns(position = "front") %>% 
  rename(`Controls No. (%)`=Control,
         `Cases No. (%)`=Case, 
         value=Breed) %>% 
  mutate(effect="Breed") %>% 
  select(effect,everything())

# Casos y controles por isla
tbi <- study %>% tabyl(Island,group) %>% 
  adorn_percentages("col") %>% 
  adorn_pct_formatting(affix_sign=FALSE) %>% 
  adorn_ns(position = "front") %>% 
  rename(`Controls No. (%)`=Control,
         `Cases No. (%)`=Case, 
         value=Island) %>% 
  mutate(effect="Island") %>% 
  select(effect,everything())

# Casos y controles
tb <- bind_rows(tbi,tbr)
tbref <- tb %>% 
  filter(value=="Gran Canaria"|value=="Crossbreed") %>% 
  mutate(OR=1)

# Función para ordenar los datos de la tabla en orden descendente de OR
ordena <- function(data) arrange(data,desc(OR))

# Construcción de la tabla
ORCI <- exp(confint.default(modelo))
sm <- summary(modelo)$coefficients
effect <- ifelse(grepl("Island",rownames(sm)),"Island",
                    ifelse(grepl("Breed",rownames(sm)),"Breed","--"))
value <- gsub("Breed","",rownames(sm))
value <- gsub("Island","",value)
rownames(sm) <- rownames(ORCI) <- NULL
ORtable <- cbind(data.frame(effect=effect,value=value),sm,ORCI) %>% 
  rename(p=`Pr(>|z|)`) %>% 
  slice(-1) %>% 
  #filter(effect=="Island"| (effect=="Breed"&p<0.05)) %>% 
  mutate("P-value"=formatPval(p),
         OR=round(exp(Estimate),2),
         CI=sprintf("[%.2f,%.2f]", `2.5 %`,`97.5 %`)) %>% 
  select(effect,value,OR,CI,"P-value") %>% 
  group_by(effect) %>% 
  nest() %>% 
  mutate(data=map(data,ordena)) %>% 
  unnest(cols=data) %>% 
  left_join(tb) %>% 
  full_join(tbref) 

ORtable[c(44,1:6,45,7:43),] %>% 
  flextable() %>% 
  autofit() %>% 
  hline(i=7) %>% 
  merge_v(j=1) %>% 
  valign(j=1,valign="top")
```

```{r}
# Forestplot con los OR significativos`
#signif <- which(summary(modelo)$coefficients[,4]<0.05)
OR=exp(coef(modelo)) #[signif]
ci.OR=exp(confint.default(modelo)) #[signif,]
tbOR=data.frame(coef=OR[-1], lower=ci.OR[-1,1], upper=ci.OR[-1,2])
tbOR <- tbOR[grep("Breed",rownames(tbOR)),]
rownames(tbOR)=gsub("Breed","",rownames(tbOR))
o=order(tbOR$coef,decreasing=TRUE)
tbOR=tbOR[o,]
ref <- data.frame(coef=1,lower=NA,upper=NA)
rownames(ref) <- "Crossbreed"
##tbOR <- rbind(ref,tbOR)
tiff("Fig6.tiff",width=7,height=10,units = "in", res=300)
fp <- forestplot(rownames(tbOR), 
           tbOR,
           fn.ci_norm = fpDrawCircleCI,
           boxsize=0.3,
           col=fpColors(box="royalblue",line="darkblue", 
                        summary="royalblue", zero="grey60"),
           clip=c(-3,15),
           lty.ci = 7,
           zero=1,
           txt_gp = fpTxtGp(ticks=gpar(cex=0.9)),
           xlab="OR")

dev.off()



```