---
title: "Seeded GMM Demo"
author: "C. Durso"
date: "8/5/2024"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(stats)
library(tidyverse)
library(readxl)
library(ggplot2)

library(ggpubr)
library(kdensity)
library(mclust)
source("em_fit_and_visualization_basic.R")

```

## Overview

This analysis uses Gaussian mixture models as a starting point for identifying instars in the width data. The default choice of seed for MClust doesn't always give biologically reasonable mixtures. Using a multiplicative sequence of proposed means offers one lightweight method for pushing the models towards biologically reasonable solutions.

Check 
https://www.webpages.uidaho.edu/~stevel/519/literatures/MCLUST%20for%20R.pdf

https://sites.stat.washington.edu/raftery/Research/PDF/Scrucca2015ADAC.pdf

https://mclust-org.github.io/mclust/reference/hc.html



## Load data

Restrict to columns `Individual` through `Notes`

```{r}
# Generate example from observations
# dat<-read_xlsx("headcapsuledata.xlsx",range="A1:G6320",sheet="Raw Data")
# names(dat)<-c("Individual", "Year.Mom","year","Days.Old","width","diet","County")
# dat<-dat[dat$diet=="BW",]
# set.seed(234567)
# dat<-dat[sample(1:nrow(dat),replace=TRUE),]
# dat<-filter(dat, !is.na(year))
# names(dat)<-str_to_lower(names(dat))
# save(dat,file="data_products/example_data.RData")
load("data_products/example_data.RData")
table(dat$diet)


```

## Plots

Create coarse and fine histograms.

```{r}
ggplot(dat,aes(x=width))+geom_histogram(bins=30)+
  labs(x="head capsule width")
ggplot(dat,aes(x=width))+geom_histogram(binwidth=0.01)+
  labs(x="head capsule width")

```

# Estimate parameters from GMM

As a less ad hoc way of estimating distinct instars than the visual inspection option, try GMM.
So far, this just reports the number of Gaussians used for the maximum likelihood best fit. 


Functions to visualize the GMM fit given the fit and the data 

```{r}
# Density of mixture model
## scaled by the count of observations
f.nls<-function(x,ct,mu,stdev){
  f.vec<-sapply(1:length(ct),function(i){ct[i]*dnorm(x,mu[i],stdev[i])})
  return(rowSums(f.vec))
}

# Individual Gaussians, scaled so area under the curve is in the same units
# as the area of the histogram bars: prop.factor is the number of data items 
# times the binwidth.

f<-function(x,i,fit.this,prop.factor){
  return(fit.this$parameters$pro[i]*prop.factor*dnorm(x,fit.this$parameters$mean[i], sqrt(fit.this$parameters$variance$sigmasq[i])))
}

# histogram, individual Gaussians, and density
gmm.plot.make<-function(fit.this,dat.this,binwidth.this){
  g<-ggplot()+geom_histogram(data=dat.this, aes(x=width), binwidth=binwidth.this,fill="gray")+
    stat_function(fun=f.nls, 
                  args=list(ct=fit.this$parameters$pro*nrow(dat)*binwidth.this,
                            mu=fit.this$parameters$mean,
                            stdev=sqrt(fit.this$parameters$variance$sigmasq)),
                  color="orange")
  
  for(i in 1:fit.this$G){
    g<-g+stat_function(fun=f,args=list(i=i,fit.this=fit.this,
                                       prop.factor=nrow(dat)*binwidth.this))+
      geom_vline(xintercept = fit.this$parameters$mean[i])
  }
  return(g)
}
```


## Apply the Brooks-Dyar seed method


```{r}
fit.mclust<-em.fit.get(dat,7)
gmm.plot.make(fit.mclust,dat,binwidth.this = .01)

```

