---
title: "Farm characterization for Florence's PhD project "
subtitle: "Farm typology based on multivariate analysis"
author: "Deo"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include = F}
knitr::opts_chunk$set(echo = T)
```
A farm typology is proposed based on a Factor Analysis for Mixed Data (Data) run for the selected variables

## Load the 'clean' dataset, rename and select variables

```{r}
# set working directory ## Change accordingly
setwd('~/Harare 2023/diversification_rotation/Florence')
# load package
require(tidyverse)
rm(list = ls())
clean_survey <- readxl::read_excel('2024-10-07.CLEANED DATA FARMER SEGMENTATION.xlsx',
                                   sheet = 'Selected_data',
                                   skip = 2,
                                   trim_ws = T,
                                   .name_repair = 'universal')
# List of variables 
names(clean_survey)

# Remove unnecessary columns and anonymize respondents 
descrip_survey <- clean_survey |>
  select(!c(Household_Head_Name )) 

famd_df <- descrip_survey |>
  select(!ID.NO) |>
  na.omit()

pca_df <- famd_df |>
  mutate(highestEdu_qualification.HHH. = case_when(
    highestEdu_qualification.HHH. == 'never_attended' ~ 0,
    highestEdu_qualification.HHH. == 'primary_school' ~ 1,
    highestEdu_qualification.HHH. == 'secondary' ~ 2,
    highestEdu_qualification.HHH. == 'tertiary' ~ 3,
    .default = NA) ,
         across(where(is.character),
                function(x) ifelse(x == 'no', 0, ifelse(x == 'yes', 1, NA))))

```


## Matrix of scatterplots to see correlations for pairs of variables
## Then, compare PCA for 3 subsets of the dataset, followed by hierarchical clustering
```{r}

# Correlation matrix to visualize data
P00 <- GGally::ggpairs(famd_df)
png('correlation_matrix.png', height = 25, width = 20, units = 'cm', res = 1000)
P00
ggsave('correlation_matrix.png')
dev.off()

# Mock PCA on reduced data
res_pca <- FactoMineR::PCA(pca_df)
# choose how many dimensions to incorporate in the typology
# the more dimensions, the more variability is captured, and the more noise is introduced
res_pca$eig
# what variables are best represented by each dimension
factoextra::fviz_screeplot(res_pca, addlabels = T)
factoextra::fviz_pca_biplot(res_pca)
factoextra::fviz_pca_biplot(res_pca, axes = c(1, 3))
factoextra::fviz_pca_biplot(res_pca, axes = c(2, 3))
# choose appropriate palette for visualization
pal <- colorRampPalette(c('lightskyblue1', 'darkblue'))
# choose the most relevant variables to describe each dimension
res_pca$var$contrib[, 1:5]
for(i in list(c(1, 2), c(1, 3), c(1, 4), c(1, 5))){
  for(j in names(pca_df)){
    nb_col <- length(unique(pca_df[[j]]))
    print(factoextra::fviz_ellipses(res_pca, 
                                  habillage = j, 
                                  axes = i, 
                                  geom = 'text', 
                                  addEllipses = T, 
                                  palette = pal(nb_col)
                                  )
        )
  }
}
# Hierarchical clustering on the first 5 dimensions that capture 51% of total variability
eg_hc <- res_pca$ind$coord
dist_matrix <- dist(eg_hc)
hc <- hclust(dist_matrix, method = 'ward.D') # could be 'complete' or 'average'
plot(hc, main = 'Subset 1 - Dendogram on the first 5 components', xlab = '', sub = '', cex = 0.6)
# decide on the number of clusters to make: 2, 3, 4 or more
typo_1 <- pca_df
for(cl in 4:2){
  clusters <- cutree(hc, k = cl)
  typo_1$farm_type <- paste0('farm_', clusters)
  print(paste0('================ Number of farm types = ', cl, ' ================'))
  for(i in names(pca_df)){
    res_anova <- lm(typo_1[[i]] ~ typo_1[['farm_type']])
    print(paste0('--------------------', i, '----------------'))
    print(anova(res_anova))
    print('-----------')
    print(emmeans::emmeans(res_anova, 
                           pairwise ~ farm_type, 
                           type = 'response', 
                           adjust = 'bonferroni'))
  }
}
# Choose a few box-plots to display, or any other visual. Here is just an e.g. for per capita income
boxplot_variable <- function(var_name, var_title){
  P00 <- ggplot(typo_1, aes_string('farm_type', var_name)) +
  geom_boxplot() +
  labs(x = 'Farm type', title = var_title) + 
  theme_bw()
  png(paste0('plot_', var_name, '.png'), height = 7.5, width = 10, units = 'cm', res = 1000)
  P00 
  ggsave(paste0('plot_', var_name, '.png')) 
  dev.off()
}
boxplot_variable('HHSize', 'Size of household')
boxplot_variable('HH_income_', 'Household income, USD')
boxplot_variable('PercFarmIncome', 'Proportion of income derived from farming, %')
boxplot_variable('Ngoats', 'Number of goats')
boxplot_variable('Nchickens', 'Number of chickens')
boxplot_variable('TotalLand', 'Farm size, ha ???')
boxplot_variable('AvgMZ_yield.kgs.', 'Maize yield, kg/ha ???')
boxplot_variable('PercentageLegumeArea', 'Share of farm area occupied by legume crops, %')
boxplot_variable('Manure_application_in_Maize_field_during_2022.23_season', 'Manure application')
boxplot_variable('Count_intens_options', 'Count_intens_options')
boxplot_variable('Do_you_pracatice_intercroping.', 'Number of intensification options practised')
boxplot_variable('Do_you_own_any_vegetable_garden.', 'Presence of vegetable garden')
boxplot_variable('Count_sust_practices', 'Count_sust_practices')

```


## Compare multivariate techniques on the same subset of variables
```{r}
##############################################################################
# Mock FADM
res_famd <- FactoMineR::FAMD(famd_df)
factoextra::fviz_screeplot(res_famd, addlabels = T)
# factoextra::fviz_famd(res_famd)

for(i in list(c(1, 2), c(1, 3), c(1, 4) )){
  for(j in names(famd_df)){
    nb_col <- length(unique(famd_df[[j]]))
    print(factoextra::fviz_ellipses(res_famd, 
                                  habillage = j, 
                                  axes = i, 
                                  geom = 'text', 
                                  addEllipses = T, 
                                  palette = pal(nb_col)
                                  )
        )
  }
}
# Hierarchical clustering on the first 5 dimensions that capture 64% of total variability
eg_hc <- res_famd$ind$coord[, 1:5]
dist_matrix <- dist(eg_hc)
hc <- hclust(dist_matrix, method = 'ward.D') # could be 'complete' or 'average'
plot(hc, main = 'Subset 3 - Dendogram on the first 5 components', xlab = '', sub = '', cex = 0.6)
# decide on the number of clusters to make: 2, 3, 4 or more
typo_3 <- famd_df
for(cl in 4:2){
  clusters <- cutree(hc, k = cl)
  typo_3$farm_type <- paste0('farm_', clusters)
  print(paste0('================ Number of farm types = ', cl, '================'))
  for(i in names(famd_df)[c(1, 3:10, 12:15)]){
    res_anova <- lm(typo_3[[i]] ~ typo_3[['farm_type']])
    print(paste0('--------------------', i, '----------------'))
    print(anova(res_anova))
    print('-----------')
    print(emmeans::emmeans(res_anova, 
                           pairwise ~ farm_type, 
                           type = 'response', 
                           adjust = 'bonferroni'))
  }
}


################################################################################
# Mock hierarchical clustering on original selected variables
eg_hc <- pca_df
dist_matrix <- dist(eg_hc)
hc <- hclust(dist_matrix, method = 'ward.D') # could be 'complete' or 'avrage'
plot(hc, main = "Dendrogram for original variables", xlab = "", sub = "", cex = 0.6)

# decide on the number of clusters to make: 2, 3, 4 or more
typo_4 <- pca_df
for(cl in 4:2){
  print(paste0('================ Number of farm types = ', cl, '================'))
  clusters <- cutree(hc, k = cl)
  typo_4$farm_type <- paste0('farm_', clusters)
  for(i in names(eg_hc)){
    res_anova <- lm(typo_4[[i]] ~ typo_4[['farm_type']])
    print(paste0('--------------------', i, '----------------'))
    print(anova(res_anova))
    print('-----------')
    print(emmeans::emmeans(res_anova, 
                           pairwise ~ farm_type, 
                           type = 'response', 
                           adjust = 'bonferroni'))
  }
  
  ggplot(typo_4, aes(hh_size, cultivated_area_ha, colour = farm_type)) + 
  geom_point() +
  theme_bw() +
  labs(title = "Hierarchical Clustering")
  
  ggplot(typo_4, aes(cluster, cultivated_area_ha, colour = farm_type)) + 
  geom_boxplot() +
  theme_bw() +
  labs(title = "Hierarchical Clustering")

}
################################################################################
# Mock k-means clustering
eg_kmeans <- pca_df

# decide on the number of clusters to make: 2, 3, 4 or more
typo_5 <- eg_kmeans
for(cl in 2:4){
  kmeans_result <- kmeans(eg_kmeans, centers = cl)
  clusters <- as.factor(kmeans_result$cluster)
  typo_5$farm_type <- paste0('farm_', clusters)
  
  for(i in names(eg_kmeans)){
    res_anova <- lm(typo_5[[i]] ~ typo_5[['farm_type']])
    print(paste0('--------------------', i, '----------------'))
    print(anova(res_anova))
    print('-----------')
    print(emmeans::emmeans(res_anova, 
                           pairwise ~ farm_type, 
                           type = 'response', 
                           adjust = 'bonferroni'))
  }
  ggplot(typo_5, aes(hh_size, cultivated_area_ha, colour = farm_type)) + 
  geom_point() +
  theme_bw() +
  labs(title = paste0("kmeans clustering with ", cl, ' clusters') )
  
  ggplot(typo_5, aes(kmeans_cluster, cultivated_area_ha, colour = farm_type)) + 
  geom_boxplot() +
  theme_bw() +
  labs(title = paste0("kmeans clustering with ", cl, ' clusters') )
}

```
I hope the cryptic codes were readable; otherwise, feel free to reach out.


I think there were too many variables selected. For 44 respondents, I would not have more than 8 or 10 variables. Consider selecting the most relevant.


I proposed several methods of clustering: but my subjective preference goes to the FAMD or the PCA, followed by a clustering in 2 or 3 groups.
Anyway, whatever method, you would select the contrasting characteristics must be summarized; and I do prefer a series of boxplot for variables showing significant differences. It may be more challenging to give it a coherent explanation.


ALL THE BEST AND SEE YOU SOON!!!