---
title: "Analyses for AJBA Revision, Salmi et al."
author: "Dan Hall, Hanna Kim, UGA Statistical Consulting Center"
date: "`r format(Sys.Date(),'%B %d, %Y')`"
output:
  html_document:
    toc: TRUE
#    df_print: paged
    code_folding: hide     
  pdf_document: 
  word_document: 
    toc: TRUE
editor_options: 
  chunk_output_type: console
---


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

This file contains results from the first phase of field work conducted over the years 2004-2005. These results are based on a corrected data set sent by Roberta on July 11, 2025.



```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}

library(readxl)
library(here)
library(lme4); library(lmerTest)
library(car)
library(emmeans)
library(tidyverse)
library(ordinal)
library(nnet)
library(glmmTMB)
library(visreg)
library(ggeffects)
library(MASS)
library(mclogit)
library(ggmosaic)
library(mathjaxr) 
library(gridExtra)
library(aod)
library(vcdExtra)
library(knitr)
conflicted::conflict_prefer("select","dplyr")
conflicted::conflict_prefer("filter","dplyr")

options(contrasts=c("contr.sum","contr.sum"))

# All data
initData <- read_xlsx(here('Corrected_FinalDataset_2025.xlsx'),
                     sheet = 3,skip=1,
                     col_names=all.vars(~focal+focalNo+sex+date+month+year+monthYear+season+infant+infantAge+minsObserved+minsAbove+minsGround+minsHt0+minsHt1+minsHt2+minsHt3+minsFeedingAbove+minsRestingAbove+minsMovingAbove+minsSocializingAbove+minsFruitAbove+minsLeavesAbove+minsOtherFoodAbove),
                     col_types=c("text","numeric","numeric","date","text","numeric","text","text","numeric","text",rep("numeric",14)))
# infantAge sometimes has two ages (e.g., "1, 5") so treating it as text

# Creating factors and other variables 
allData <- initData %>% 
  mutate(
    seasonFac=factor(season,
                     levels=c("LOW","MID","HIGH"),
                     labels=c("Low Fr.","Mid Fr.","High Fr.")),
    sexFac=factor(sex, levels=c(0,1), labels = c("Male","Female")),
    focalFac=factor(focal),
    infantFac=factor(infant,levels = c(0, 1), 
                     labels = c("Absent", "Present")),
    propMinAbove=minsAbove/minsObserved,
    dateFac=factor(date),
    obsFac=factor(1:nrow(initData)))

# Create phase 1/2 data
phase1Data <- allData %>% 
  filter(year<2009) %>% 
  mutate(focalFac=fct_drop(focalFac))
phase2Data <- allData %>% 
  filter(year>=2009) %>% 
  mutate(focalFac=fct_drop(focalFac))

```


# 1. Two-way Table 

This is a 2-way table to examine the data distribution across each factors. We can see that the data are well-distributed across each factor.

```{r, echo=FALSE} 
with(phase1Data,ftable(seasonFac,focalFac,row.vars=1))
```


# 2. Arboreality 

## 2.1 Model Selection 


For Arboreality, we recommended an analysis based on a **logistic mixed-effect model**. In this model, the response is similar to a binomial proportion, and the explanatory variables include both fixed and random effects. Fixed effects would include those for season and focal (instead of sex) and possible interaction among these factors. We considered random day effects because, on some days, responses are taken on multiple individuals and such data are likely to be similar (correlated). The inclusion of random observation effects is a solution to deal with the over-dispersion problem which is common in binomial response models. 

Note. Overdispersion arises when the assumptions of the binomial distribution are not satisfied. In particular, the binomial distribution applies to responses that can be thought of the number of successes out of n trials, where the success probabilities are constant and the outcomes of the trials are mutually independent. Although number of minutes engaged in an activity out of total number of minutes observed is superficially similar to a binomial outcome, the minute-by-minute results are certainly not independent. In such situations, a model built on the binomial distribution can be used if additional variability is incorporated into the model to account for the failure of the distributional assumption. In this case, that variability is built in through observation-specific random effects. 

```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
Q1.final <- glmmTMB(cbind(minsGround,minsAbove)~1+
                 seasonFac*focalFac+
                 (1|dateFac)+(1|obsFac),
               contrasts=list(seasonFac=contr.sum,
                              focalFac=contr.sum),
               data=phase1Data,
               family=binomial())
```


Therefore, the below equation is the final model.

$$
log(\frac{\pi_{ijl}}{1-\pi_{ijl}}) = \mu + \alpha_i + \beta_j + \alpha\beta_{ij}  + b_l + c_{ijl},
$$


$$\alpha_i : \text{season effect ,   i = Low, Mid, High} $$
$$\beta_j : \text{focal effect ,   j = Beatrice, Ebuka, Kingo, Mama, Mekome} $$

$$b_l : \text{random day effect, ~} N(0,\sigma^2_{b}) $$
$$c_{ijl} : \text{random obs effect, ~} N(0,\sigma^2_{c}) $$


## 2.2 Results

The below ANOVA table represents that "main effect of focal" is significant.
Now, we present more details about the marginally significant effect. 

```{r, echo=FALSE}
Anova(Q1.final,type=3)
```

### 2.2.1 Effect of Focal (+Sex)      

We give estimated marginal means for each animal (each level of focalFac) (these are marginal estimates of the probability of being above ground) and estimate pairwise odds ratio. In addition, odds ratio comparing Kingo to the females, on average is presented. This odds ratio quantifies the sex effect.

We find that:

- The proportion of time spent above ground is lowest for Kingo (7.4%), while Beatrice has the highest proportion (17.2%).
- All odds ratios are not significant except for Beatrice/Kingo (p-value: 0.0001) and Kingo/Mama (p-value: 0.0022). In other words, The C.I for Beatrice/Kingo and Kingo/Mama do not include 1. 
- The Beatrice/Kingo odds ratio of 2.603 means that the odds of spending time above ground are estimated to be 2.603 times higher for Beatrice than for Kingo.



```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
# emmeans: focal 
Q1.final.emm <- emmeans(Q1.final,specs=list(focal=~focalFac),type="response")
```



```{r, echo=FALSE}
###############################################
# Figure 3 ####################################
###############################################

# emmeans: focal

#plot(Q1.final.emm$focal, CIs=TRUE,type="response", xlab="Estimated marginal probability") 
Q1.final.emm$focal. <- as.data.frame(Q1.final.emm$focal)
ggplot(Q1.final.emm$focal., aes(x = focalFac, y = prob, color = as.factor(focalFac))) +
  geom_point(size=2) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2) + 
  labs(x = "Focal Animal", y = "Estimated marginal probability", color = "Focal") + 
  theme_minimal() + 
  theme(legend.position = "none") -> Fig3
Fig3 
ggsave(here("Fig3.eps"),plot=Fig3,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig3.tiff"),plot=Fig3,device="tiff",units="mm",dpi=800,scale=1)
```


```{r, echo=FALSE}
# emmenas: focal-pairwise (O.R)
a<-contrast(Q1.final.emm$focal, method="pairwise",type="response")%>% summary(infer=c(T,T))
ggplot(a, aes(x = contrast, y = odds.ratio)) +
  geom_point() +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1) +
  labs(x = "Contrast", y = "Odds Ratio") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7)) +
  geom_hline(yintercept = 1, color = "red", linetype = "dashed") 



Q1.final.emm$focal
a
```


In addition, one of our interests is the comparison by sex, and the results are shown below.

We find that there is a significant difference in the probability of time spent above ground between gender. (odds of arboreality are estimated to be 2.03 times higher for females, on average, than for Kingo) 


```{r, echo=FALSE}
# emmenas: focal-sex (O.R)
a<-contrast(Q1.final.emm$focal,
            method=list(female.v.male=c(1,1,-4,1,1)/4), type="response") %>% summary(infer=c(T,T)); a

```


# 3. Height 


## 3.1 Model Selection--Generalized logit Model 

For Height, we originally suggested an analysis based on a cumulative logistic mixed-effect model. However, based on the mosaic plots below, we can see that only seasonFac satisfies proportional odds assumption. This can be seen from the fact that season has a similar effect on the probability of being at a low height as it does on the probability of being at a low or middle height. In contrast, focal has a different effect on being at a low height than it does no being at a low or middle height. Also, although the response variable is divided into three levels of height, which looks like an ordinal variable, the reasons why gorillas choose certain heights may be due to factors other than height itself. For example, gorillas may prefer the lowest and highest locations because more food is available there. Based on both empirical and theoretical reasons (according to Dr. Salmi), it does not appear that the heights are ordered from low to middle to high in terms of preference or food type, and proportional odds (i.e., parallel slopes) assumptions fail to hold. Moreover, we initially tried using proportional odds models and the hypothesis of proportional odds failed to hold then formally tested based on these models. For these reasons, we instead based our analysis on generalized logit models, a class of models typically used for nominal response variables. As with arboreality, we consider a mixed-effect version of this sort of model, to account for heterogeneity from day to day and overdispersion.

Now, we will proceed with a **generalized logistic mixed-effect model** rather than the previously proposed cumulative model.


```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
long1_data_ht <- phase1Data %>%
  pivot_longer(cols = 15:17, names_to = "Category", values_to = "Count") 


long1_data_ht %>% as.data.frame() %>% head()

long1_data_ht2 <- long1_data_ht %>%
  mutate(Category = factor(Category, 
                           levels = c("minsHt1", "minsHt2", "minsHt3"), labels=c("Low Ht.", "Mid Ht.", "High Ht."),
                           ordered = TRUE))
```



```{r, warning=FALSE}
ggplot(long1_data_ht2) +
  geom_mosaic(aes(x=product(seasonFac),fill=Category,weight=Count))
ggplot(long1_data_ht2) +
  geom_mosaic(aes(x=product(focalFac),fill=Category,weight=Count))
```


In this model, the response would be where the gorilla is in the tree (low, medium, high) given that it is above ground, and the explanatory variables include both fixed and random effects. Fixed effects would include those for season and sex (instead of focal animal) and possible interaction among these factors. We considered random day and random observation effects like previous model. Especially, Ebuka did not spend time at high height across all seasons. The fact that the explanatory factors in this study are not fully crossed places limitation on the models that can be entertained. Thus, sex (instead of focal animal) is included in the model. 

The below equation is the final model.


$$ 
log(\frac{\pi_{ijl,mid}}{\pi_{ijl,low}}) = \mu_{mid} + \alpha_{i,mid} + \beta_{j,mid}   + \alpha\beta_{ij,mid}  + b_{l,mid} + c_{ijl,med}
$$ 
$$ 
log(\frac{\pi_{ijl,high}}{1-\pi_{ijl,low}}) = \mu_{high} + \alpha_{i,high} + \beta_{j,high}  + \alpha\beta_{ij,high}  +  b_{l,high} + c_{ijl,high}
$$ 


$$\alpha_i : \text{season effect ,   i = Low, Mid, High} $$
$$\beta_j : \text{sex effect ,   j = Male, Female} $$


$$b_l : \text{random day effect, ~} N(0,\sigma^2_{b}) $$
$$c_{ijl} : \text{random obs effect, ~} N(0,\sigma^2_{c}) $$




## 3.2 Results

```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
fit1.1 <- mblogit(Category~1+seasonFac*sexFac
                  ,
                  weights=Count,
                  random=list(~1|dateFac, ~1|obsFac),
                  data=long1_data_ht2,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=FALSE,
                  link=logit) #@ update for adding random effect v3
summary(fit1.1)


coef(fit1.1)
length(coef(fit1.1))


# 11/12 new
lM_season <- cbind(matrix(0,4,2),diag(1,4,4),matrix(0,4,6))
wald.test(Sigma=vcov(fit1.1), b=coef(fit1.1), L=lM_season) %>%
  print(digits=4) # seas sig 
lM_sex <- cbind(matrix(0,2,6),diag(1,2,2),matrix(0,2,4))
wald.test(Sigma=vcov(fit1.1), b=coef(fit1.1), L=lM_sex) %>%
  print(digits=4) # sex NOT sig 
lM_sex_seas <- cbind(matrix(0,4,8),diag(1,4,4))
wald.test(Sigma=vcov(fit1.1), b=coef(fit1.1), L=lM_sex_seas) %>%
  print(digits=4) # sex:seas  sig 


```



The below table represents that the interaction between seas:sex is significant. Now, we present more details about the significant interaction. 



```{r, echo=FALSE}
data <- data.frame(
  Term = c("Season", "Sex",  "season:sex"),
  D.f. = c("4", "2", "4"),
  Chi.sq = c("63.81","10.31", "11.13"),
  `P-value` = c("4.58e-13 ***", "0.0058 **", "0.02518 *")
)

kable(data, col.names = c("Term", "D.F", "Chi.sq", "P-value"))
```


```{r, echo=FALSE}
# emmeans
fit1.1.emm <- emmeans(fit1.1,specs=list(seas=~seasonFac,
                                        sex=~sexFac,
                                        seasXsex=~seasonFac*sexFac), mode="prob", by="Category")


```


### 3.2.3 Effect of Season:Sex (by Sex)


We give estimated marginal means for the probability of being at each height for each combination of season by sex and estimate a pairwise difference between season for each sex.


We find that:

- Season effects on the probability distribution across the three heights depend on Sex (and vice versa). 

- Specifically, For male, the probability of being at low height is low during LOW and MID seasons, but increases sharply in the HIGH season. Meanwhile, the probability of being in the mid height gradually decreases as the season progresses. The probability of being in the high height is high during the LOW and MID seasons, but drops sharply in the HIGH season.

- For females, the probability of being in the low height gradually increases with the season, while the probability of being in the mid/high height steadily decreases.

- For male, the difference in the probability of being in the low height is significant between the LOW vs HIGH seasons and between the MID vs HIGH seasons. Similarly, the difference in the probability of being in the high height is significant between the LOW vs HIGH season and the MID vs HIGH seasons.

- For female, the difference in the probability of being in the low height is significant between the LOW vs MID seasons and LOW vs HIGH seasons. Similarly, the difference in the probability of being in the mid height is significant between the LOW vs HIGH season. In addition, the difference in the probability of being in the high height is significant between the LOW vs MID season and LOW vs HIGH season.




```{r, echo=FALSE}
# emmeans: seas:sex by sex

a <- as.data.frame(fit1.1.emm$seasXsex)
ggplot(a, aes(x = seasonFac, y = prob, color = sexFac, group = sexFac)) +
  geom_point(position = position_dodge(width = 0.3),size=2) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.5,position = position_dodge(width = 0.3)) + 
  facet_wrap(~ Category) + 
  labs(x = "Season",  y = "Estimated marginal probability", color="Sex") + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7))

```

```{r, echo=FALSE}
# emmeans: seas:sex by sex - pairwise Diff.
a<-contrast(fit1.1.emm$seasXsex,method="pairwise",by=c("sexFac","Category"),type="response") %>% summary(infer=c(T,T))
ggplot(a, aes(x = contrast, y = estimate, color = as.factor(sexFac), group = sexFac)) +
  geom_point(position = position_dodge(width = 0.2),size=2) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.5,position = position_dodge(width = 0.2)) +  
  facet_wrap(~ Category) + 
  labs(x = "Season", y = "Estimated contrast in marginal probabilities", color = "Sex") +
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7))+
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") 

fit1.1.emm$seasXsex
a
```



Based on the discussion during the meeting on November 26, the revised interaction plot for season:sex is shown below. Each panel represents a season (LOW, MID, HIGH), with the x-axis corresponding to Height (Low, Mid, High) and the plots differentiated by Sex.

```{r, echo=FALSE}
###############################################
# Figure 4 ####################################
###############################################
# @@1126 
a <- as.data.frame(fit1.1.emm$seasXsex)
ggplot(a, aes(x = factor(Category,
                         levels=c("Low Ht.","Mid Ht.","High Ht."),
                         labels=c("Low","Mid","High")), 
              y = prob, color = sexFac, group = sexFac)) +
  geom_point(position = position_dodge(width = 0.3),size=2) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.5,position = position_dodge(width = 0.3)) + 
  facet_wrap(~ seasonFac) + 
  labs(x = "Tree Height",  y = "Estimated marginal probability", color="Sex") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7)) -> Fig4
Fig4
ggsave(here("Fig4Alt.eps"),plot=Fig4,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig4Alt.tiff"),plot=Fig4,device="tiff",units="mm",dpi=800,scale=1)

```




### 3.2.4 Effect of Season:Sex (by Season)

We give estimated marginal means of the probability of being at each height for each combination of sex by season and estimate a pairwise differences in these probabilities between sexes for each season.



We find that:

- For all seasons, the probability of a female being in the low height is higher than that of a male. Similarly, the probability of a female being in the mid height is higher than that of a male (except for HIGH season). However, the probability of being in the high height is higher for males.

- In the MID season, the probabilities of being in low/high height differ significantly between sex.

- In the LOW and HIGH seasons, the probabilities of being in each height differ very little between sex.


```{r, echo=FALSE}
# emmeans: seas:sex by seas

a <- as.data.frame(fit1.1.emm$seasXsex)
ggplot(a, aes(x = sexFac, y = prob, color = seasonFac, group = seasonFac)) +
  geom_point(position = position_dodge(width = 0.3),size=2) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.5,position = position_dodge(width = 0.3)) + 
  facet_wrap(~ Category) + 
  labs(x = "Sex",  y = "Estimated marginal probability", color="Season") + 
  theme_minimal()

```


```{r, echo=FALSE}
# emmenas: seas:sex by seas - pairwise Diff.
a<-contrast(fit1.1.emm$seasXsex,method="pairwise",by=c("seasonFac","Category"),type="response") %>% summary(infer=c(T,T))
ggplot(a, aes(x = contrast, y = estimate, color = seasonFac, group = seasonFac)) +
  geom_point(position = position_dodge(width = 0.3),size=2) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.5,position = position_dodge(width = 0.3)) +  
  facet_wrap(~ Category) + 
  labs(x = "Sex", y = "Estimated contrast in marginal probabilities", color = "Season") +
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7))+
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") 

fit1.1.emm$seasXsex
a
```

## 3.3 Model Selection--Cumulative logit Model 

I double-checked whether a cumulative logit model would suffice here now that the data have changed. The answer still seems to be no.

```{r, echo=FALSE, results="hide", message=FALSE, warning=FALSE,eval=FALSE}
m1.Ht.clogit <- clmm2(Category~1+
                      seasonFac*sexFac,  
                    random=dateFac , # since response is not binomial, no need to overdispersion problem 
                    weights=Count,
                    link="logistic", 
                    control = clmm2.control(method = "nlminb"), Hess = TRUE, 
                    data=long1_data_ht2)
summary(m1.Ht.clogit)
AIC(m1.Ht.clogit) # cum logit
AIC(fit1.1) # generalized logit
nominal_test(m1.Ht.clogit)

m1a.Ht.clogit <- clmm(Category~1+
                      seasonFac*sexFac+(1|dateFac),  
                    weights=Count,
                    link="logit", 
                    control = clmm.control(method = "nlminb"), Hess = TRUE, 
                    data=long1_data_ht2)
AIC(m1.Ht.clogit,m1a.Ht.clogit) # same
nominal_test(m1a.Ht.clogit)

m2.Ht.clogit <- clmm2(Category~1+
                      seasonFac*sexFac,  
                    random=dateFac , # since response is not binomial, no need to overdispersion problem 
                    weights=Count,
                    nominal=~seasonFac*sexFac,
                    link="logistic", 
                    control = clmm2.control(method = "nlminb",maxIter=100), Hess = TRUE, 
                    data=long1_data_ht2)


m3.Ht.clogit <- clm(Category~1+
                      seasonFac*focalFac,  
                    weights=Count,
                    link="logit", 
                    data=long1_data_ht2)
nominal_test(m3.Ht.clogit) # seems that prop odds assumption does not hold
```

# 4. Activity 

## 4.1 Model Selection

For Activity, we suggested an analysis based on a **generalized logistic mixed-effect model**. In this model, the nominal response would be which type of activity (feeding, resting, moving, social) and the explanatory variables include both fixed and random effects. Fixed effects would include those for season, focal and possible interaction among these factors. We considered random day and random observation effects like previous model. Here, we plan to collapse the response category levels from (feeding, resting, moving, social) to (feeding, resting (including social), moving).

The below equation is the final model.


$$ 
log(\frac{\pi_{ijl,move}}{\pi_{ijl,feed}}) = \mu_{move} + \alpha_{i,move} + \beta_{j,move} +\alpha\beta_{ij,move}  +  b_{l,move} + c_{ijl,move}
$$ 
$$ 
log(\frac{\pi_{ijl,rest}}{\pi_{ijl,feed}}) = \mu_{rest} + \alpha_{i,rest} + \beta_{j,rest}  +\alpha\beta_{ij,rest}  + b_{l,rest} + c_{ijl,rest}
$$ 





$$\alpha_i : \text{season effect ,   i = Low, Mid, High} $$
$$\beta_j : \text{focal effect ,   j = Beatrice, Ebuka, Kingo, Mama, Mekome}$$

$$b_l : \text{random day effect, ~} N(0,\sigma^2_{b}) $$
$$c_{ijl} : \text{random obs effect, ~} N(0,\sigma^2_{c}) $$

## 4.2 Results

The below table represents that focal effect is significant.
Now, we present more details about the marginally significant effect. 

```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
long1_data_act <- phase1Data %>%
  pivot_longer(cols = 18:21, names_to = "Category", values_to = "Count") 

long1_data_act %>% as.data.frame() %>% head()

long1_data_act2 <- long1_data_act %>%
  mutate(Category = factor(Category, 
                           levels = c("minsFeedingAbove", "minsRestingAbove", "minsMovingAbove", "minsSocializingAbove"), labels=c("Feeding", "Resting", "Moving", "Socializing"), 
                           ordered = TRUE))

###################### @@ 10/1 OPTION 1. collapse response 

# OPTION 1. collapse response 
long1_data_act3 <- long1_data_act %>%
  mutate(long1_data_act2,newcat=fct_collapse(Category,Feeding="Feeding", Moving="Moving", other_level="Resting/Socializing"))


xtabs(Count~newcat+seasonFac+focalFac,data=long1_data_act3)
xtabs(Count~newcat+seasonFac+sexFac,data=long1_data_act3)


fit2.2 <- mblogit(newcat~1+
                      seasonFac*focalFac 
                   , 
                  weights=Count,
                  random=list(~1|dateFac, ~1|obsFac),
                  data=long1_data_act3,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=FALSE) #@ update for adding random effect v3
summary(fit2.2)

###########@10/26 ANOVA table RE-TRY!!!! ###########
coef(fit2.2)
length(coef(fit2.2))


# 11/12
lM_season <- cbind(matrix(0,4,2),diag(1,4,4),matrix(0,4,24))
wald.test(Sigma=vcov(fit2.2), b=coef(fit2.2), L=lM_season) %>%
  print(digits=4) # seas not sig 
lM_focal <- cbind(matrix(0,8,6),diag(1,8,8), matrix(0,8,16))
wald.test(Sigma=vcov(fit2.2), b=coef(fit2.2), L=lM_focal) %>%
  print(digits=4) # focal sig
lM_seas_focal <- cbind(matrix(0,16,14),diag(1,16,16))
wald.test(Sigma=vcov(fit2.2), b=coef(fit2.2), L=lM_seas_focal) %>%
  print(digits=4) # seas:foc not  sig


# emmeans
fit2.2.emm <- emmeans(fit2.2,specs=list(focal=~focalFac), mode="prob", by="newcat")

```

```{r, echo=FALSE}
data <- data.frame(
  Term = c("Season", "Focal", "season:focal"),
  D.f. = c("4", "8", "16"),
  Chi.sq = c("13.56","42.02", "15.20"),
  `P-value` = c("0.0088", "1.34e-06 *", "0.5101")
)

kable(data, col.names = c("Term", "D.F", "Chi.sq", "P-value"))
```


### 4.2.1 Effect of Focal (+Sex)

We give estimated marginal means for Focal and estimate a pairwise difference of proportion.

- The probability of feeding is estimated to be higher in Kingo than others while, for both moving and other activities, the probability of that activity is estimated to be lower in Kingo than females.

- Some differences in the estimated probabilities of feeding/moving/resting between Focal are significant.

```{r, echo=FALSE}
###############################################
# Figure 6 ####################################
###############################################
# emmeans: focal

#plot(fit1.1.emm$focal, CIs=TRUE,type="response", xlab="Probability") 
fit2.2.emm$focal. <- as.data.frame(fit2.2.emm$focal)
ggplot(fit2.2.emm$focal., aes(x = focalFac, y = prob, color = as.factor(focalFac))) +
  geom_point(size=1.8) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.3) + 
  labs(x = "Focal animal", y = "Estimated marginal probability", color = "Focal") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7),
        legend.position = "none") +
  facet_wrap(~ newcat) -> Fig6
Fig6
ggsave(here("Fig6Alt.eps"),plot=Fig6,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig6Alt.tiff"),plot=Fig6,device="tiff",units="mm",dpi=800,scale=1)


# emmeans: focal-pairwise (diff)
a<-contrast(fit2.2.emm$focal, method="pairwise",type="response")%>% summary(infer=c(T,T));a
ggplot(a, aes(x = contrast, y = estimate)) +
  geom_point() +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1) +
  labs(x = "Contrast", y = "Estimated contrast in marginal probabilities") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7)) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  facet_wrap(~  newcat) 


fit2.2.emm$sex
a
```



In addition, one of our interests is the comparison by sex, and the results are shown below.

We find that there is a significant difference in the probability of feeding/moving/resting between gender.  


```{r, echo=FALSE}
# emmenas: focal-sex (prop diff)
a<-contrast(fit2.2.emm$focal,
            method=list(female.v.male=c(1,1,-4,1,1)/4), type="response") %>% summary(infer=c(T,T)); a # @ update v3

```

### 4.2.2 Effect of Season

We give estimated marginal probabilities and contrasts between those probabilities for season.

```{r, echo=FALSE}
###############################################
# Figure 7 ####################################
###############################################
# emmeans: season

fit2.2.seas <- emmeans(fit2.2,specs=list(season=~seasonFac), 
                      mode="prob", by="newcat")

fit2.2.seas %>% as.data.frame() %>% 
  mutate(seasonFac=factor(seasonFac,
                          levels=c("Low Fr.","Mid Fr.","High Fr."))) %>%
  ggplot(aes(x = seasonFac, y = prob, color = seasonFac)) +
  geom_point(size=2) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.3) + 
  labs(x = "Season", y = "Estimated marginal probability", color = "Season") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 0, hjust = 1, size=7),
        legend.position="none") +
  facet_wrap(~ newcat) -> Fig7
Fig7
ggsave(here("Fig7Alt.eps"),plot=Fig7,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig7Alt.tiff"),plot=Fig7,device="tiff",units="mm",dpi=800,scale=1)

# emmeans: season-pairwise (diff)
a <- contrast(fit2.2.seas$season, method="pairwise",
              by="newcat", type="response") %>% summary(infer=c(T,T))
ggplot(a, aes(x = contrast, y = estimate)) +
  geom_point() +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1) +
  labs(x = "Contrast", y = "Estimated contrast in marginal probabilities") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7)) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  facet_wrap(~  newcat) 


fit2.2.seas$season
a
```


- The probability of feeding is estimated to be significantly higher in low season than in mid season and significantly higher in low season than in high season.

- The probability of moving is estimated to be significantly lower in low season than in mid season.


# 5. Food 

## 5.1 Model Selection 

For Food, we suggested an analysis based on a **generalized logistic mixed-effect model**. In this model, the nominal response would be which type of food (leaves/bark, fruit/flower, herb, other) and the explanatory variables include both fixed and random effects. Fixed effects would include those for season, focal and possible interaction among these factors. We considered random day and random obs effects like previous model. Here, we plan to collapse the response category levels from (leaves/bark, fruit/flower, herb, other) to (leaves/bark, fruit/flower, other (including herb)).

The below equation is the final model.


$$ 
log(\frac{\pi_{ijl,leav}}{\pi_{ijl,fruit}}) = \mu_{leav} + \alpha_{i,leav} + \beta_{j,leav}   + \alpha\beta_{ij,leav} + b_{l,leav} + c_{ijl,leav}
$$ 
$$ 
log(\frac{\pi_{ijl,other}}{\pi_{ijl,fruit}}) = \mu_{other} + \alpha_{i,other} + \beta_{j,other}  + \alpha\beta_{ij,other} + b_{l,other} + c_{ijl,other}
$$ 


$$\alpha_i : \text{season effect ,   i = Low, Mid, High} $$

$$\beta_j : \text{sex effect ,    j = 1,2 for males, females}$$




$$b_l : \text{random day effect, ~} N(0,\sigma^2_{b}) $$
$$c_{ijl} : \text{random obs effect, ~} N(0,\sigma^2_{c}) $$

* Note that a model of similar form with focal effects rather than sex effects did not converge, so a simpler model that replaced focal by sex was used.

## 5.2 Results



```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
long1_data_food <- phase1Data %>%
  pivot_longer(cols = 22:24, names_to = "Category", values_to = "Count") 

long1_data_food %>% as.data.frame() %>% head()

long1_data_food2 <- long1_data_food %>%
  mutate(Category = factor(Category, 
                           levels = c("minsFruitAbove", "minsLeavesAbove", "minsOtherFoodAbove"), labels=c("Fruit", "Leaves", "Other"),
                           ordered = TRUE))


###################### @@ 10/1 OPTION 1. collapse response 

# OPTION 1. collapse response 
long1_data_food3 <- long1_data_food2 %>%
  mutate(newcat=Category)

xtabs(Count~newcat+seasonFac+focalFac,data=long1_data_food3)
xtabs(Count~newcat+seasonFac+sexFac,data=long1_data_food3)


# fit3.2 <- mblogit(newcat~1+
#                     seasonFac*focalFac,
#                   weights=Count,
#                   random=list(~1|dateFac, ~1|obsFac),
#                   data=long1_data_food3,
#                   control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
#                   dispersion=TRUE) #@ update for adding random effect v3
# Did not converge

fit3.2a <- mblogit(newcat~1+
                    seasonFac*sexFac, 
                  weights=Count,
                  random=list(~1|dateFac, ~1|obsFac),
                  data=long1_data_food3,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=TRUE) #@ update for adding random effect v3
summary(fit3.2a)


###########@10/26 ANOVA table RE-TRY!!!! ###########
coef(fit3.2a)
length(coef(fit3.2a))

## @ 1112
lM_season <- cbind(matrix(0,4,2),diag(1,4,4),matrix(0,4,6))
wald.test(Sigma=vcov(fit3.2a), b=coef(fit3.2a), L=lM_season) %>%
  print(digits=4)# seas sig 
lM_sex <- cbind(matrix(0,2,6),diag(1,2,2), matrix(0,2,4))
wald.test(Sigma=vcov(fit3.2a), b=coef(fit3.2a), L=lM_sex) %>%
  print(digits=4) # sex nonsig
lM_seas_sex <- cbind(matrix(0,4,8),diag(1,4,4))
wald.test(Sigma=vcov(fit3.2a), b=coef(fit3.2a), L=lM_seas_sex) %>%
  print(digits=4) # sex nonsig # seas:sex not  sig


# emmeans
fit3.2a.emm <- emmeans(fit3.2a,specs=list(seas=~seasonFac, seasXsex=~seasonFac*sexFac), mode="prob", by="newcat")

```

The table below shows that the main effect of season and season:focal effect are significant. Now, we present more details about the marginally significant effects. 


```{r, echo=FALSE}
data <- data.frame(
  Term = c("season", "sex", "season:sex"),
  D.f. = c("4", "8", "16"),
  Chi.sq = c("137.9","0.0100", "7.39"),
  `P-value` = c("0 ***", "0.995", "0.117")
)

kable(data, col.names = c("Term", "D.F", "Chi.sq", "P-value"))
```



### 5.2.1 Effect of Season

We give estimated marginal mean probabilities of eating each food type for each season and estimate pairwise differences in these probabilities across the three seasons. 

- The probability of eating Fruit increases across Season, while the probability of eating Leaves decreases across Season.

- All pairwise differences between seasons in the probabilities of eating Fruit are significant. The same is true for Leaves.

```{r, echo=FALSE}
###############################################
# Figure 8 ####################################
###############################################
# emmenas: seas

#plot(fit1.1.emm$seas, CIs=TRUE,type="response", xlab="Probability") 
fit3.2a.emm$seas. <- as.data.frame(fit3.2a.emm$seas)
ggplot(fit3.2a.emm$seas., aes(x = seasonFac, y = prob, color = as.factor(seasonFac))) +
  geom_point(size=2) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2) + 
  labs(x = "Season", y = "Estimated marginal probability", color = "Season") + 
  theme_minimal()+
  theme(legend.position="none") +
  theme(axis.text.x = element_text(angle = 0, hjust = 1, size=7)) +
  facet_wrap(~ newcat) -> Fig8
Fig8
ggsave(here("Fig8Alt.eps"),plot=Fig8,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig8Alt.tiff"),plot=Fig8,device="tiff",units="mm",dpi=800,scale=1)

# emmeans: seas-pairwise (diff)
a<-contrast(fit3.2a.emm$seas, method="pairwise",type="response")%>% summary(infer=c(T,T))
ggplot(a, aes(x = contrast, y = estimate)) +
  geom_point() +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1) +
  labs(x = "Contrast", y = "Estimated contrast in marginal probabilities") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7)) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  facet_wrap(~  newcat) 



fit3.2a.emm$seas
a
```

### 5.2.2 Effect of Season:Sex

The test of the interaction between season and sex reported above is a joint test of the interaction across both response functions. That is, it is a test of whether these factors interact conducted with respect to the joint distribution of the 3-category response. Since this test is not significant, in the interest of avoiding an inflated Type I error rate, I advise not conducting further tests like those below that examine the interaction effects of these factors on each response function of the model and that examine sex differences separately by season. 

However, it is possible to look more specifically at a test of interaction in each response function. Such tests are conducted below. The first one tests whether season and sex interact in their effects on the log odds of eating leaves rather than fruit, and the second one tests whether season and sex interact on the odds of eating other food rather than fruit. Both are non-significant.  

```{r, echo=FALSE}
lM_seas_sex1 <- cbind(matrix(0,2,8),diag(1,2,2),diag(0,2,2))
wald.test(Sigma=vcov(fit3.2a), b=coef(fit3.2a), L=lM_seas_sex1) %>%
  print(digits=4) # sex nonsig # seas:sex not  sig
lM_seas_sex2 <- cbind(matrix(0,2,8),diag(0,2,2),diag(1,2,2))
wald.test(Sigma=vcov(fit3.2a), b=coef(fit3.2a), L=lM_seas_sex2) %>%
  print(digits=4) # sex nonsig # seas:sex not  sig

```

It is also possible to test for equal probabilities of eating each food type across sex, separately in each season. Those tests are conducted below.

```{r, echo=FALSE}
contrast(fit3.2a.emm$seasXsex,method="pairwise",type="response",by=c("newcat","seasonFac"))

```

The season by sex interaction plot is given below.

```{r}
fit3.2a.emm$seasXsex %>% as.data.frame() %>%
  ggplot(aes(x = sexFac, y = prob, color = seasonFac, group = seasonFac)) +
  geom_point(position = position_dodge(width = 0.3)) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1,position = position_dodge(width = 0.3)) + 
  facet_wrap(~ newcat) + 
  labs(x = "Sex",  y = "Estimated marginal probability", color="Season") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7))

```


# 6. Arboreality by time

Below, conditional density plots are constructed showing the distribution of arboreality over time separately by sex.

```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
###############################################
# Figure 2 ####################################
###############################################

timeData <- read_xlsx(here('Corrected_FinalDataset_2025.xlsx'),sheet = 1)
names(timeData)

timeData$Time %>% hour() %>% head(30)
timeData$Time %>% minute() %>% head(30)

timeData2 <- mutate(timeData,
                  dayTime = hour(Time)+minute(Time)/60,
                  sexFac=factor(Focal=="KI",levels=c(TRUE,FALSE),
                                labels=c("Male","Female")),
                  yFac=factor(`On/Above`,levels=0:1,
                              labels=c("On ground","In trees")))
timeData2%>%as.data.frame()%>%head()

#cdplot(yFac~dayTime,data=timeData2 , xlab="Time", ylab="Probability")

plotData <- timeData2%>%select(dayTime,yFac,sexFac) %>% na.omit()
## Plot
ggplot(plotData, aes(dayTime, after_stat(count), fill = yFac)) + 
  geom_density(position='fill', alpha = .9) + 
  xlab("Time (hrs)") + labs(y="Estimated marginal probability",fill='') +
  theme(legend.text=element_text(size=12), 
        axis.title=element_text(size=14)) +
  facet_wrap(~sexFac,nrow=2) -> Fig2
Fig2
ggsave(here("Fig2.eps"),plot=Fig2,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig2.tiff"),plot=Fig2,device="tiff",units="mm",dpi=800,scale=1)

```


# 7. Comparison of Activity Distribution On Ground vs In Trees

```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
groundData <- read_xlsx(here('Corrected_FinalDataset_2025.xlsx'),
                     sheet = 2,skip=1,
                     col_names=all.vars(~focal+focalNo+sex+date+month+year+monthYear+season+infant+infantAge+minsObserved+minsGround+minsAbove+minsHt0+minsHt1+minsHt2+minsHt3+minsFeeding_Above+minsResting_Above+minsMoving_Above+minsSocializing_Above+minsFruit_Above+minsLeaves_Above+minsOtherFood_Above+minsFeeding_Ground+minsResting_Ground+minsMoving_Ground+minsSocializing_Ground+minsFruit_Ground+minsLeaves_Ground+minsOtherFood_Ground),
                     col_types=c("text","numeric","numeric","date","text","numeric","text","text","numeric","text",rep("numeric",21))) %>% 
  mutate(
    seasonFac=factor(season,levels=c("LOW","MID","HIGH"),
                     labels=c("Low Fr.","Mid Fr.","High Fr.")),
    sexFac=factor(sex, levels=c(0,1), labels = c("Male","Female")),
    focalFac=factor(focal),
    infantFac=factor(infant,levels = c(0, 1), 
                     labels = c("Absent", "Present")),
    minsOther_Above=minsResting_Above+minsSocializing_Above,
    minsOther_Ground=minsResting_Ground+minsSocializing_Ground,
    propMinAbove=minsAbove/minsObserved,
    dateFac=factor(date),
    obsFac=factor(1:length(focal)))


```

## 7.1 Model 

To compare the activity distribution between on the ground versus in the trees, we use a **generalized logistic mixed-effect model** similar to that of section 4, but with random effects for each focal by day combination. This accomplishes the pairing of the comparison between activity distribution on and off the ground for the same animal on the same day. In this model, the nominal response would be activity type (coarsened to 3 categories: feeding, other=resting/socializing, moving) and the explanatory variables include both fixed and random effects. Fixed effects are included for a factor indicating whether the measurements are taken on the ground or in the trees. In addition to random effects for each focal by day combination we also include observation-level random effects to account for overdispersion. 

The model equation is shown below. Here, $i,j,k$ are indices for focal animal, day, and ground vs. tree, respectively.   

$$ 
log(\frac{\pi_{ijk,move}}{\pi_{ijk,feed}}) = \mu_{move} + \alpha_{k,move} + b_{ij,move} + c_{ijk,move}
$$ 
$$ 
log(\frac{\pi_{ijk,other}}{\pi_{ijk,feed}}) = \mu_{other} + \alpha_{k,other} + b_{ij,other} + c_{ijk,other}
$$ 


$$\alpha_i : \mbox{ground vs tree effect ,   $i = 1,2$ for ground,tree, respectively} $$

$$b_{ij,move},b_{ij,other} : \text{random animal by day effects, jointly} \sim N(0,D_{b}) $$
$$c_{ijl,move},c_{ijl,other} : \text{random observation effects, jointly} \sim N(0,D_{c}) $$

## 7.2 Results

```{r, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
names(groundData)
groundData %>% select(c(1:9,18,20,36,25,27,37)) %>% names() 

long1 <- groundData %>% select(c(1:9,18,20,36,25,27,37)) %>%  
  pivot_longer(cols = 10:15, names_to = c("activity","location"), 
               values_to = "count",
               names_sep="_") %>%
  mutate(locationFac=factor(location,
                            levels=c("Ground","Above"),
                            labels=c("On ground","In trees")),
         dateFac=factor(date),
         focalFac=factor(focal),
         obsFac=factor(1:length(focal)),
         focalDateFac=factor(focalFac:dateFac),
         activityFac=factor(activity,
                            levels=c("minsFeeding","minsMoving","minsOther"),
                            labels=c("Feeding","Moving","Other"))) %>%
  arrange(date,focal,activity,location)
long1 %>% as.data.frame() %>% head()

long2 <- rbind(
  groundData %>% select(c(1:9,18,20,36)) %>%
    rename(minsFeeding=minsFeeding_Above,
           minsOther=minsOther_Above,
           minsMoving=minsMoving_Above) %>%
    mutate(ground=0),
    groundData %>% select(c(1:9,25,27,37)) %>%
    rename(minsFeeding=minsFeeding_Ground,
           minsOther=minsOther_Ground,
           minsMoving=minsMoving_Ground) %>%
    mutate(ground=1)
) %>%
  mutate(groundFac=factor(ground,levels=1:0,
                          labels=c("On Ground","In Trees")),
         dateFac=factor(date),
         focalFac=factor(focal),
         focalDateFac=factor(focalFac:dateFac),
         obsFac=factor(1:length(focal))) %>%
  arrange(date,focal,ground)
long2 %>% as.data.frame() %>% head()

xtabs(count~locationFac+activityFac,data=long1)
cmhTab <- xtabs(count~locationFac+activityFac+dateFac+focalFac,data=long1)

#CMHtest(cmhTab,overall=TRUE)$ALL


m1.ground <- mblogit(cbind(minsFeeding,minsMoving,minsOther)~1+groundFac, 
                  random=list(~1|focalDateFac, ~1|obsFac),
                  data=long2,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=FALSE) #@ update for adding random effect v3
summary(m1.ground)

m0.ground <- mblogit(cbind(minsFeeding,minsMoving,minsOther)~1, 
                  random=list(~1|focalDateFac, ~1|obsFac),
                  data=long2,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=FALSE) #@ update for adding random effect v3

coef(m1.ground)

# 2-df test of groundFac:
# wald.test(Sigma=vcov(m1.ground), 
#           b=coef(m1.ground), 
#           L=cbind(matrix(0,2,2),diag(1,2,2))) %>%
#   print(digits=5) # groundFac is significant 

# emmeans
m1.ground.emm <- emmeans(m1.ground,specs=list(ground=~groundFac),
                         mult.names="category", 
                         mode="prob", by="category")

# m1.ground.emm <- emmeans(m1.ground,specs=list(ground=~groundFac),
#                          mult.names="category", 
#                          mode="prob", by="category")

```

A test of the overall effect of location (on ground vs above ground) on the distribution of activity is given below.

```{r}
wald.test(Sigma=vcov(m1.ground), 
          b=coef(m1.ground), 
          L=cbind(matrix(0,2,2),diag(1,2,2))) %>%
  print(digits=5) # groundFac is significant 
```

Estimates of the probability of each activity type by location are given below. Note that "Category" refers to activity type (1=Feeding, 2=Moving, 3=Resting or Socializing).

```{r}
m1.ground.emm$ground

```

Tests of the effects of location on the probability of each activity type are given below.

```{r}
contrast(m1.ground.emm$ground, method="pairwise",by="category", type="response")%>% summary(infer=c(T,T))
```

A plot of the estimated probabilities of each activity by location is given below.

```{r, echo=FALSE}
###############################################
# Figure 5 ####################################
###############################################
as.data.frame(m1.ground.emm$ground) %>% mutate(catFac=factor(category,levels=1:3,labels=c("Feeding","Moving","Resting/Socializing"))) %>%
  ggplot(aes(x = groundFac, y = prob, color = groundFac)) +
  geom_point(size=1.5) +  
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.15) + 
  labs(x = "Location", y = "Estimated marginal probability", color = "Location") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=0)) +
  ylim(0,1) +
  facet_wrap(~ catFac) -> Fig5
Fig5
ggsave(here("Fig5Alt.eps"),plot=Fig5,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig5Alt.tiff"),plot=Fig5,device="tiff",units="mm",dpi=800,scale=1)
```



# 8. Comparison of Food Choice On Ground vs In Trees

```{r, echo=FALSE, results="hide",, message=FALSE, warning=FALSE}
groundData <- read_xlsx(here('Corrected_FinalDataset_2025.xlsx'),
                     sheet = 2,skip=1,
                     col_names=all.vars(~focal+focalNo+sex+date+month+year+monthYear+season+infant+infantAge+minsObserved+minsGround+minsAbove+minsHt0+minsHt1+minsHt2+minsHt3+minsFeeding_Above+minsResting_Above+minsMoving_Above+minsSocializing_Above+minsFruit_Above+minsLeaves_Above+minsOtherFood_Above+minsFeeding_Ground+minsResting_Ground+minsMoving_Ground+minsSocializing_Ground+minsFruit_Ground+minsLeaves_Ground+minsOtherFood_Ground),
                     col_types=c("text","numeric","numeric","date","text","numeric","text","text","numeric","text",rep("numeric",21))) %>%
  mutate(
    seasonFac=factor(season,levels=c("LOW","MID","HIGH"),
                     labels=c("Low Fr.","Mid Fr.","High Fr.")),
    sexFac=factor(sex, levels=c(0,1), labels = c("Male","Female")),
    focalFac=factor(focal),
    infantFac=factor(infant,levels = c(0, 1),
                     labels = c("Absent", "Present")),
    minsOther_Above=minsResting_Above+minsSocializing_Above,
    minsOther_Ground=minsResting_Ground+minsSocializing_Ground,
    propMinAbove=minsAbove/minsObserved,
    dateFac=factor(date),
    obsFac=factor(1:length(focal)))


# minsFeeding_Above
# minsFruit_Above, minsLeaves_Above, minsOtherFood_Above
# minsFeeding_Ground
# minsFruit_Ground, minsLeaves_Ground, minsOtherFood_Ground
```

## 8.1 Model

To compare the food type distribution when on the ground versus when in the trees, we use a **generalized logistic mixed-effect model** similar to that of section 4, but with random effects for each focal by day combination. This accomplishes the pairing of the comparison between activity distribution on and off the ground for the same animal on the same day. In this model, the nominal response would be food type (with 3 categories: fruit, leaves, other food) and the explanatory variables include both fixed and random effects. Fixed effects are included for a factor indicating whether the measurements are taken on the ground or in the trees. In addition to random effects for each focal by day combination we also include observation-level random effects to account for overdispersion.

The model equation is shown below. Here, $i,j,k$ are indices for focal animal, day, and ground vs. tree, respectively.

$$
log(\frac{\pi_{ijk,fruit}}{\pi_{ijk,other}}) = \mu_{fruit} + \alpha_{k,fruit} + b_{ij,fruit} + c_{ijk,fruit}
$$
$$
log(\frac{\pi_{ijk,leaves}}{\pi_{ijk,other}}) = \mu_{leaves} + \alpha_{k,leaves} + b_{ij,leaves} + c_{ijk,leaves}
$$


$$\alpha_i : \mbox{ground vs tree effect ,   $i = 1,2$ for ground,tree, respectively} $$

$$b_{ij,fruit},b_{ij,leaves} : \text{random animal by day effects, jointly} \sim N(0,D_{b}) $$
$$c_{ijl,fruit},c_{ijl,leaves} : \text{random observation effects, jointly} \sim N(0,D_{c}) $$

## 8.2 Results

```{r, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
# minsFeeding_Above
# minsFruit_Above, minsLeaves_Above, minsOtherFood_Above
# minsFeeding_Ground
# minsFruit_Ground, minsLeaves_Ground, minsOtherFood_Ground

names(groundData)
groundData %>% select(c(1:9,22:24,29:31)) %>% names()


long3 <- rbind(
  groundData %>% select(c(1:9,22:24)) %>%
    rename(minsFruit=minsFruit_Above,
           minsLeaves=minsLeaves_Above,
           minsOtherFood=minsOtherFood_Above) %>%
    mutate(ground=0),
    groundData %>% select(c(1:9,29:31)) %>%
    rename(minsFruit=minsFruit_Ground,
           minsLeaves=minsLeaves_Ground,
           minsOtherFood=minsOtherFood_Ground) %>%
    mutate(ground=1)
) %>%
  mutate(groundFac=factor(ground,levels=1:0,
                          labels=c("On Ground","In Trees")),
         dateFac=factor(date),
         focalFac=factor(focal),
         focalDateFac=factor(focalFac:dateFac),
         obsFac=factor(1:length(focal))) %>%
  arrange(date,focal,ground)
long3 %>% as.data.frame() %>% head()


m11.ground <- mblogit(cbind(minsFruit,minsLeaves,minsOtherFood)~1+groundFac,
                  random=list(~1|focalDateFac, ~1|obsFac),
                  data=long3,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=FALSE) #@ update for adding random effect v3
summary(m11.ground)

m10.ground <- mblogit(cbind(minsFruit,minsLeaves,minsOtherFood)~1,
                  random=list(~1|focalDateFac, ~1|obsFac),
                  data=long3,
                  control=mmclogit.control(maxit=50),  #@ update for adding random effect v3
                  dispersion=FALSE) #@ update for adding random effect v3

coef(m11.ground)

# emmeans
m11.ground.emm <- emmeans(m11.ground,specs=list(ground=~groundFac),
                         mult.names="category",
                         mode="prob", by="category")


```

A test of the overall effect of location (on ground vs above ground) on the distribution of food choice is given below.

```{r}
wald.test(Sigma=vcov(m11.ground),
          b=coef(m11.ground),
          L=cbind(matrix(0,2,2),diag(1,2,2))) %>%
  print(digits=5) # groundFac is significant
```

Estimates of the probability of each activity type by location are given below. Note that "Category" refers to food type (1=Fruit, 2=Leaves, 3=Other food).

```{r}
m11.ground.emm$ground

```

Tests of the effects of location on the probability of each activity type are given below.

```{r}
contrast(m11.ground.emm$ground, method="pairwise",by="category", type="response")%>% summary(infer=c(T,T))
```

A plot of the estimated probabilities of each activity by location is given below.

```{r, echo=FALSE}
###############################################
# Figure 9 ####################################
###############################################

as.data.frame(m11.ground.emm$ground) %>% mutate(catFac=factor(category,levels=1:3,labels=c("Fruit","Leaves","Other food"))) %>%
  ggplot(aes(x = groundFac, y = prob, color = groundFac)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.15) +
  labs(x = "Location", y = "Estimated marginal probability", color = "Location") +
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=07),legend.position="none") +
  ylim(0,1) +
  facet_wrap(~ catFac) -> Fig9
Fig9
ggsave(here("Fig9.eps"),plot=Fig9,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig9.tiff"),plot=Fig9,device="tiff",units="mm",dpi=800,scale=1)

as.data.frame(m11.ground.emm$ground) %>% mutate(catFac=factor(category,levels=1:3,labels=c("Fruit","Leaves","Other food"))) %>%
  ggplot(aes(x = groundFac, y = prob, color = groundFac)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.15) +
  labs(x = "Location", y = "Estimated marginal probability", color = "Location") +
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=0)) +
  ylim(0,1) +
  facet_wrap(~ catFac) -> Fig9Alt
Fig9Alt
ggsave(here("Fig9Alt.eps"),plot=Fig9Alt,device="eps",units="mm",dpi=800,scale=1)
ggsave(here("Fig9Alt.tiff"),plot=Fig9Alt,device="tiff",units="mm",dpi=800,scale=1)

```