# Supporting the codesign of agronomic innovations with data

Frédéric Baudron, Telma Sibanda, Gatien Falconnier

► To cite this version:

Frédéric Baudron, Telma Sibanda, Gatien Falconnier. Supporting the codesign of agronomic innovations with data. CIRAD; CIMMYT, 130 p., 2025. hal-05397083

HAL Id: hal-05397083

https://hal.science/hal-05397083v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Supporting the Codesign of
Agronomic Innovations
with Data

<figure><img src="image_2.png" /></figure>

Frédéric Baudron (CIRAD, Montpellier, France)
Telma Sibanda (CIMMYT, Harare, Zimbabwe)
Gatien Falconnier (CIRAD, Harare, Zimbabwe)

November 26, 2025
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>1.1</td><td>Purpose of this book</td><td>4</td></tr><tr><td>1.2</td><td>KoboToolBox and R</td><td>5</td></tr><tr><td>1.3</td><td>Structure of this book</td><td>6</td></tr><tr><td>1.4</td><td>Mbire District as a case study</td><td>8</td></tr><tr><td>2</td><td>Learning KoboToolBox and R</td><td>11</td></tr><tr><td>2.1</td><td>A quick tutorial on KoboToolbox</td><td>11</td></tr><tr><td>2.2</td><td>Resources to learn R</td><td>15</td></tr><tr><td>3</td><td>Baselining</td><td>16</td></tr><tr><td>3.1</td><td>Data preparation</td><td>16</td></tr><tr><td>3.2</td><td>Preliminary analyses</td><td>29</td></tr><tr><td>3.3</td><td>R Markdown</td><td>49</td></tr><tr><td>4</td><td>Farm typologies</td><td>50</td></tr><tr><td>4.1</td><td>Concept of farm typology</td><td>50</td></tr><tr><td>4.2</td><td>Data preparation</td><td>51</td></tr><tr><td>4.3</td><td>Delineation of farm types</td><td>58</td></tr><tr><td>4.4</td><td>Interpretation of the farm typology</td><td>68</td></tr><tr><td>4.5</td><td>Selection of stratified samples</td><td>75</td></tr><tr><td>5</td><td>Managing on-farm experiments</td><td>77</td></tr><tr><td>5.1</td><td>Generating QR codes</td><td>77</td></tr><tr><td>5.2</td><td>Integrating QR codes in KoboToolBox surveys</td><td>83</td></tr></table>
<table><tr><td>CONTENTS</td><td>3</td></tr><tr><td>6 Evaluating on-farm experiments</td><td>90</td></tr><tr><td>6.1 Visualization of treatments' performance</td><td>91</td></tr><tr><td>6.2 Shiny app</td><td>99</td></tr><tr><td>7 Scaling</td><td>108</td></tr><tr><td>7.1 Data preparation</td><td>108</td></tr><tr><td>7.2 Transferability potential</td><td>117</td></tr><tr><td>8 Conclusions</td><td>124</td></tr><tr><td>Acknowledgments</td><td>125</td></tr></table>
# Chapter 1

# Introduction

## 1.1 Purpose of this book

Living labs - defined as physical or virtual spaces that bring together multiple stakeholders to collaboratively tackle societal challenges, fostering co-creation and real-world testing of solutions (Hossain et al., 2019) - have become increasingly popular in development-oriented agronomy. However, we observe that such codesign processes nurtured in living labs often favor local knowledge over academic knowledge and are rarely informed by systematically collected data (Baudron, 2024).

While local knowledge is invaluable, relying exclusively on it can inadvertently limit the creativity and problem-solving potential of the codesign process (Asten et al., 2009). In practice, this can lead to innovations that reflect existing practices (e.g., manure management, intercropping, minimum tillage) rather than exploring novel solutions to the key bottlenecks in the systems under study.

We argue that researchers have a crucial role in generating, organizing, analyzing, and presenting data at every stage of the codesign process — from diagnostic and needs assessment, to prototyping, testing, and evaluation (Giller et al., 2008). Diagnostics and analytical tools from systems agronomy can help make challenges, opportunities, and trade-offs more explicit, supporting evidence-based decision making and enhancing the creativity of the codesign process.

At the same time, data cannot simply be presented to Living Lab members in raw form. To be useful in codesign, data must be adapted to the audience, for example through visualizations, simplified tables, interactive tools, or role-playing exercises, while still allowing participants to draw their own conclusions. This approach aligns with the principles of discovery learning, where participants engage actively with information, explore patterns, and derive in-
sights themselves, thereby improving understanding, ownership, and the quality
of the codesign outcomes.

One example of a process where academic knowledge is essential and cannot be
substituted by local knowledge is de novo design — the creation of novel crop-
ping or farming systems that break away from existing practices and explore new
possibilities (Meynard et al., 2012). Data-driven and model-based approaches
allow researchers to explore alternative designs systematically, which can then
feed into the codesign process in Living Labs.

This book presents a set of reproducible, data-driven workflows integrating
KoBoToolBox with R. With the use of examples, the book provides clear, step-
by-step guidance for researchers and practitioners on how data can be used to
enrich codesign processes. The book was compiled using the bookdown package
(Xie, 2025), which enables the seamless integration of R code, analytical results,
and narrative into a single, fully reproducible document.

A web version of this book is also available here: https://fbaudron.github.io/
book_supporting codesign/.

1.2 KoboToolBox and R

For the workflows in this book to be usable by the largest number, our ap-
proach is based on two open-source softwares: **KoboToolBox** (https://www.kobotoolbox.org/) and **R** (https://www.r-project.org/) operated through **R-Studio Desktop** (https://posit.co/download/rstudio-desktop/). We use these softwares because they are open-source, used by a large community of practice, and under constant development.

Digital surveys (e.g., using KoboToolBox) are much faster than surveys us-
ing pen and paper, reduce errors (no transcription and error trapping during
recording) and allow for the data to be immediately available right after it's
collected. Digital surveys also offer the possibility of associating geolocations
and multimedia data.

KoboToolBox is free and open access, and includes a free secured server, with
summary analysis and mapping functions (a major advantage of KoboToolBox
over alternatives, Figure 1.1). KoboToolBox offers a set of easy to use tools,
but complex logics can be built with them. KoboToolBox produces sharable
forms and is XLSForm compatible (XLSForm is a simple spreadsheet format
for designing digital data-collection forms). Data collection can be done offline.
Progress of data collection can be monitored without being on-site. KoboTool-
Box also provides a question library, and has a growing community of practice
and online help.

R is an open-source and free software environment (both a language and a
program). As a program language software, it offers replicability and shareabil-
ity (a major difference with ‘point-and-click’ softwares). R has a large online
<figure><img src="image_3.png" /><figcaption>Figure 1.1: Screenshot of the mapping functions of the KoboToolBox server</figcaption></figure>

community, is probably one of the most powerful computer environment for
data analysis currently available, and became very popular across many dis-
ciplines. Amongst other functions, R reads and writes files, manipulates and
summarizes data, runs statistical tests and models, and produces well-designed
publication-quality plots. The core functionalities of R are extended by hun-
dreds of packages (and growing). R interacts with other softwares, and the size
of datasets analyzed is only limited by the machine used.

1.3 Structure of this book

This book provides a collection of XLSForms (input to KoboToolBox) and well-
documented R scripts developed during the implementation of the CGIAR ini-
tiative 'Transformational Agroecology across Food, Land, and Water Systems'
(https://www.cgiar.org/initiative/agroecology/, referred to as 'CGIAR Agroe-
cology Initiative' in the rest of the document) in Zimbabwe, and that can be
adapted to other contexts and projects (Figure 1.2).
<figure><img src="image_4.png" /><figcaption>Figure 1.2: Workflow used in the CGIAR Agroecology initiative in Zimbabwe to support the codesign process with data, using two open-source softwares: KoboToolBox and R</figcaption></figure>

As such, the book is structured around five main chapters:

* Baselining: collecting data from farm surveys and subjecting them to summary statistics (Chapter 3)

* Understanding farm diversity and selecting representative farms: farm typologies and stratified samples (Chapter 4)

* Managing on-farm experiments: collecting and handling data efficiently (Chapter 5)

* Evaluating on-farm experiments in a participatory way: use of Shiny apps (Chapter 6)

* Scaling: technology extrapolation domains (Chapter 7)
1.4 Mbire District as a case study

This book uses data from Mbire District as a case study. Mbire District is located in Northern Zimbabwe (Figure 1.3) and is part of the Mid-Zambezi Valley, which is formed by the former floodplains of the Zambezi River between the Victoria Falls and Cabora Bassa Lake, at an average elevation of 400 m above sea level. The Mid-Zambezi Valley has a dry tropical climate, with low and very variable annual rainfall (on average between 450 and 650 mm per year) and a mean annual temperature of about 25°C. Two seasons are clearly defined: a rainy season from December to March and a long dry season from April to November.

Although the study area is communal land (i.e., state-owned land designated for small-scale family farming), it hosts significant wildlife populations (Gaidet et al., 2003), and is characterized by high biodiversity: over 40 large mammal, 200 bird and 700 plant species have been recorded (Coid et al., 2001). It is an area of great importance for the conservation of many emblematic African mammals including the African savanna elephant (Loxodonta Africana), the Cape buffalo (Syncerus caffer caffer), the African lion (Panthera leo), and the African wild dog (Lycaon pictus).

Mbire was one of the first districts in Zimbabwe to implement the world-
renowned Communal Area Management Programme for Indigenous Resources
(CAMPFIRE) in 1989 (Taylor, 2009). CAMPFIRE enables local communities,
through their local government, to manage wildlife as an economic asset for
rural development, and is considered one of the earliest programmes worldwide
expressing the ‘new conservation approach’ (Hulme and Murphree, 1999).

The rich biodiversity of Mbire is threatened by rapid land cover change that was earlier driven by cotton production (Baudron et al., 2011), and more recently by livestock production (Baudron et al., 2022). Wildlife imposes significant costs on local people in the form of crop destruction, livestock predation and human death and injuries, calling for urgent solutions (including landscape designs) that reduce these costs while maximizing benefits from nature.

Several projects implemented by CIRAD and/or CIMMYT have aimed at increasing agricultural productivity in the District with minimum negative consequences for biodiversity, including the project ‘Conservation de la Biodiversité dans la Vallée du Zambèze après Eradication de la Mouche Tse Tse’ (https://agritrop.cirad.fr/314708/1/ID314708.pdf) from 1997 to 2003, the project ‘Public-Private-Community Partnerships to improve food security and livelihoods in the South East Lowveld and Mid Zambezi Valley’ (https://umr-astre.cirad.fr/recherche/projets/parsel) from 2009 to 2011, the research program ‘Competing Claims on Natural Resources’ (https://www.competingclaims.nl/CC/Home.html) from 2006 to 2011, the CGIAR Agroecology Initiative from 2022 to 2024, and the CGIAR ‘Multifunc- tional Landscapes’ program (https://www.cgiar.org/cgiar-research-porfolio- 2025-2030/multifunctional-landscapes/) from 2025 to date.
<figure><img src="image_5.png" /><figcaption>Figure 1.3: Location of Mbire District on an orthographic projection of elevation and country boundaries</figcaption></figure>
Codes to produce the map of Figure 1.3 on R can be found here: https://github.com/FBaudron/Supporting-codesign/blob/main/01-introduction/map_site.R.
Chapter 2

Learning KoboToolBox and
R

2.1 A quick tutorial on KoboToolbox

To get familiar with KoboToolBox, follow the steps below.

2.1.1 Building a simple form manually

* Sign in: https://kf.kobotoolbox.org/accounts/login/?next=/#/;

* Click 'New' -> 'Build from scratch' -> Enter project name ('test') -> click 'Create project' (Sector and Country are compulsory: you can enter 'Humanitarian - Food Security' and 'Zimbabwe')

* Click '+' -> click 'Add question' -> click 'Text' -> replace 'New question' by 'Name of the head of the household'

* Click on 'Setting' (cog on the right side of the question box) -> click 'Yes' in 'Mandatory response'

* Click on the lower '+' sign -> click 'Add question' -> click 'Number -> replace 'Enter a number' by 'Age of the head of the household'. Make the response mandatory.

* Add a validation criterion. Click on 'Settings' -> click on 'Validation criteria' (left part of the frame) -> make the criteria 'Less than or equal to' 100.

* Click on the lower '+' sign -> click 'Add question' -> click 'Select one' -> replace 'New question' by 'Did you grow maize during the last season?' -> replace 'Option 1' by 'Yes' and 'Option 2' by 'No'. Make the response mandatory.

* Click on the lower '+' sign -> click 'Add question' -> click 'Date' -> replace 'Enter a date' by 'Planting date'.
* Add a skip logic. Click on 'Settings' -> click on 'Skip logic' (left part of the frame) -> click on 'Add a condition' -> select 'Did you grow maize during the last season?' -> keep the logic 'Yes'

* Click on the lower '+' sign -> click 'Add question' -> click 'Ranking' -> replace 'New question' by 'Reason to grow maize' -> add a second choice -> Make 3 items to be ranked: food, income, and feed. Add the same skip logic as for 'Planting date' (only if the response to 'Did you grow maize during the last season?' is 'Yes'

* Click on the lower '+' sign -> click 'Add question' -> click 'Point'

* Click on the lower '+' sign -> click 'Add question' -> click 'Acknowledge' -> replace 'Acknowledge' by 'Do you agree to be part of this survey?'. Drag this question so it becomes the first question.

* Click 'Save' (top right corner)

* Preview the form (eye icon, top left corner)

* Close the form. And deploy the form (if you don't see the function to deploy, retrieve the test form in 'Drafts')

### 2.1.2 Filling forms on your android device

* Follow the steps here: https://support.kobotoolbox.org/data_collection_kobocollect.html

* Use this URL: https://kc.kobotoolbox.org, your user name and your password

* Go to 'Get a Blank Form', and get the form 'test'

* Go to 'Fill a Blank Form', select your form, fill it up, and save

### 2.1.3 Visualizing and downloading the collected data on your computer

* Go to 'Summary'. Check that your entry(ies) is(are) there.

* Go to 'Data'. Inspect the data.

* Go to 'Map' (top left of the 'Data' tab, bottom of the menu)

* Go to 'Download' -> Click on 'Export' -> Click on 'Download' in the section 'Exports'. You have your data ready for analysis!

### 2.1.4 Saving a PDF of your questionnaire to share with colleagues

* Click on form preview (eye icon, top left corner)

* Click on 'Print'

* Change 'Destination' into 'Save as PDF'
### 2.1.5 Working with an XLSForm

* Click on 'Form' -> Click on 'More action' ('It') -> Click on 'Download XLS'.

* Take a moment to study the content of the tab 'survey' and the tab 'choices'. The tab 'survey' contains all your questions (in the order they appear in your form). The tab 'choices' contains all the factorial choices for your questions 'Select one'

* Save the XLS Form as 'test2.xls'

* Change the list name for the question 'Did you grow maize last season' into 'maize', in both the 'survey' and the 'choices' tabs.

* Change the list name for the question 'Reason to grow maize' into 'reason', in both the 'survey' and the 'choices' tabs.

* In the 'survey' tab, make all the answers as 'required' (except for geopoint) by changing 'false' into 'true'

* In the 'survey' tab, column 'name', simplify variable names (be careful spaces are not allowed!) into: 'Name', 'Age', 'Maize', 'Pdate', 'Reason', 'Reason_label'. Remember to make the same changes in the 'relevant' column (e.g., change 'Did_you_grow_maize_during_the_' to 'Maize')

* Create a question 'Sex of the head of the household' under the question 'Name of the head of the household'. In the 'survey' tab, insert a row, enter 'select_one sex' in 'type', 'Sex' in 'name', 'Sex of the head of the household' in 'label', 'true' in 'required'. In the 'choices' tab, insert two rows with 'sex' in 'list_name' for both, and 'male' - 'Male' for 'name' and 'label' for the first one, and 'female' - 'Female' for 'name' and 'label' for the second one.

* Save

* Go on your KoBoToolbox account and click on 'New'

* Click on 'Upload an XLSForm'

* Upload the file you created 'test2.xls'

* Click on 'Create project' (keep the name 'test2')

* Preview the form. You can also go on 'Edit' (pencil icon) to see the changes that have been made to the form

### 2.1.6 Adding cascading questions

* Save 'test2.xls' as 'test3.xls'.

* In the 'survey' tab, add three rows before the question 'Name of the head of the household'

* In the first row, enter 'select_one country' (as type), 'Country' (as name), 'Country' (as label), 'true' (in required). Be careful with the designation of the list_name (no capital C) and the designations of the name and the label (capital C).

* In the second row, enter 'select_one province' (as type), 'Province' (as name), 'Province' (as label), 'true' (in required). Be careful with the
designation of the list_name (no capital P) and the designations of the
name and the label (capital P).

* In the third row, enter 'select_one district' (as type), 'District' (as name), 'District' (as label), 'true' (in required). Be careful with the designation of the list_name (no capital D) and the designations of the name and the label (capital D).

* Add a column 'choice_filter' (after the column 'required'). In this column, enter `Country=${Country}` for the row on Province, and `Province=${Province}` for the row on District (nothing for the row on Country).

* In the 'choices' tab, add two columns 'Country' and 'Province' (after label)

* Add 14 rows corresponding to the 2 countries (Zimbabwe, Zambia), the 4 provinces (Manicaland, Matebeleland, Southern, Northern) and the 8 districts (Chipinge, Nyanga, Hwange, Beitbridge, Mazabuka, Monze, Mbala, Senga). Enter the right list_name for each (country, province, district), names (with capital letters), and label (with capital letters). For Province names, enter the corresponding Country in the column 'Country' you created. For District names, enter the corresponding Province name in the column 'Province' you created.

* Save

* Go on your KoBoToolbox account and click on 'New'

* Click on 'Upload an XLSForm'

* Upload the file you created 'test3.xls'

* Click on 'Create project' (keep the name 'test3')

* Preview the form. You can also go on 'Edit' (pencil icon) to see the changes that have been made to the form

2.1.7 An example of survey

Once you are familiar with the basics of KoboToolbox, you can analyse a more
complex XLSForm.

The XLSForm used to build the survey questionnaire in KoboToolBox for the
baseline survey of the CGIAR Agroecology Initiative in Zimbabwe can be found
here: https://github.com/FBaudron/Supporting-codesign/blob/main/02-
learning-kobo-r/BASELINE%20SURVEY%20TOOL.xlsx. The PDF of the
questionnaire is also provided: https://github.com/FBaudron/Supporting-
codesign/blob/main/02-learning-kobo-r/BASELINE%20SURVEY%20TOOL.
pdf. Feel free to adapt this XLSForm to your own use.

We will use a sample of the data generated from this questionnaire in the next
chapter.
## 2.2 Resources to learn R

To work with R (https://www.r-project.org/) on your own computer, you need to download the program and install it, and install R-Studio Desktop (https://posit.co/download/rstudio-desktop/), which is a separate program that makes R easier to use.

There are many free resources on the web to learn R:

* R for Beginners by Emmanuel Paradis (https://cran.r-project.org/doc/contrib/Paradis-rdebuts_en.pdf)

* R tutorial by Kelly Black (https://www.cyclismo.org/tutorial/R/)

* Introduction to R by Ross Ihaka (https://www.stat.auckland.ac.nz/~ihaka/120/Notes/ch02.pdf)

* The Art of R Programming (https://nostarch.com/artofr.htm)

* YouTube Playlist (https://www.youtube.com/ playlist?list=PLOU2XLYxmsIK9qQfztXeybpHvru-TrqAP)
# Chapter 3

## Baselining

In this chapter, we will go through steps to prepare the data collected from a KoboToolBox survey for analysis, and run preliminary analyses (mainly figures and tables of summary statistics) using R.

Download a sample of the data generated from the questionnaire presented in the previous chapter here: https://github.com/FBaudron/Supporting-codesign/blob/main/03-baselining/BASELINE%20DATA.xlsx. This dataset will be used in the following workflow.

A standalone R script containing the full workflow is available at the following link: https://github.com/FBaudron/Supporting-codesign/blob/main/03-baselining/baseline.R.

### 3.1 Data preparation

We begin by clearing the environment, ensuring all necessary packages are installed and loaded, and importing the data.

```r
rm(list = ls())
```

```r
if (!require("openxlsex")) install.packages("openxlsex")
if (!require("janitor")) install.packages("janitor")
if (!require("dplyr")) install.packages("dplyr")
if (!require("geodata")) install.packages("geodata")
if (!require("terra")) install.packages("terra")
if (!require("ggspatial")) install.packages("ggspatial")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("tidyterra")) install.packages("tidyterra")
if (!require("data.table")) install.packages("data.table")
```
```r
if (!require("tidyr")) install.packages("tidyr")
if (!require("gtsummary")) install.packages("gtsummary")
if (!require("kableExtra")) install.packages("kableExtra")
```

```r
library(openxlxs)
library(janitor)
library(dplyr)
library(geodata)
library(terra)
library(ggspatial)
library(ggplot2)
library(tidyterra)
library(data.table)
library(tidyrr)
library gtsummary
library(kableExtra)
```

```r
data = read.xlsx("BASELINE DATA.xlsx")
```

### 3.1.1 Data cleaning

This section provides a number of steps to clean the data extracted from KoboToolBox and produce a dataset that can be used for analysis.

We first pass the dataframe through the function *clean_names()* of the R package **janitor** (Firke, 2024), which makes names unique, converts names to lower cases, replaces spaces, special characters, and punctuation with underscores, etc.

```r
data = clean_names(data)
```

District and wards were recorded separately for farms previously interviewed and for farms newly interviewed. Therefore, we merge the 2 District variables and the 2 Ward variables, keep the merged columns only, and remove the row 'Has this farm already been interviewed earlier or is this a new farm?'

```r
data$district = ifelse(is.na(data$district), data$district_2,
  data$district)
data$ward = ifelse(is.na(data$ward), data$ward_2, data$ward)
data = data[, -c(3, 6:7)]
```

We merge the factors 'Ward 3' and 'Ward 3 9 12 17' of the variable 'ward'.
```R
data$ward = ifelse(data$ward == "Ward 3", "Ward 3 9 12 17", data$ward)
table(data$ward)
```

Ward 2 Ward 3 9 12 17
134 136

For some variables, we extract the last part of the variable name and make it the new variable name.

```R
# names(data)[c(7:10, 54:62, 68:85, 95:112)]
names(data)[c(7:10, 54:62, 68:85, 95:112)] = sub(".*_", "",
  names(data)[c(7:10, 54:62,
    68:85, 95:112)])
names(data)[c(7:10, 54:62, 68:85, 95:112)]
```

```
[1] "latitude" "longitude" "altitude" "precision" "bean"
[6] "cowpea" "groundnut" "lablab" "pigeonpea" "bean"
[11] "pumpkin" "watermelon" "cucumber" "maize" "bean"
[16] "potato" "tomato" "onion" "garlic" "cabbage"
[21] "rape" "kovo" "carrot" "occra" "watermelon"
[26] "cucumber" "butternut" "cane" "banana" "pawpaw"
[31] "nartje" "maize" "bean" "potato" "tomato"
[36] "onion" "garlic" "cabbage" "rape" "kovo"
[41] "carrot" "occra" "watermelon" "cucumber" "butternut"
[46] "cane" "banana" "pawpaw" "nartje"
```

We correct the names for those variables with a composed name, to avoid problems later.

```R
names(data)[54] = "common"
names(data)[59] = "velvet bean"
names(data)[69] = "common bean"
names(data)[70] = "sweet potato"
names(data)[82] = "sugar cane"
names(data)[96] = "common bean"
names(data)[97] = "sweet potato"
names(data)[109] = "sugar cane"
```

We rename manually other long variables.
## 3.1. DATA PREPARATION

19

```r
names(data)[c(33:49, 65, 92, 173, 177, 187:188, 210, 243)] =
  c("cropped_area", "fallow",
    "maize", "sorghum", "pearl_millet", "finger_millet", "sesame",
    "groundnut", "tobacco",
    "cotton", "potato", "sweet_potato", "sunflower", "soybean",
    "common bean", "cowpea",
    "roundnut", "ind_garden", "com_garden", "donkeys", "pigs",
    "bee_hives", "fish_pond",
    "main_food_source", "main_income_source")
```

We turn No/Yes columns into 0/1.

```r
# str(data[c(65, 92, 154:168, 188, 192:208, 225:241)])
data[c(65, 92, 154:168, 188, 192:208, 225:241)] = lapply(data[c(65, 92,
  154:168,
  187, 192:208, 225:241)], function(x) ifelse(x == "Yes", 1, 0))
# str(data[c(65, 92, 154:168, 188, 192:208, 225:241)])
```

We turn some continuous variables into binary ones (0/1).

```r
data$bee_hives = ifelse(data$bee_hives > 0, 1, 0)

table(data$bee_hives)
```

```
0   1
242  28
```

### 3.1.2 Calculation of new variables

After cleaning the dataset, we calculate a number of variables we will use in further analyses.

First, we calculate the number of elderly family members, the number of adults, the number of young family members, and the family size.

```r
data$elderly = data$total_number_of_adult_males_aged_61_or_more +
  data$total_number_of_adult_females_aged_61_or_more

data$adult = data$total_number_of_adult_males_aged_25_60 +
  data$total_number_of_adult_females_aged_25_60

data$youth = data$total_number_of_young_males_aged_15_25 +
  data$total_number_of_young_females_aged_15_25
```
```
data$family_size = data$elderly + data$adult + data$youth +
  data$total_number_of_children_and_teens_aged_3_14 +
  data$total_number_of_infants_of_age_0_to_2
```

We then calculate the total equipment value, assuming a value of 10,000 USD for a tractor, 2,000 USD for a two-wheel tractor, 100 USD for a plough, 150 USD for a cultivator, 500 USD for a scotch cart, 50 USD for a wheelbarrow, 50 USD for a solar panel, and 300 USD for a motorized grinding mill.

```
data$eq_value = 10000 * data$nb_of_tractors + 2000 *
  data$nb_of_twoWheel_tractors +
  100 * data$nb_of_ploughs + 150 * data$nb_of_cultivators + 500 *
  data$nb_of_scotchcarts +
  50 * data$nb_of_wheelbarows + 50 * data$nb_of_solar_panels + 300 *
  data$nb_ofmotorized_grinding_mill

summary(data$eq_value)
```

```
Min. 1st Qu.  Median   Mean 3rd Qu.  Max.
0.0    50.0   100.0   360.4   337.5 20000.0
```

We calculate the area cropped in small grains, in legume crops and in other crops.

```
data$small_grains_area = data$sorghum + data$pearl_millet +
  data$finger_millet

data$legumes = data$groundnut + data$soybean + data$common_bean +
  data$cowpea + data$roundnut

data$other = data$sesame + data$tobacco + data$cotton + data$potato +
  data$sweet_potato +
  data$sunflower
```

We ensure that the sum of the areas cropped in maize, in small grains, in legumes crops and in other crops area doesn't exceed the total cropped area.

```
data$other = ifelse(data$other > data$cropped_area - data$maize -
  data$small_grains_area -
  data$legumes, data$other, data$cropped_area - data$maize -
  data$small_grains_area -
  data$legumes)
```

We calculate the percentage of total land operated represented by fallow land.
data$prop_fallow = data$fallow * 100/(data$cropped_area + data$fallow + 
\rightarrow 0.001)

summary(data$prop_fallow)

<table><thead><tr><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th></tr></thead><tbody><tr><td>0.00</td><td>0.00</td><td>0.00</td><td>20.01</td><td>39.99</td><td>99.98</td></tr></tbody></table>

We calculate the percentage of cropped area occupied by non-cereal crops.

data$prop_non_cereals = (data$legumes + data$other) * 100/(data$maize + 
\rightarrow data$small_grains_area + 
data$legumes + data$other + 0.001)

summary(data$prop_non_cereals)

<table><thead><tr><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th></tr></thead><tbody><tr><td>0.00</td><td>0.00</td><td>24.97</td><td>29.45</td><td>49.98</td><td>99.97</td></tr></tbody></table>

We calculate the total cereal production as the sum of maize production and sorghum production, which are the main cereal crops in the area.

data[c(123, 136)][is.na(data[c(123, 136)])] = 0

data$cerealProd = data$total_m maize production in the season 2021_22_kg + 
\rightarrow data$total_sorghum production in the season 2021_22_kg

summary(data$cerealprod)

<table><thead><tr><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th></tr></thead><tbody><tr><td>0.0</td><td>250.0</td><td>600.0</td><td>859.2</td><td>1087.5</td><td>6000.0</td></tr></tbody></table>

We calculate the total quantity of fertilizer used, and the total quantity of organic amendments - manure and compost - used.

data[c(88:91, 119:123, 132:153)][is.na(data[c(88:91, 119:123, 132:153)])] 
\rightarrow 0

data$fertilizers =

\rightarrow data$quantity_kg_of_basal_fertilizer_applied_to_cotton_in_theSeason_2021_22_kg +

\rightarrow data$quantity_kg_of_basal_fertilizer_applied_toroundnut_in_theSeason_2021_22_kg +
→ data$quantity_kg_of_basal_fertilizer_applied_to_maize_in_theSeason_2021_22
+

→ data$quantity_kg_of_basal_fertilizer_applied_to_sesame_in_theSeason_2021_22
+

→ data$quantity_kg_of_basal_fertilizer_applied_to_sorghum_in_theSeason_2021_22
+

→ data$quantity_kg_of_basal_fertilizer_applied_to_the_garden_in_the_last_12Months
+

→ data$quantity_kg_of_topdressing_fertilizer_applied_to_cotton_in_theSeason_2021_22
+

→ data$quantity_kg_of_topdressing_fertilizer_applied_to_groundnut_in_theSeason_2021_22
+

→ data$quantity_kg_of_topdressing_fertilizer_applied_to_maize_in_theSeason_2021_22
+

→ data$quantity_kg_of_topdressing_fertilizer_applied_to_sesame_in_theSeason_2021_22
+

→ data$quantity_kg_of_topdressing_fertilizer_applied_to_sorghum_in_theSeason_2021_22
+

→ data$quantity_kg_of_topdressing_fertilizer_applied_to_the_garden_in_the_last_12Months

data$manure =

→ data$quantity_kg_of_manure_applied_to_cotton_in_theSeason_2021_22 +
  data$quantity_kg_of_manure_applied_to_groundnut_in_theSeason_2021_22

→ + data$quantity_kg_of_manure_applied_to_maize_in_theSeason_2021_22 +
  data$quantity_kg_of_manure_applied_to_sesame_in_theSeason_2021_22 +

→ data$quantity_kg_of_manure_applied_to_sorghum_in_theSeason_2021_22 +

→ data$quantity_kg_of_manure_applied_to_the_garden_in_the_last_12Months

data$compost =

→ data$quantity_kg_of_compost_applied_to_cotton_in_theSeason_2021_22 +

→ data$quantity_kg_of_compost_applied_toonde nut_in_theSeason_2021_22

→ + data$quantity_kg_of_compost_applied_to_maize_in_theSeason_2021_22

→ +
  data$quantity_kg_of_compost_applied_to_sesame_in_theSeason_2021_22 +

→ data$quantity_kg_of_compost_applied_to_sorghum_in_theSeason_2021_22

+

→ data$quantity_kg_of_compost_applied_to_the_garden_in_the_last_12Months
```
data$amendment = data$manure + data$compost
```

We calculate the total number of cattle, goats, small ruminants (sheep and goats) and poultry owned.

```
data$cattle = data$howmany_adult_indigenous_cattle_do_you_currently_own
  + data$howmany_improved_beef_cattle_do_you_currently_own +
  data$howmany_improved_dairy_cattle_do_you_currently_own +
  data$howmany_juvenile_cattle_3_yr_old_indigenous_and_improved_do_you_currently_own

data$goats = data$howmany_indigenous-goats_do_you_currently_own +
  data$howmany_improved-goats_do_you_currently_own

data$small_rum = data$goats + data$howmany_sheep_do_you_currently_own

data$poultry = data$howmany_indigenous_chickens_do_you_currently_own +
  data$howmanyсыlayer_chickens_do_you_currently_own +
  data$howmany_broiler_chickens_do_you_currently_own +
  data$howmany_guinea_fowls_do_you_currently_own +
  data$howmany_turkeys_do_you_currently_own +
  data$howmany_ducks_do_you_currently_own +
  data$howmany_geese_do_you_currently_own
```

We convert livestock numbers into Tropical Livestock Units, using a 250 kg live weight value for one TLU (Houérou and Hoste, 1977), and assuming poultry to be equivalent to 0.01 TLU, sheep and goats to 0.1 TLU, pigs to 0.2 TLU, donkeys to 0.5 TLU, and all types of cattle to 0.7 TLU (Jahnke, 1982).

```
data$TLU = 0.7 * data$cattle + 0.5 * data$donkeys + 0.1 * data$small_rum
  + 0.2 *
  data$pigs + 0.01 * data$poultry
```

```
summary(data$TLU)
```

<table><thead><tr><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th></tr></thead><tbody><tr><td>0.000</td><td>0.200</td><td>1.280</td><td>4.264</td><td>5.588</td><td>69.590</td></tr></tbody></table>

Similarly, we calculate the total livestock offtake in TLU.

```
data[c(189:191)][is.na(data[c(189:191)])] = 0

data$offtake = 0.7 *
  data$number_of_cattle_sold_or_slaughtered_in_the_past_12Months +
  0.1 * data$number_of-goats_sold_or_slaughtered_in_the_past_12months
  + 0.01 *
```
```
-> data$number_of_indigenous_chicken_sold_or_slaughtered_in_the_past_12Months
summary(data$offtake)
```

```
Min. 1st Qu.  Median  Mean 3rd Qu.  Max.
0.0000  0.0000  0.0300  0.3817  0.3000 16.0000
```

We calculate the consumption (0/1) of vegetables, fruits and meat in the 24 hours preceding the interview, and calculate a Household Dietary Diversity Score (HDDS) based on 12 food groups: (1) cereals; (2) white tubers and roots; (3) vegetables; (4) fruits; (5) meat; (6) eggs; (7) fish and other seafood; (8) legumes, nuts and seeds; (9) milk and milk products; (10) oils and fats; (11) sweats; and (12) spices, condiments and beverages (Kennedy et al., 2010).

We replace zero values for HDDS (an impossibility) with missing values.

```
data$vegetables = data$vitamin_a_rich_vegetables_and_tubers +
  data$dark_green_leafy_vegetables +
  data$other_vegetables
data$vegetables = ifelse(data$vegetables > 0, 1, 0)
```

```
data$fruits = data$vitamin_a_rich_fruits + data$other_fruits
data$fruits = ifelse(data$fruits > 0, 1, 0)
```

```
data$meat = data$organ_meat + data$flesh_meat
data$meat = ifelse(data$meat > 0, 1, 0)
```

```
data$hdds = data$cereals + data$white_roots_and_tubers + data$vegetables
  + data$fruits +
  data$meat + data$eggs + data$fish_and_seafood +
  data$legumes_nuts_and_seeds +
  data$milk_and_milk_products + data$oils_and_fats + data$sweets +
  data$spices_condiments_beverages
```

```
data[which(data$hdds == 0), ] = NA
summary(data$hdds)
```

```
Min. 1st Qu.  Median  Mean 3rd Qu.  Max.  NA's
1.000  2.000  3.000  3.147  4.000  8.000  4
```

We calculate the proportion of the year the monthly food security is moderate or high, as a proxy of food security.
```
data$food_security = rowMeans(data[212:223] != "Low")
```

```
summary(data$food_security)
```

<table><thead><tr><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th><th>NA's</th></tr></thead><tbody><tr><td>0.0000</td><td>0.4167</td><td>0.5833</td><td>0.5583</td><td>0.7500</td><td>1.0000</td><td>4</td></tr></tbody></table>

We then calculate crop diversity. We need to ensure that crop species that may have been grown as a main crop and an intercrop, and/or an individual garden crop, and/or a communal garden crop are only counted once.

We first add columns for crops that were listed as additional, in case they were not already included as crop variables.

```
data$cowpea.w =
```

```
  ifelse(data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Cowpeas" |
```

```
  data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Cowpeas, lablab", 1, 0)
```

```
data$lablab.w =
```

```
  ifelse(data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Cowpeas, lablab", 1, 0)
```

```
data$runinga.w =
```

```
  ifelse(data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Runinga, rossella", 1, 0)
```

```
data$rosella.w =
```

```
  ifelse(data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"rosella" |
```

```
  data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Rosella" |
```

```
  data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"rosella" |
```

```
  data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Rosella" |
```

```
  data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```

```
"Rossella" |
```

```
  data$list_other_main_crops_you_grew_during_theSeason_2021_22_separated_by_commas
```

```
==
```
```
"Rotella" |
    data$list_other_main_crops_you_grew_during_the season_2021_22_separated_by_commas
==
"Runinga, rossella", 1, 0)
```

```
data$maize.x = ifelse(data$list_other_intercrops_separated_by_commas ==
    "Maize,rosella",
    1, 0)
```

```
data$rosella.x = ifelse(data$list_other_intercrops_separated_by_commas ==
    "Maize,rosella",
    1, 0)
```

```
data$sweet_sorghum.x =
    ifelse(data$list_other_intercrops_separated_by_commas ==
        "Sweetsoughum,milet",
        1, 0)
```

```
data$pearl_millet.x =
    ifelse(data$list_other_intercrops_separated_by_commas ==
        "Sweetsoughum,milet",
        1, 0)
```

```
data$chilli.y =
    ifelse(data$list_other_crops_you_grew_in_your Brady's garden separated by commas
        ==
        "Chill,potato", 1, 0)
```

```
data$potato.y =
    ifelse(data$list_other_crops_you_grew_in_your Brady's garden separated by commas
        ==
        "Chill,potato", 1, 0)
```

```
data$katarina.y =
    ifelse(data$list_other_crops_you_grew_in_your Brady's garden separated by commas
        ==
        "Katarina" |
        data$list_other_crops_you_grew_in_your Brady's garden separated by commas
        ==
        "Sugarloaf, katarina", 1, 0)
```

```
data$muboora.y =
    ifelse(data$list_other_crops_you_grew_in_your Brady's garden separated by commas
        ==
        "Muboora", 1, 0)
```

```
data$shoeshine.y =
    ifelse(data$list_other_crops_you_grew_in_your Brady's garden separated by commas
        ==
```
"Shoeshine", 1, 0)

data$sugarloaf.y =
  ifelse(data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf, katarina" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf, tsunga" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Tsunga,sugarloaf", 1, 0)

data$tsunga.y =
  ifelse(data$list_other_crops_you_grew_in_your Brady
  ==
  "Sugarloaf, tsunga" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Tsunga" |
  data$list_other_crops_you_grew_in_your Brady
  ==
  "Tsunga,sugarloaf", 1, 0)

We replace missing values by zeros.

We then clean crop species names to strip off dots and numbers at the end.

We collapse duplicates and retain the maximum value across columns of the
same crop.

We finally compute the crop diversity as the number of crop species present.

data[c(35:49, 54:62, 68:85, 95:112, 280:294)][is.na(data[c(35:49, 54:62,
  68:85, 95:112,
  280:294)]] = 0
```
base_names = sub("\\..*", "", names(data)[c(35:49, 54:62, 68:85, 95:112,
  280:294)])
```

```
collapsed = data[c(35:49, 54:62, 68:85, 95:112, 280:294)] %>%
  mutate(row_id = row_number()) %>%
  group_by(row_id) %>%
  summarise(across(everything(), ~.x, .names = "{.col}") %>%
  select(-row_id)
```

```
collapsed = as.data.frame(lapplysplit.default(data[c(35:49, 54:62,
  68:85, 95:112,
  280:294)], base_names), function(df) rowSums(df > 0)))
```

```
data$crop_diversity = rowSums(collapsed > 0)
```

```
summary(data$crop_diversity)
```

```
Min. 1st Qu.  Median   Mean 3rd Qu.  Max.
0.000    2.000    4.000    4.348    6.000   14.000
```

Livestock diversity can be calculated in a similar, but simpler way (as there is only one variable per livestock species)

```
data[c(269:270, 272, 173, 177, 181:188)][is.na(data[c(269:270, 272, 173,
  177, 181:188)])] = 0
```

```
data$livestock_diversity = rowSums(data[c(269:270, 272, 173, 177,
  181:188)] > 0)
```

```
summary(data$livestock_diversity)
```

```
Min. 1st Qu.  Median   Mean 3rd Qu.  Max.
0.000    1.000    2.000    2.396    3.000    8.000
```

We calculate the total diversity of species present on the farm, as the sum of crop diversity and livestock diversity.

```
data$tot_div = data$crop_diversity + data$livestock_diversity
```

```
summary(data$tot_div)
```

```
Min. 1st Qu.  Median   Mean 3rd Qu.  Max.
0.000    4.000    6.000    6.744    9.000   22.000
```

Create a small dataset
small_dataset = data[, c(253, 4, 3, 11:13, 257, 33:34, 263, 265, 268:269,
  271:272,
  187, 65, 92, 210, 243, 279, 278, 155, 157, 159:164, 167:168, 202:205,
  297, 264,
  274, 258)]

saveRDS(small_dataset, "small_dataset.rds")

3.2 Preliminary analyses

3.2.1 Map of surveyed farms

First, we download the shapefile of Zimbabwe, at the administrative level 3 (i.e.,
at the ward-level). Such shapefiles can be downloaded through the function
gadm() of the R package geodata (Hijmans et al., 2024).

zwe3 = gadm(country = "ZWE", level = 3, path = "Geodata")
plot(zwe3)

<figure><img src="image_6.png" /></figure>

We then subset this shapefile to the district of Mbire.

mbire = zwe3[zwe3$NAME_2 == "Mbire", ]
plot(mbire)
<figure><img src="image_7.png" /></figure>

Next, we download a raster of elevation, which we will use as background for
the map, again from the R package `geodata`.

We download the tile corresponding to the study area, using the coordinates of
the center of the extent of Mbire.

We give the name ‘elevation’ to that raster.

```
elevation = elevation_3s(lon = (30.04626 + 31.14537)/2, lat = (-16.41686
  ↪ + -15.61071)/2,
  path = "Geodata")
```

```
names(elevation) = "elevation"
```

```
plot(elevation)
```
<figure><img src="image_8.png" /></figure>

We use the R package **terra** (Hijmans, 2024) for all spatial calculations below.

We use EPSG:4326 coordinate system, and crop our raster to the extent of the
area that was surveyed.

elevation = project(elevation, "EPSG:4326")

elevation = crop(elevation, ext(30.2, 30.65, -16.23, -15.97))

plot(elevation)

<figure><img src="image_9.png" /></figure>
We use the coordinates in the dataframe to produce a shapefile of farm locations.

We ensure to use the same coordinate system as for the elevation raster.

```r
dataVect = vect(data, geom = c("longitude", "latitude"), crs =
  "EPSG:4326")
```

We extract a shapefile of the 5 wards included in the survey.

```r
wards = mbire[mbire$NAME_3 == "Ward 2" | mbire$NAME_3 == "Ward 3" |
  mbire$NAME_3 ==
    "Ward 9" | mbire$NAME_3 == "Ward 12" | mbire$NAME_3 == "Ward 17", ]
plot(wards)
```

<figure><img src="image_10.png" /></figure>

We get the coordinates of the centroid of each of these 5 wards, to use them in our map for the labels of the wards.

We jitter the coordinates of the centroid of Ward 2 (to avoid overlap with surveyed points).

```r
ward_centroids = centroids(wards)
ward_centroids = data.frame(x = crds(ward_centroids)[, 1], y =
  crds(ward_centroids)[,
    2], ward = wards$NAME_3)
ward_centroids$y[ward_centroids$ward == "Ward 2"] =
  ward_centroids$y[ward_centroids$ward ==
```
"Ward 2"] - 0.05

ward_centroids$x[ward_centroids$ward == "Ward 2"] =
  ward_centroids$x[ward_centroids$ward ==
    "Ward 2"] + 0.05

We enter manually coordinates for a label "MOZAMBIQUE".

mozambique = data.frame(x = 30.55, y = -15.99, label = "MOZAMBIQUE")

We plot the map using the R package **ggplot2** (Wickham, 2016). Before that,
we convert the elevation raster into a dataframe.

We use `annotation_scale()` and `annotation_northArrow()` from the R package `ggspatial` (Dunnington, 2023).

elevation DF = as.data.frame(elevation, xy = TRUE)
ggplot() + themeBW() + geom_raster(data = na.omit(elevation DF), aes(x =
  x, y = y,
  fill = elevation)) + geom_spatvector(data = mbire, fill = NA,
  linewidth = 1,
  color = "grey30") + geom_spatvector(data = data vect, shape = 16,
  color = "#781C6DFF",
  alpha = 0.5, size = 3) + geom_text(data = ward_centroids, aes(x = x,
  y = y, label = ward),
  color = "black", size = 4, fontface = "bold") + geom_text(data =
  mozambique,
  aes(x = x, y = y, label = label), color = "black", size = 5, fontface
  = "italic") +
  scale_fill_distiller(palette = "RdYlGn", limits = c(350, 650), breaks
  = seq(350,
    650, by = 75)) + scale_y continuous(breaks = seq(-16.22, -15.97,
  by = 0.1)) +
  coord.sf(xlim = c(30.2, 30.65), ylim = c(-15.97, -16.23), expand =
  FALSE) + labs(fill = "Elevation (m.a.s.l.)") +
  theme plot. title = elementDA blank(), axis. title = elementDA blank(),
  axis. text.y = element_text(angle = 90,
    hjust = 0.5, size = 12), axis.text.x = element_text(size = 12),
  legend. title = element_text(size = 14),
  legend.text = element_text(size = 12), legend. position = "right",
  plot. margin = margin(10,
    10, 10, 10)) + annotation_scale(location = "bl", width. hint =
  0.5, height = unit(0.25,
    "cm"), text_cex = 1, pad_x = unit(0.15, "cm"), pad_y = unit(0.15,
    "cm")) + annotation_northArrow(location = "tl",
    which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"),
    pad_x = unit(0.15,
```
"cm"), pad_y = unit(0.15, "cm"), style =
  northArrow_fancy_orienteering(text_size = 10))
```

<figure><img src="image_11.png" /></figure>

Next, we present key visualizations of survey data: stacked bar chart, lollipop diagram, bubble chart, and pie chart, all using the R package `ggplot2`.

### 3.2.2 Stacked bar chart

We will plot perceptions of trends in soil compaction, soil erosion, soil fertility and tree cover.

We create a subset of the dataset, rename variables to have nice labels for the plot, and remove empty rows.

```
land_mgt = data[, c(115:118)]
```

```
names(land_mgt) = c("Soil fertility", "Soil erosion", "Soil compaction",
  "Tree cover")
```

```
land_mgt = na.omit(land_mgt)
```

We transform the dataset from wide to long format using the function `gather()` from the R package `tidyr` (Wickham et al., 2024), calculate percentages for each characteristic in each status ('deteriorating', 'stable', 'improving'), and reorder the levels of the factor 'status' (from 'deteriorating' to 'improving').
```r
land_mgt = gather(land_mgt, Characteristic, Status, "Soil
  -> fertility":"Tree cover")
```

```r
land_mgt = land_mgt %>%
  group_by(Characteristic, Status) %>%
  summarise(n = n()) %>%
  mutate(Percentage = n/sum(n) * 100) %>%
  select(-n) %>%
  ungroup()
```

```r
land_mgt$Status = factor(land_mgt$Status, levels = c("Improving",
  -> "Stable", "Deteriorating"))
```

Finally, we plot the stack bar chart.

Note the use of `coord flip()`.

```r
ggplot(data = land_mgt, aes(x = Characteristic, y = Percentage, fill =
  -> Status)) +
  geom_bar(stat = "identity", width = 0.7, color = "black") +
  -> coord flip() + theme_bw() +
  scale_fill_manual(values = c("#FCFFA4FF", "#ED6925FF", "#781C6DFF"))
  -> theme(axis.title.x = element_text(size = 14,
  face = "bold", margin = margin(10, 0, 0, 0)), axis.title.y =
  -> elementBlank(),
  axis.text.x = element_text(size = 10, face = "bold"), axis.text.y =
  -> element_text(size = 12,
    face = "bold"), legend.title = element_text(size = 14, face =
  -> "bold"), legend.text = element_text(size = 12))
```

<figure><img src="image_12.png" /></figure>
### 3.2.3 More complex stacked bar chart

We will plot crop area distribution per farm, for all farm surveyed, ordering farms by decreasing total cropped area.

We create a subset of the dataset, correct an abnormal value for area cropped in legume, calculate total cropped area, order the dataset by decreasing cropped area and add an ID.

```r
crop = data[, c(35, 259:261)]
crop$legumes = ifelse(crop$legumes > 5, 0, crop$legumes)
crop$cropped_area = crop$maize + crop$small_grains_area + crop$legumes +
  ↪️ crop$other
crop = crop[order(-crop$cropped_area), ]
crop$ID = 1:nrow(crop)
```

Next, we apply a rolling average to the different variables we will use in the plot with subsets of 10 farms to smooth the curves for better visualization/easier interpretation.

```r
crop$maize = frollmean(crop$maize, 10)
crop$small_grains_area = frollmean(crop$small_grains_area, 10)
crop$legumes = frollmean(crop$legumes, 10)
crop$other = frollmean(crop$other, 10)
```

Next, we transform the dataset from wide to long format, and we reorder the levels of the factor crop, for 'maize' to appear at the bottom of our plot, and 'other crops' at the top.

```r
crop = gather(crop, crop, area, "maize":"other")
crop$crop = factor(crop$crop, levels = c("maize", "small_grains_area",
  "legumes", "other"))
```

Finally, we plot the stack bar chart.

```r
crop %>%
  ggplot(aes(x = ID, y = area, fill = as.factor(crop))) + geom_bar(stat
  = "identity",
  width = 1, alpha = 0.9, position = position_stack(reverse = TRUE)) +
  scale_fill_manual(values = c("#FCA50AFF",
```
```r
"#DD513AFF", "#932667FF", "#420A68FF"), labels = c("Maize", "Small
grains", "Legumes",
"Other crops")) + xlab("Farms") + ylab("Area (ha)") + themeBW() +
scale_x continuous(expand = c(0,
0)) + scale_y continuous(expand = c(0, 0)) + theme(axial. text.x =
elementقد(),
axial. text.y = element_text(size = 12, face = "bold"), axial. title.x =
element_text(size = 14,
face = "bold", margin = margin(10, 0, 0, 0)), axial. title.y =
element_text(size = 14,
face = "bold", margin = margin(0, 10, 0, 0)), axial. ticks.x =
elementurre(),
legend. background = elementrect(color = "black"), legend. title =
elementurre(),
legend. text = element_text(size = 12), legend. position = c(0.9, 0.9),
legend. justification = c(0.9,
0.9))
```

<figure><img src="image_13.png" /></figure>

### 3.2.4 Lollipop diagram

We will plot the frequency of species used as intercrops.

We create a subset of the dataset, and rename variables to have nice labels for the plot.

```r
inter = data[, c(54:62, 284:287)]

names(inter) = c("Common bean", "Cowpea", "Groundnut", "Lablab",
  "Pigeonpea", "Velvet bean",
```
```
"Pumpkin", "Watermelon", "Local cucumber", "Maize", "Rosella", "Sweet
sorghum",
"Pearl millet")
```

We transform our dataset from wide to long format, sum the number of occurrences by intercrops, and remove zeros.

```r
to = gather(inter, intercrop, presence, "Common bean":"Pearl millet")
to = aggregate(presence ~ intercrop, FUN = sum, na.rm = TRUE, data =
  inter)
to = subset(inter, presence > 0)
```

Finally, we plot the lollipop diagram.

Note the use of `arrange()` and `coord flip()`.

```r
inter %>
  arrange(presence) %>
  mutate(order = factor(intercrop, intercrop)) %>
  ggplot(aes(x = intercrop, y = presence)) + geom_segment(aes(x =
  order, xend = intercrop,
  y = 0, yend = presence), color = "#ED6925FF", linewidth = 1) +
  geom_point(color = "#ED6925FF",
  size = 8, alpha = 0.5) + xlab("Intercrop species") + ylab("Nb of
  farms") + themeBW() +
  coord flip() + theme(axistitle.x = element_text(size = 14, face =
  "bold", margin = margin(10,
  0, 0, 0)), axis.title.y = element_text(size = 14, face = "bold",
  margin = margin(0,
  10, 0, 0)), axis.text.x = element_text(size = 12, face = "bold"),
  axis.text.y = element_text(size = 10,
  face = "bold"))
```
<figure><img src="image_14.png" /></figure>

### 3.2.5 Bubble chart

We will plot the food security status by month of the population of farms surveyed.

We create a subset of the dataset, rename variables to have nice labels for the plot, and remove empty rows.

```r
fsec = data[, c(212:223)]
```

```r
names(fsec) = c("Nov 2022", "Oct 2022", "Sep 2022", "Aug 2022", "Jul 2022", "Jun 2022", "May 2022", "Apr 2022", "Mar 2022", "Feb 2022", "Jan 2022", "Dec 2021")
```

```r
fsec = na.omit(fsec)
```

We transform the dataset from wide to long format, and reorder the levels of the factor 'status' and 'months' (the former for the data to be displayed from low to high food security, and the latter for data to be displayed in chronological order).

```r
fsec = gather(fsec, month, status, "Nov 2022":"Dec 2021")
```

```r
fsec$status = factor(fsec$status, levels = c("Low", "Medium", "High"))
```

```r
fsec$month = factor(fsec$month, levels = c("Dec 2021", "Jan 2022", "Feb 2022", "Mar 2022"))
```
<figure><img src="image_15.png" /></figure>

Finally, we plot the bubble chart.

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

### 3.2.6 Pie chart

We will plot the average use of maize residues.

We create a subset of the dataset, and rename variables to have nice labels for the plot.
maiz_res = data[, c(125:131)]

names(maiz_res) = c("Burnt", "Retained in the field", "Grazed in the
  → field", "Harvested for feed",
  "Harvested for compost", "Harvested for fuel", "Harvested for
  → fencing")

We transform the dataset from wide to long format, turn percentage values
(then characters) as numeric, calculate the mean percentage for each use, and
remove uses for which the percentage is zero.

maiz_res = gather(maiz_res, use, perc, "Burnt": "Harvested for fencing")

maiz_res$perc = as.numeric(sub("%", "", maiz_res$perc))

maiz_res = aggregate(perc ~ use, FUN = mean, na.rm = TRUE, data =
  maiz_res)

maiz_res = maiz_res[maiz_res$perc > 0, ]

Finally, we plot the pie chart.

Note the use of `coord_polar()` and `theme{void}`.

ggplot(maiz_res, aes(x = "", y = perc, fill = use)) + geom_bar(alpha =
  0.5, width = 1,
  stat = "identity", color = "black") + coord_polar("y", start = 0) +
  theme{'void'} +
  scale_fill_viridis_doption = "B", name = "Maize residue use") +
  theme(legend.text = element_text(size = 12),
  legend.title = element_text(size = 14, face = "bold"),
  legend(position = "right")
<figure><img src="image_18.png" /></figure>

Now that we have built a number of exploratory graphs, we will generate a summary table.

### 3.2.7 Summary table

Many R packages can be used to generate tables. Here, we will use the R package `gtsummary` (Sjoberg et al., 2021), in part for the nice rendering on RMarkdown (see below)

We will produce a summary table comparing the two wards from the small dataset created above, using the function `tbl_summary()`.

We display means and standard deviations for continuous variables, with one decimal only, and percentages for categorical variables.

For proper rendering as PDF, we convert the table as a kable object using the R package `kableExtra` (Zhu, 2024), and display it as landscape.

```r
small_dataset[, -c(1, 3)] %>

tbl_summary(by = ward, statistic = list(all continuous() ~ "{mean}")
           ({sd})", all categorial() ~ "{p}%")
           digits = all continuous() ~ 1) %>

add_overall() %>

as kable(format = "latex", booktabs = TRUE, longtable = TRUE) %>

kable styling(latex_options = c("scale_down", "repeat_header"),
  font_size = 6) %>

landscape()
```
<table><thead><tr><th>**Characteristic**</th><th>**Overall** N = 266</th><th>**Ward 2** N = 132</th><th>**Ward 3 9 12 17** N = 134</th></tr></thead><tbody><tr><td>how_old_is_the_head_of_theHousehold</td><td>51.4 (15.7)</td><td>49.7 (16.0)</td><td>53.1 (15.3)</td></tr><tr><td>what_is_the sex_of_the_head_of_theHousehold</td><td></td><td></td><td></td></tr><tr><td>Female</td><td>23%</td><td>21%</td><td>25%</td></tr><tr><td>Male</td><td>77%</td><td>79%</td><td>75%</td></tr><tr><td>what_is_the_education_level_of_the_head_of_theHousehold</td><td></td><td></td><td></td></tr><tr><td>Primary level</td><td>60%</td><td>55%</td><td>65%</td></tr><tr><td>Secondary level</td><td>36%</td><td>40%</td><td>32%</td></tr><tr><td>Tertiary level</td><td>1.1%</td><td>0%</td><td>2.2%</td></tr><tr><td>Vocational school</td><td>2.6%</td><td>4.5%</td><td>0.7%</td></tr><tr><td>family_size</td><td>5.6 (3.2)</td><td>6.0 (3.6)</td><td>5.3 (2.6)</td></tr><tr><td>cropped_area</td><td>1.8 (1.5)</td><td>1.5 (1.3)</td><td>2.1 (1.6)</td></tr><tr><td>fallow</td><td>0.7 (1.4)</td><td>0.6 (1.0)</td><td>0.9 (1.6)</td></tr><tr><td>prop_non_cereals</td><td>29.9 (29.1)</td><td>30.0 (30.8)</td><td>29.8 (27.5)</td></tr><tr><td>fertilizers</td><td>87.0 (115.2)</td><td>72.9 (90.9)</td><td>101.0 (133.9)</td></tr><tr><td>amendment</td><td>304.5 (626.1)</td><td>294.1 (653.2)</td><td>314.8 (600.5)</td></tr><tr><td>cattle</td><td>4.3 (8.9)</td><td>3.5 (8.7)</td><td>5.2 (9.1)</td></tr><tr><td>small_rum</td><td>10.4 (20.1)</td><td>9.4 (17.4)</td><td>11.4 (22.5)</td></tr><tr><td>poultry</td><td>11.7 (29.9)</td><td>10.0 (11.3)</td><td>13.3 (40.6)</td></tr><tr><td>bee_hives</td><td>11%</td><td>7.6%</td><td>13%</td></tr><tr><td>ind_garden</td><td>48%</td><td>55%</td><td>40%</td></tr><tr><td>com_garden</td><td>2.3%</td><td>3.8%</td><td>0.7%</td></tr><tr><td>main_food_source</td><td></td><td></td><td></td></tr><tr><td>Cash purchase from income</td><td>3.4%</td><td>3.8%</td><td>3.0%</td></tr><tr><td>Casual labour for food</td><td>18%</td><td>30%</td><td>6.0%</td></tr><tr><td>Food aid</td><td>3.0%</td><td>5.3%</td><td>0.7%</td></tr><tr><td>Own production</td><td>75%</td><td>61%</td><td>90%</td></tr><tr><td>Purchase from cash transfer</td><td>0.4%</td><td>0%</td><td>0.7%</td></tr><tr><td>Remittances</td><td>0.4%</td><td>0.8%</td><td>0%</td></tr><tr><td>main_income_source</td><td></td><td></td><td></td></tr><tr><td>Artisan</td><td>0.8%</td><td>1.5%</td><td>0%</td></tr><tr><td>Casual labour</td><td>35%</td><td>43%</td><td>28%</td></tr><tr><td>Crop sales</td><td>36%</td><td>28%</td><td>44%</td></tr><tr><td>Livestock sales</td><td>15%</td><td>9.8%</td><td>20%</td></tr><tr><td>Pension</td><td>0.4%</td><td>0.8%</td><td>0%</td></tr><tr><td>Remittances</td><td>5.3%</td><td>7.6%</td><td>3.0%</td></tr><tr><td>Salary or wages</td><td>3.4%</td><td>5.3%</td><td>1.5%</td></tr></tbody></table>
(continued)

<table><thead><tr><th>**Characteristic**</th><th>**Overall** N = 266</th><th>**Ward 2** N = 132</th><th>**Ward 3** 9 12 17** N = 134</th></tr></thead><tbody><tr><td>Smallscale mining</td><td>0.4%</td><td>0%</td><td>0.7%</td></tr><tr><td>Trade</td><td>3.4%</td><td>3.8%</td><td>3.0%</td></tr><tr><td>food_security</td><td>0.6 (0.3)</td><td>0.5 (0.3)</td><td>0.6 (0.3)</td></tr><tr><td>hdds</td><td></td><td></td><td></td></tr><tr><td>1</td><td>1.1%</td><td>2.3%</td><td>0%</td></tr><tr><td>2</td><td>30%</td><td>12%</td><td>47%</td></tr><tr><td>3</td><td>36%</td><td>47%</td><td>25%</td></tr><tr><td>4</td><td>23%</td><td>28%</td><td>18%</td></tr><tr><td>5</td><td>8.6%</td><td>9.8%</td><td>7.5%</td></tr><tr><td>6</td><td>1.5%</td><td>0.8%</td><td>2.2%</td></tr><tr><td>8</td><td>0.4%</td><td>0%</td><td>0.7%</td></tr><tr><td>community_seed_banks</td><td>12%</td><td>1.5%</td><td>22%</td></tr><tr><td>small_grains</td><td>87%</td><td>86%</td><td>87%</td></tr><tr><td>crop_rotate</td><td>51%</td><td>49%</td><td>53%</td></tr><tr><td>intercropping</td><td>17%</td><td>11%</td><td>23%</td></tr><tr><td>cover_crops_for_erosion_control_and_or_soil_fertility</td><td>8.6%</td><td>6.8%</td><td>10%</td></tr><tr><td>mulching</td><td>9.8%</td><td>9.8%</td><td>9.7%</td></tr><tr><td>integrated_pest_management_i_e_scouting_etc</td><td>50%</td><td>43%</td><td>57%</td></tr><tr><td>compost_manure</td><td>17%</td><td>17%</td><td>17%</td></tr><tr><td>tree_planting_on_farm</td><td>14%</td><td>9.1%</td><td>19%</td></tr><tr><td>treeRetention_natural_regeneration_on_farm</td><td>37%</td><td>44%</td><td>30%</td></tr><tr><td>homemade_animalfeeds_made_with_localy_available_ingredients</td><td>26%</td><td>28%</td><td>23%</td></tr><tr><td>fodder_production_for_ruminants_e_g_velvet_bean_lablab</td><td>4.5%</td><td>4.5%</td><td>4.5%</td></tr><tr><td>fodder_preservation_for_ruminants_e_g_silage_making</td><td>4.1%</td><td>1.5%</td><td>6.7%</td></tr><tr><td>survivaleron_feeding_of_productive_livestock_in</td><td>4.5%</td><td>4.5%</td><td>4.5%</td></tr><tr><td>tot_div</td><td>6.8 (3.7)</td><td>6.7 (3.4)</td><td>7.0 (3.9)</td></tr><tr><td>cerealprod</td><td>867.0 (962.0)</td><td>548.5 (464.6)</td><td>1,180.8 (1,196.4)</td></tr><tr><td>offtake</td><td>0.4 (1.3)</td><td>0.3 (1.0)</td><td>0.4 (1.5)</td></tr><tr><td>eq_value</td><td>290.0 (689.5)</td><td>206.4 (308.1)</td><td>372.4 (916.5)</td></tr></tbody></table>
We add a line 'type' to force hdds to be read as a continuous variable.

We add 'label' lines to recode variable labels into nicer labels.

We add a column with a statistical test comparing values between the two
wards - through the function *add_p()* - and a column with values for the overall
population - with the function *add_overall*().

small_dataset[, -c(1, 3)] %>
  tbl_summary(by = ward, statistic = list(all continuous() ~ "{mean}")
  ({sd}), all categorica l() ~
  "{p}%"), type = list(hdds ~ "continuous"), digits =
  all continuous() ~ 1,
  label = list(how_old_is the head of the household ~ "Age of the
  head of the household",
  what is the sex of the head of the household ~ "Female-headed
  households",
  what is the education level of the head of the household ~
  "Education of the head of the household higher than primary",
  what is the education level of the head of the household ~
  "Family size",
  eq_value ~ "Equipment value (USD)", cropped area ~ "Total
  cropped area (ha)",
  fallow ~ "Fallow land (ha)", prop_non cereals ~ "Non-cereal
  crops (% total cropped area)",
  fertilizers ~ "Total fertilizer used (kg)", amendment ~
  "Total manure + compost used (kg)",
  cattle ~ "Cattle (n)", small rum ~ "Small ruminants (n)",
  poultry ~ "Poultry (n)",
  bee hives ~ "Bee keeping", ind garden ~ "Owning an individual
  garden",
  com garden ~ "Having access to a communal garden",
  main food source ~
  "Own production as main source of food",
  main income source ~ "Farming as main source of income",
  food security ~ "Proportion of the year being food secured",
  hdds ~ "24H household dietary diversity score (0-12)",
  community seed banks ~ "Using a community seed bank",
  small grains ~
  "Using small grains", crop rotation ~ "Using crop
  rotation", intercropping ~
  "Using intercropping",
  cover crops for erosion control and or soil fertility ~
  "Using cover crops", mulching ~ "Using mulching",
  integrated pest management i e scouting etc ~
  "Using integrated pest management", compost manure ~
  "Using compost and manure",
  tree planting on farm ~ "Planting trees on-farm",
  tree retention natural regeneration on farm ~
  "Retaining naturally regenerated trees on-farm",
  homemade animal feeds made with locally available ingredients ~
```
"Using homemade animal feed",
-> fodder_production_for_ruminants_e_g_velvet_bean_lablab ~
    "Produccing fodder",
-> fodder_preservation_for_ruminants_e_g_silage_making ~
    "Preserving fodder",
-> survival_offering_offering_of_productive_livestock_inUsing survival offering", tot_div ~ "Total farm diversity
(n species)",
    cerealProd ~ "Total cereal production (kg/yr)", offtake ~
-> "Total livestock offtake (TLU/yr)") %>
add_p(test = list(all continuous() ~ "kruskal.test",
-> all categorial() ~ "chisq.test")) %>
add_overall() %>
as kable(format = "latex", booktabs = TRUE, longtable = TRUE) %>
kable_styling(latex_options = c("scale_down", "repeat_header"),
-> font_size = 6) %>
landscape()
```
<table><thead><tr><th>**Characteristic**</th><th>**Overall** N = 266</th><th>**Ward 2** N = 132</th><th>**Ward 3 9 12 17** N = 134</th><th>**p-value**</th></tr></thead><tbody><tr><td>Age of the head of the household</td><td>51.4 (15.7)</td><td>49.7 (16.0)</td><td>53.1 (15.3)</td><td>0.036</td></tr><tr><td>Female-headed households</td><td></td><td></td><td></td><td>0.6</td></tr><tr><td>Female</td><td>23%</td><td>21%</td><td>25%</td><td></td></tr><tr><td>Male</td><td>77%</td><td>79%</td><td>75%</td><td></td></tr><tr><td>Family size</td><td></td><td></td><td></td><td>0.032</td></tr><tr><td>Primary level</td><td>60%</td><td>55%</td><td>65%</td><td></td></tr><tr><td>Secondary level</td><td>36%</td><td>40%</td><td>32%</td><td></td></tr><tr><td>Tertiary level</td><td>1.1%</td><td>0%</td><td>2.2%</td><td></td></tr><tr><td>Vocational school</td><td>2.6%</td><td>4.5%</td><td>0.7%</td><td></td></tr><tr><td>family_size</td><td>5.6 (3.2)</td><td>6.0 (3.6)</td><td>5.3 (2.6)</td><td>0.10</td></tr><tr><td>Total cropped area (ha)</td><td>1.8 (1.5)</td><td>1.5 (1.3)</td><td>2.1 (1.6)</td><td>&lt;0.001</td></tr><tr><td>Fallow land (ha)</td><td>0.7 (1.4)</td><td>0.6 (1.0)</td><td>0.9 (1.6)</td><td>0.2</td></tr><tr><td>Non-cereal crops (% total cropped area)</td><td>29.9 (29.1)</td><td>30.0 (30.8)</td><td>29.8 (27.5)</td><td>0.7</td></tr><tr><td>Total fertilizer used (kg)</td><td>87.0 (115.2)</td><td>72.9 (90.9)</td><td>101.0 (133.9)</td><td>0.14</td></tr><tr><td>Total manure + compost used (kg)</td><td>304.5 (626.1)</td><td>294.1 (653.2)</td><td>314.8 (600.5)</td><td>0.6</td></tr><tr><td>Cattle (n)</td><td>4.3 (8.9)</td><td>3.5 (8.7)</td><td>5.2 (9.1)</td><td>&lt;0.001</td></tr><tr><td>Small ruminants (n)</td><td>10.4 (20.1)</td><td>9.4 (17.4)</td><td>11.4 (22.5)</td><td>0.7</td></tr><tr><td>Poultry (n)</td><td>11.7 (29.9)</td><td>10.0 (11.3)</td><td>13.3 (40.6)</td><td>0.4</td></tr><tr><td>Bee keeping</td><td>11%</td><td>7.6%</td><td>13%</td><td>0.2</td></tr><tr><td>Owning an individual garden</td><td>48%</td><td>55%</td><td>40%</td><td>0.020</td></tr><tr><td>Having access to a communal garden</td><td>2.3%</td><td>3.8%</td><td>0.7%</td><td>0.2</td></tr><tr><td>Own production as main source of food</td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>Cash purchase from income</td><td>3.4%</td><td>3.8%</td><td>3.0%</td><td></td></tr><tr><td>Casual labour for food</td><td>18%</td><td>30%</td><td>6.0%</td><td></td></tr><tr><td>Food aid</td><td>3.0%</td><td>5.3%</td><td>0.7%</td><td></td></tr><tr><td>Own production</td><td>75%</td><td>61%</td><td>90%</td><td></td></tr><tr><td>Purchase from cash transfer</td><td>0.4%</td><td>0%</td><td>0.7%</td><td></td></tr><tr><td>Remittances</td><td>0.4%</td><td>0.8%</td><td>0%</td><td></td></tr><tr><td>Farming as main source of income</td><td></td><td></td><td></td><td>0.003</td></tr><tr><td>Artisan</td><td>0.8%</td><td>1.5%</td><td>0%</td><td></td></tr><tr><td>Casual labour</td><td>35%</td><td>43%</td><td>28%</td><td></td></tr><tr><td>Crop sales</td><td>36%</td><td>28%</td><td>44%</td><td></td></tr><tr><td>Livestock sales</td><td>15%</td><td>9.8%</td><td>20%</td><td></td></tr><tr><td>Pension</td><td>0.4%</td><td>0.8%</td><td>0%</td><td></td></tr><tr><td>Remittances</td><td>5.3%</td><td>7.6%</td><td>3.0%</td><td></td></tr><tr><td>Salary or wages</td><td>3.4%</td><td>5.3%</td><td>1.5%</td><td></td></tr></tbody></table>
(continued)

<table><thead><tr><th>**Characteristic**</th><th>**Overall** N = 266</th><th>**Ward 2** N = 132</th><th>**Ward 3 9 12 17** N = 134</th><th>**p-value**</th></tr></thead><tbody><tr><td>Smallscale mining</td><td>0.4%</td><td>0%</td><td>0.7%</td><td></td></tr><tr><td>Trade</td><td>3.4%</td><td>3.8%</td><td>3.0%</td><td></td></tr><tr><td>Proportion of the year being food secured</td><td>0.6 (0.3)</td><td>0.5 (0.3)</td><td>0.6 (0.3)</td><td>0.022</td></tr><tr><td>24H household dietary diversity score (0-12)</td><td>3.1 (1.1)</td><td>3.3 (0.9)</td><td>3.0 (1.2)</td><td>&lt;0.001</td></tr><tr><td>Using a community seed bank</td><td>12%</td><td>1.5%</td><td>22%</td><td>&lt;0.001</td></tr><tr><td>Using small grains</td><td>87%</td><td>86%</td><td>87%</td><td>&gt;0.9</td></tr><tr><td>Using crop rotation</td><td>51%</td><td>49%</td><td>53%</td><td>0.6</td></tr><tr><td>Using intercropping</td><td>17%</td><td>11%</td><td>23%</td><td>0.018</td></tr><tr><td>Using cover crops</td><td>8.6%</td><td>6.8%</td><td>10%</td><td>0.4</td></tr><tr><td>Using mulching</td><td>9.8%</td><td>9.8%</td><td>9.7%</td><td>&gt;0.9</td></tr><tr><td>Using integrated pest management</td><td>50%</td><td>43%</td><td>57%</td><td>0.027</td></tr><tr><td>Using compost and manure</td><td>17%</td><td>17%</td><td>17%</td><td>&gt;0.9</td></tr><tr><td>Planting trees on-farm</td><td>14%</td><td>9.1%</td><td>19%</td><td>0.026</td></tr><tr><td>Retaining naturally regenerated trees on-farm</td><td>37%</td><td>44%</td><td>30%</td><td>0.024</td></tr><tr><td>Using homemade animal feed</td><td>26%</td><td>28%</td><td>23%</td><td>0.4</td></tr><tr><td>Producing fodder</td><td>4.5%</td><td>4.5%</td><td>4.5%</td><td>&gt;0.9</td></tr><tr><td>Preserving fodder</td><td>4.1%</td><td>1.5%</td><td>6.7%</td><td>0.068</td></tr><tr><td>Using survival feeding</td><td>4.5%</td><td>4.5%</td><td>4.5%</td><td>&gt;0.9</td></tr><tr><td>Total farm diversity (n species)</td><td>6.8 (3.7)</td><td>6.7 (3.4)</td><td>7.0 (3.9)</td><td>0.7</td></tr><tr><td>Total cereal production (kg/yr)</td><td>867.0 (962.0)</td><td>548.5 (464.6)</td><td>1,180.8 (1,196.4)</td><td>&lt;0.001</td></tr><tr><td>Total livestock offtake (TLU/yr)</td><td>0.4 (1.3)</td><td>0.3 (1.0)</td><td>0.4 (1.5)</td><td>0.6</td></tr><tr><td>Equipment value (USD)</td><td>290.0 (689.5)</td><td>206.4 (308.1)</td><td>372.4 (916.5)</td><td>0.004</td></tr></tbody></table>
3.3 RMarkdown

RMarkdown is an authoring framework supported by the R package rmark-
down (Allaire et al., 2024) that combines text, code, and outputs (e.g., tables,
plots, and results) in a single document. It executes codes (R, Python, etc) and
generates high quality reports (pdf, word, html, etc).

You can develop a RMarkdown baseline report while your survey is ongoing:
tables, plots and results will update as the dataset gets populated. And you
can produce your final baseline report and share it with your partners and
donors the same day the last farm is being interviewed!

Here are the codes to generate a simple RMarkdown baseline report with the fig-
ures we generated earlier: https://github.com/FBaudron/Supporting-codesign/
blob/main/03-baselining/simple_rmarkdown_baseline_report.Rmd.

And here is the pdf rendered by these codes: https://github.com/FBaudron/Supporting-codesign/blob/main/03-baselining/simple_rmarkdown_baseline_report.pdf.
# Chapter 4

## Farm typologies

This chapter proposes a step-by-step method to understand farm diversity and select representative farms for e.g., on-farm experiments.

### 4.1 Concept of farm typology

Farm typologies aim at characterizing the diversity of farming systems and their distribution in heterogenous communities, as a basis for prioritizing interventions (Hammond et al., 2020) and/or to select representative farms, or ‘prototype farms’ (Alvarez et al., 2018).

Various methodologies, all with their strengths and their weaknesses, are commonly used to delineate typologies, from (quantitative) statistical typology (Tittonell et al., 2010) to (qualitative) expert-based typology involving farmers and other knowledgeable stakeholders (Kebede, 2009).

Here, we are using a method of statistical typology delineation developed by Hassall et al. (2023). One of the key strength of the method is that non-metric multidimensional scaling can be used, allowing for the selection of qualitative variables, not only quantitative ones, for typology delineation.

In the example below, we will use the small dataset we created in the previous chapter (and that was saved as a rds file). You can also download it here (and load it with the function *read.csv*()): https://github.com/FBaudron/Supporting-codesign/blob/main/04-typologies/small_dataset.csv

A standalone R script containing the full workflow is available at the following link: https://github.com/FBaudron/Supporting-codesign/blob/main/04-typologies/typology.R.
## 4.2 Data preparation

We begin by clearing the environment, ensuring all necessary packages are installed and loaded, and importing the data.

```r
rm(list = ls())
```

```r
if (!require("openxlsx")) install.packages("openxlsx")
if (!require("vegan")) install.packages("vegan")
if (!require("ade4")) install.packages("ade4")
if (!require("dendextend")) install.packages("dendextend")
if (!require("gtsummary")) install.packages("gtsummary")
if (!require("kableExtra")) install.packages("kableExtra")
if (!require("splitstackshape")) install.packages("splitstackshape")
if (!require("NbClust")) install.packages("NbClust")
if (!require("randomForest")) install.packages("randomForest")
```

```r
library(openxlsx)
library(vegan)
library(ade4)
library(dendextend)
library gtsummary
library(kableExtra)
library splitstackshape
library(NbClust)
library(randomForest)
```

```r
data = readRDS("small_dataset.rds")
```

```r
# data = read.xlsx('small_dataset.xlsx', sheet = 1)
```

We first rename some variables, recode some categorical variables into binary ones, and convert NAs into zeros for some variables.

```r
names(data)[c(1, 4:6, 27, 29, 33:36)] = c("ID", "age_hhh", "sex_hhh",
  "edu_hhh",
  "cover_crops", "integrated_pest_management", "homemadefeeds",
  "fodder_production",
  "fodder_preservation", "survival_offering")
```

```r
data$sex_hhh = as.numeric(data$sex_hhhh == "Female")
```

```r
data$edu_hhhh = as.numeric(data$edu_hhhh != "Primary level")
```

```r
data$main_food_source = as.numeric(data$main_food_source == "Own production")
```
```r
data$main_income_source = as.numeric(data$main_income_source == "Crop sales" | data$main_income_source == "Livestock sales")

data[c(23:36)][is.na(data[c(23:36)])] = 0
```

### 4.2.1 Removing outliers

We remove extreme values to avoid skewing the analysis. For this, we plot and examine boxplots. Points beyond the whiskers of the boxplot have values higher than 1.5 times the interquartile range and are considered potential outliers. When such extreme points are observable, we remove them from the dataset.

Age of the head of the household doesn't have any true marked outlier.

```r
boxplot(data$age_hhh)
```

<figure><img src="image_19.png" /></figure>

But family size does.

```r
boxplot(data$family_size)
```
<figure><img src="image_20.png" /></figure>

Thus, we remove all records with a family size greater than 20.

```r
data = data[$family_size < 20, ]
boxplot(data$family_size)
```

<figure><img src="image_21.png" /></figure>

We inspect all other continuous variables, and find outliers for the number of small ruminants, the number of poultry, livestock offtake, and equipment value. We remove these extreme values from the dataset that we will use to delineate the farm typology.
```r
# boxplot(data$cropped_area)
# boxplot(data$fallow)
# boxplot(data$prop_non_cereals)
# boxplot(data$fertilizers)
# boxplot(data$amendment)
# boxplot(data$cattle)
# boxplot(data$small_rum)
data = data(data$small_rum < 100, ]
# boxplot(data$small_rum)
# boxplot(data$poultry)
data = data(data$poultry < 400, ]
# boxplot(data$poultry)
# boxplot(data$food_security)
# boxplot(data$hdds)
# boxplot(data$tot_div)
# boxplot(data$cereal ProD)
# boxplot(data$offtake)
data = data(data$offtake < 6, ]
# boxplot(data$offtake)
# boxplot(data$eq_value)
data = data(data$eq_value < 8000, ]
# boxplot(data$eq_value)
```

### 4.2.2 Normalization of continuous variables

To reduce the influence of extreme values, continuous variables are log-transformed.

We create first a copy of the dataset before data transformation.

```
typo = data
```

We plot histograms of continuous variables, and when distribution departs too much from a Gaussian distribution, we logtransform the variable. In that case,
we take care of adding a term - half of the minimum positive value - to avoid
passing zeros to log10.

Age of the head of the household has a relatively Gaussian distribution.

hist(typo$age_hhh)

### Histogram of typo$age_hhh

<figure><img src="image_22.png" /></figure>

But family size has a distribution that departs markedly from a Gaussian distribution.

hist(typo$family_size)

### Histogram of typo$family_size

<figure><img src="image_23.png" /></figure>
Thus, we logtransform the variable family size.

```r
typo$family_size = log10(typo$family_size + (0.5 *
  min(typo$family_size[typo$family_size >
    0], na.rm = TRUE)))
hist(typo$family_size)
```

<figure><img src="image_24.png" /><figcaption>Histogram of typo$family_size</figcaption></figure>

We inspect all other continuous variables, and find the total cropped area, the fallow area, the total quantity of fertilizer used, the total quantity of organic amendments used, the number of cattle, the number of small ruminants, the number of poultry, the total quantity of cereal produced, the livestock offtake, and the equipment value to depart markedly from a Gaussian distribution. We therefore log-transformed these variables.

```r
# hist(typo$cropped_area)
typo$cropped_area = log10(typo$cropped_area + (0.5 *
  min(typo$cropped_area[typo$cropped_area >
    0], na.rm = TRUE)))
# hist(typo$cropped_area)

# hist(typo$fallow)
typo$fallow = log10(typo$fallow + (0.5 * min(typo$fallow[typo$fallow >
  0], na.rm = TRUE)))
# hist(typo$fallow)

# hist(typo$prop_non_cereals)

# hist(typo$fertilizers)
typo$fertilizers = log10(typo$fertilizers + (0.5 *
  min(typo$fertilizers[typo$fertilizers >
    0], na.rm = TRUE)))
```
```
0], na.rm = TRUE)))
# hist(typo$fertilizers)

# hist(typo$amendment)
typo$amendment = log10(typo$amendment + (0.5 *
  min(typo$amendment[typo$amendment >
  0], na.rm = TRUE)))
# hist(typo$amendment)

# hist(typo$cattle)
typo$cattle = log10(typo$cattle + (0.5 * min(typo$cattle[typo$cattle >
  0], na.rm = TRUE)))
# hist(typo$cattle)

# hist(typo$small_rum)
typo$small_rum = log10(typo$small_rum + (0.5 *
  min(typo$small_rum[typo$small_rum >
  0], na.rm = TRUE)))
# hist(typo$small_rum)

# hist(typo$poultry)
typo$poultry = log10(typo$poultry + (0.5 * min(typo$poultry[typo$poultry
  > 0], na.rm = TRUE)))
# hist(typo$poultry)

# hist(typo$food_security)

# hist(typo$hdds)

# hist(typo$tot_div)

# hist(typo$cereal ProD)
typo$cereal ProD = log10(typo$cereal ProD + (0.5 *
  min(typo$cereal ProD[typo$cereal ProD >
  0], na.rm = TRUE)))
# hist(typo$cereal ProD)

# hist(typo$offtake)
typo$offtake = log10(typo$offtake + (0.5 * min(typo$offtake[typo$offtake
  > 0], na.rm = TRUE)))
# hist(typo$offtake)

# hist(typo$eq_value)
typo$eq_value = log10(typo$eq_value + (0.5 *
  min(typo$eq_value[typo$eq_value > 0],
  na.rm = TRUE)))
# hist(typo$eq_value)
```
## 4.3 Delineation of farm types

Next, we calculate dissimilarities and distances, reduce the dimensionality of the dataset through principal coordinate analysis, and delineate farm types through hierarchical cluster analysis.

We first remove records with missing values.

```
typo = na.omit(typo)
```

### 4.3.1 Calculation of dissimilarities and distances

We calculate dissimilarities for 3 groups of variables: structural variables, functional variables, and adoption variables. Structural and functional variables are further divided into continuous and binary variables, for which different methods are used in the calculation of dissimilarities.

We first calculate the scaled euclidean dissimilarity between continuous structural variables using the function *vegdist()* of the R package **vegan** (Oksanen et al., 2025).

```
dEucs = vegdist(typo[, c("age_hhh", "family_size", "cropped_area",
  "fallow", "prop_non_cereals",
  "cattle", "small_rum", "poultry", "eq_value")], method = "gower")
```

Similarly, we calculate the scaled euclidean dissimilarity between continuous functional variables.

```
dEucf = vegdist(typo[, c("fertilizers", "amendment", "food_security",
  "hdds", "tot_div",
  "cerealProd", "offtake")], method = "gower")
```

We then calculate the dissimilarity between binary structural variables using the function *dist.binary()* of the R package **ade4** (Thioulouse et al., 2018).

```
dBin_ppls = dist.binary(typo[, c("sex_hhh", "edu_hhh", "bee_hives",
  "ind_garden",
  "com_garden")], method = 2)
```

Similarly, we calculate the dissimilarity between binary functional variables.

```
dBin_pplf = dist.binary(typo[, c("main_food_source",
  "main_income_source")], method = 2)
```
Finally, we calculate the dissimilarity between binary adoption variables.

```
d_practices = dist-binary(typo[, c("community_seed_banks",
  "small_grains", "crop, "intercropping", "cover_crops", "mulching",
  "integrated_pest_management", "compost_manure",
  "tree_planting_on_farm",
  "tree_retention_natural_regeneration_on_farm", "homemade feeds",
  "fodder_production", "fodder_preservation", "survival feeding")],
  method = 2)
```

We combine the 2 dissimilarities of structural variables, weighing them by the number of variables in each dissimilarity.

We do the same for the 2 dissimilarities of functional variables.

```
d Structural = (9 * dEucs^2 + 5 * dBin_ppls^2)/14
```

```
d Functional = (7 * dEucf^2 + 2 * dBin_pplf^2)/9
```

We combine the three dissimilarity matrices — structural, functional, and adoption variables — into a single overall dissimilarity. Because the adoption component includes many binary variables, we downweight its contribution to avoid letting it dominate the resulting typology. Specifically, we multiply the adoption dissimilarity by 0.5 before combining, and then normalize by dividing by the sum of the weights (2.5) to keep the overall scale consistent. The factor of 0.5 used is arbitrary, but could be fine-tuned through iterations.

```
dAll = (dStructural + dFunctional + 0.5 * dPractices^2)/2.5
```

Finally, we transform the matrix of dissimilarities to a matrix of distances.

```
distAll = sqrt(2 * dAll)
```

### 4.3.2 Principal coordinate analysis & hierarchical cluster analysis

We reduce dimensionality while preserving distances using using principal coordinate analysis (PCA) implemented via the base R function `cmdscale()`.

We select k axes to explain at least 50% of the variability.
```
pco = cmdscale(distAll, eig = TRUE, k = 10)
barplot(pco$eig[1:20])
```

<figure><img src="image_25.png" /></figure>

```
cum_var = cumsum(pco$eig)/sum(pco$eig)
k_axes = which(cum_var >= 0.5)[1]
k_axes
```

[1] 4

We cluster farms based on PCO axes through hierarchical cluster analysis (HCA) carried out using the base R function *hclust()*.

We plot the dendrogram.

```
pco_var = pco$points[, 1:k_axes]

hc_pco = hclust(dist(pco_var), method = "complete")
plot(hc_pco, hang = -1)
```
<figure><img src="image_26.png" /></figure>

4.3.3 Selection of the number of clusters

We select the optimum number of clusters using the function *NbClust()* of the R package **NbClust** (Charrad et al., 2014), which computes multiple cluster validity indices (30 in our case) and compares them. The final number of clusters is determined using a majority rule, where each index 'votes' for a number of clusters and the number supported by the largest number of indices is chosen.

NbClust(pco_var, diss = dist(pco_var), distance = NULL, method =
  "complete")
<figure><img src="image_27.png" /></figure>

\*\*\* : The Hubert index is a graphical method of determining the number of clusters. In the plot of Hubert index, we seek a significant knee that corresponds to a significant increase of the value of the measure i.e the significant point in the second differences plot.

<figure><img src="image_28.png" /></figure>

\*\*\* : The D index is a graphical method of determining the number of clusters. In the plot of D index, we seek a significant knee (the significant point in the second differences plot) that corresponds to a significant increase of the measure.
* Among all indices:

* 4 proposed 2 as the best number of clusters

* 2 proposed 3 as the best number of clusters

* 5 proposed 4 as the best number of clusters

* 1 proposed 6 as the best number of clusters

* 2 proposed 7 as the best number of clusters

* 4 proposed 10 as the best number of clusters

* 1 proposed 13 as the best number of clusters

* 4 proposed 15 as the best number of clusters

***** Conclusion *****

* According to the majority rule, the best number of clusters is 4

$All.index

<table><thead><tr><th></th><th>KL</th><th>CH</th><th>Hartigan</th><th>CCC</th><th>Scott</th><th>Marriot</th><th>TrCovW</th><th>TraceW</th></tr></thead><tbody><tr><td>2</td><td>1.9625</td><td>72.0646</td><td>43.1709</td><td>-5.1974</td><td>236.6688</td><td>1596.8350</td><td>30.0213</td><td>19.0591</td></tr><tr><td>3</td><td>0.4950</td><td>63.4680</td><td>68.8217</td><td>-9.0010</td><td>353.6366</td><td>2279.2036</td><td>18.3870</td><td>16.2997</td></tr><tr><td>4</td><td>7.4435</td><td>76.4151</td><td>33.7409</td><td>-5.7550</td><td>633.6448</td><td>1362.9670</td><td>10.2851</td><td>12.8248</td></tr><tr><td>5</td><td>0.1917</td><td>73.0997</td><td>45.3120</td><td>-6.1273</td><td>761.7904</td><td>1293.4722</td><td>8.3896</td><td>11.3157</td></tr><tr><td>6</td><td>1.3434</td><td>77.7476</td><td>36.5987</td><td>-3.1357</td><td>944.7205</td><td>914.0985</td><td>5.8329</td><td>9.5911</td></tr><tr><td>7</td><td>8.0015</td><td>80.0165</td><td>15.2453</td><td>-1.1029</td><td>1053.3798</td><td>815.2061</td><td>4.7069</td><td>8.3706</td></tr><tr><td>8</td><td>0.1962</td><td>74.6461</td><td>25.6635</td><td>-1.6139</td><td>1100.9807</td><td>884.7331</td><td>4.0523</td><td>7.8895</td></tr><tr><td>9</td><td>1.0510</td><td>74.9528</td><td>23.7527</td><td>-0.5060</td><td>1206.6020</td><td>742.3897</td><td>3.5915</td><td>7.1523</td></tr><tr><td>10</td><td>34.9270</td><td>75.3400</td><td>12.3494</td><td>0.5180</td><td>1302.6660</td><td>630.6844</td><td>3.1240</td><td>6.5272</td></tr><tr><td>11</td><td>0.6694</td><td>72.1379</td><td>11.3744</td><td>0.2748</td><td>1346.8089</td><td>642.6907</td><td>2.7626</td><td>6.2163</td></tr><tr><td>12</td><td>0.1881</td><td>69.3630</td><td>7.7251</td><td>0.0600</td><td>1392.9727</td><td>639.0998</td><td>2.5024</td><td>5.9416</td></tr><tr><td>13</td><td>35.2137</td><td>65.9610</td><td>9.2750</td><td>-0.4788</td><td>1427.3519</td><td>656.1405</td><td>2.3591</td><td>5.7600</td></tr><tr><td>14</td><td>0.0045</td><td>63.6531</td><td>23.5850</td><td>-0.7189</td><td>1473.8089</td><td>635.1265</td><td>2.2128</td><td>5.5491</td></tr><tr><td>15</td><td>8.1803</td><td>66.2540</td><td>9.7789</td><td>0.7913</td><td>1589.2913</td><td>465.1975</td><td>1.8209</td><td>5.0581</td></tr></tbody></table><table><thead><tr><th></th><th>Friedman</th><th>Rubin</th><th>Cindex</th><th>DB Silhouette</th><th>Duda</th><th>Pseudot2</th><th>Beale</th><th>Ratkowsky</th></tr></thead><tbody><tr><td>2</td><td>1.1437</td><td>1.2826</td><td>0.4627</td><td>1.8797</td><td>0.2032</td><td>0.7318</td><td>46.1796</td><td>0.8779</td><td>0.1985</td></tr><tr><td>3</td><td>1.9536</td><td>1.4997</td><td>0.4461</td><td>1.6481</td><td>0.1752</td><td>0.6038</td><td>83.3362</td><td>1.5718</td><td>0.2320</td></tr><tr><td>4</td><td>4.5921</td><td>1.9061</td><td>0.4219</td><td>1.5696</td><td>0.2122</td><td>0.5819</td><td>43.8218</td><td>1.7064</td><td>0.2608</td></tr><tr><td>5</td><td>5.3597</td><td>2.1603</td><td>0.4035</td><td>1.3629</td><td>0.2229</td><td>0.5600</td><td>49.5005</td><td>1.8673</td><td>0.2900</td></tr><tr><td>6</td><td>7.3190</td><td>2.5488</td><td>0.4416</td><td>1.3070</td><td>0.2533</td><td>0.6521</td><td>35.7503</td><td>1.2692</td><td>0.2938</td></tr><tr><td>7</td><td>8.2513</td><td>2.9204</td><td>0.4134</td><td>1.3268</td><td>0.2640</td><td>0.6361</td><td>21.1712</td><td>1.3450</td><td>0.2906</td></tr><tr><td>8</td><td>8.8637</td><td>3.0985</td><td>0.4172</td><td>1.2797</td><td>0.2509</td><td>0.5048</td><td>33.3575</td><td>2.3009</td><td>0.2762</td></tr><tr><td>9</td><td>10.2412</td><td>3.4178</td><td>0.4039</td><td>1.1882</td><td>0.2611</td><td>0.6503</td><td>31.1912</td><td>1.2763</td><td>0.2686</td></tr><tr><td>10</td><td>11.4810</td><td>3.7452</td><td>0.3929</td><td>1.1489</td><td>0.2719</td><td>0.6977</td><td>21.6634</td><td>1.0255</td><td>0.2618</td></tr><tr><td>11</td><td>12.1036</td><td>3.9324</td><td>0.3994</td><td>1.2175</td><td>0.2483</td><td>0.7228</td><td>10.7362</td><td>0.8938</td><td>0.2523</td></tr></tbody></table>
4.3. DELINEATION OF FARM TYPES

65

Value_Index 46.1796 0.8779 0.2938 4.0964 0.531 NA 0.7892
Dunn Hubert SDindex Dindex SDbw
Number Clusters 15.0000 0 10.0000 0 15.0000
Value Index 0.1068 0 9.8162 0 0.1942

<table><thead><tr><th colspan="20">$Best-partition</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><td>1</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>3</td><td>2</td><td>4</td><td>2</td><td>3</td><td>2</td><td>2</td><td>2</td><td>2</td><td>1</td><td>1</td><td>2</td></tr><tr><td>21</td><td>23</td><td>24</td><td>25</td><td>26</td><td>27</td><td>28</td><td>29</td><td>30</td><td>31</td><td>32</td><td>34</td><td>35</td><td>36</td><td>37</td><td>38</td><td>39</td><td>40</td><td>41</td><td>42</td></tr><tr><td>2</td><td>2</td><td>3</td><td>2</td><td>1</td><td>3</td><td>3</td><td>1</td><td>3</td><td>4</td><td>3</td><td>4</td><td>3</td><td>3</td><td>4</td><td>3</td><td>4</td><td>2</td><td>4</td><td>1</td></tr><tr><td>43</td><td>44</td><td>45</td><td>46</td><td>47</td><td>48</td><td>49</td><td>51</td><td>52</td><td>53</td><td>54</td><td>55</td><td>56</td><td>57</td><td>58</td><td>59</td><td>60</td><td>61</td><td>63</td><td>64</td></tr><tr><td>3</td><td>2</td><td>2</td><td>1</td><td>4</td><td>1</td><td>2</td><td>1</td><td>3</td><td>3</td><td>2</td><td>4</td><td>3</td><td>4</td><td>1</td><td>2</td><td>4</td><td>3</td><td>2</td><td>4</td></tr><tr><td>65</td><td>66</td><td>67</td><td>68</td><td>69</td><td>70</td><td>71</td><td>72</td><td>74</td><td>75</td><td>76</td><td>77</td><td>78</td><td>80</td><td>81</td><td>82</td><td>83</td><td>84</td><td>85</td><td>86</td></tr><tr><td>4</td><td>1</td><td>1</td><td>4</td><td>4</td><td>1</td><td>2</td><td>4</td><td>1</td><td>3</td><td>4</td><td>4</td><td>4</td><td>1</td><td>1</td><td>1</td><td>2</td><td>1</td><td>2</td><td>1</td></tr><tr><td>87</td><td>88</td><td>89</td><td>90</td><td>91</td><td>92</td><td>93</td><td>94</td><td>95</td><td>96</td><td>97</td><td>98</td><td>99</td><td>100</td><td>101</td><td>102</td><td>103</td><td>104</td><td>105</td><td>106</td></tr><tr><td>2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>4</td><td>1</td><td>2</td><td>2</td><td>4</td><td>1</td><td>3</td><td>2</td><td>3</td><td>4</td><td>1</td><td>1</td><td>2</td><td>3</td><td>3</td></tr><tr><td>107</td><td>108</td><td>109</td><td>110</td><td>111</td><td>112</td><td>113</td><td>114</td><td>115</td><td>116</td><td>117</td><td>118</td><td>119</td><td>121</td><td>123</td><td>124</td><td>126</td><td>127</td><td>128</td><td>129</td></tr><tr><td>2</td><td>2</td><td>3</td><td>2</td><td>4</td><td>4</td><td>3</td><td>4</td><td>3</td><td>4</td><td>3</td><td>4</td><td>4</td><td>1</td><td>4</td><td>1</td><td>1</td><td>3</td><td>2</td><td>1</td></tr><tr><td>130</td><td>131</td><td>132</td><td>133</td><td>134</td><td>135</td><td>136</td><td>137</td><td>138</td><td>139</td><td>140</td><td>141</td><td>142</td><td>143</td><td>146</td><td>147</td><td>148</td><td>149</td><td>150</td><td>151</td></tr><tr><td>4</td><td>4</td><td>4</td><td>2</td><td>1</td><td>1</td><td>4</td><td>4</td><td>1</td><td>1</td><td>1</td><td>4</td><td>4</td><td>4</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>3</td></tr><tr><td>152</td><td>153</td><td>154</td><td>155</td><td>156</td><td>157</td><td>158</td><td>159</td><td>160</td><td>161</td><td>162</td><td>163</td><td>164</td><td>165</td><td>166</td><td>167</td><td>168</td><td>169</td><td>170</td><td>171</td></tr><tr><td>4</td><td>2</td><td>1</td><td>3</td><td>4</td><td>1</td><td>4</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>1</td><td>4</td><td>4</td><td>1</td><td>4</td><td>1</td><td>2</td></tr><tr><td>172</td><td>173</td><td>174</td><td>175</td><td>176</td><td>177</td><td>178</td><td>179</td><td>180</td><td>181</td><td>182</td><td>183</td><td>184</td><td>185</td><td>186</td><td>187</td><td>188</td><td>189</td><td>190</td><td>191</td></tr><tr><td>1</td><td>1</td><td>2</td><td>3</td><td>4</td><td>3</td><td>2</td><td>2</td><td>3</td><td>4</td><td>3</td><td>3</td><td>3</td><td>4</td><td>3</td><td>3</td><td>2</td><td>4</td><td>4</td><td>4</td></tr><tr><td>192</td><td>193</td><td>194</td><td>196</td><td>197</td><td>198</td><td>199</td><td>200</td><td>201</td><td>202</td><td>203</td><td>204</td><td>205</td><td>206</td><td>207</td><td>208</td><td>209</td><td>210</td><td>211</td><td>212</td></tr><tr><td>4</td><td>2</td><td>4</td><td>4</td><td>4</td><td>3</td><td>3</td><td>3</td><td>1</td><td>2</td><td>2</td><td>2</td><td>1</td><td>1</td><td>3</td><td>4</td><td>3</td><td>3</td><td>4</td><td>3</td></tr><tr><td>213</td><td>214</td><td>215</td><td>216</td><td>217</td><td>218</td><td>219</td><td>220</td><td>221</td><td>222</td><td>223</td><td>224</td><td>225</td><td>226</td><td>227</td><td>228</td><td>229</td><td>230</td><td>231</td><td>232</td></tr><tr><td>3</td><td>3</td><td>3</td><td>4</td><td>3</td><td>3</td><td>2</td><td>3</td><td>2</td><td>2</td><td>2</td><td>2</td><td>1</td><td>4</td><td>2</td><td>2</td><td>3</td><td>2</td><td>2</td><td>3</td></tr><tr><td>233</td><td>234</td><td>235</td><td>236</td><td>237</td><td>238</td><td>239</td><td>240</td><td>241</td><td>242</td><td>243</td><td>244</td><td>245</td><td>246</td><td>247</td><td>248</td><td>249</td><td>250</td><td>251</td><td>252</td></tr><tr><td>3</td><td>2</td><td>1</td><td>2</td><td>4</td><td>4</td><td>4</td><td>1</td><td>1</td><td>3</td><td>1</td><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td><td>3</td><td>4</td><td>1</td></tr><tr><td>253</td><td>254</td><td>255</td><td>256</td><td>257</td><td>258</td><td>259</td><td>260</td><td>261</td><td>262</td><td>263</td><td>264</td><td>265</td><td>266</td><td>267</td><td>268</td><td>269</td><td></td><td></td><td></td></tr><tr><td>3</td><td>3</td><td>2</td><td>3</td><td>1</td><td>3</td><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>4</td><td>4</td><td></td><td></td><td></td></tr></tbody></table>

We plot the dendrogram using the suggested optimum number of clusters of 4, visualized with functions from the R package **dendextend** (Galili, 2015). We reset the plotting window to default before that (single plot)

```r
par(mfrow = c(1, 1))
```

```r
hdend = as.dendrogram(hc_pco)
```

```r
hdend = color Fetures(hdend, k = 4)
```

```r
hdend = color labels(hdend, k = 4)
```

```r
plot(hdend)
```
<figure><img src="image_29.png" /></figure>

grpPCO = cutree(hc_pco, k = 4)

We plot the clusters on PCO planes.

plot(pco$points[, 1], pco$points[, 2], col = grpPCO)
<figure><img src="image_30.png" /></figure>

plot(pco$points[, 1], pco$points[, 3], col = grpPCo)

<figure><img src="image_31.png" /></figure>

plot(pco$points[, 2], pco$points[, 3], col = grpPCo)
<figure><img src="image_32.png" /></figure>

We add 'Type' in our dataset.

```
typo$type = grpPCO
```

```
data = merge(data, typo[, c(1, 41)], by = "ID")
```

## 4.4 Interpretation of the farm typology

### 4.4.1 Random forest

To understand which variables are the most segregating in our typology, we fit a random forest classification model and assess variable importance through the Mean Decrease in Gini index using the R package `randomForest` (Liaw and Wiener, 2002).

In our case, the prediction is good, with a low out of bage estimate of error rate of 10.75%, and class errors for the different types between 5.19% and 18.18%.

```
data$type = as.factor(data$type)
```

```
rf_type = randomForest(type ~ ., data = data[, -c(1:3)], ntree = 1500)
print(rf_type)
```

Call:
randomForest(formula = type ~ ., data = data[, -c(1:3)], ntree = 1500)
  Type of random forest: classification
  Number of trees: 1500
No. of variables tried at each split: 6

OOB estimate of error rate: 15.95%

Confusion matrix:
  1  2  3  4 class.error
1 43 10  1  9  0.31746032
2  3 63  3  0  0.08695652
3  0  0 55  5  0.08333333
4  3  0  7 55  0.15384615

Next, we plot variable importance.

The plot shows that the most segregating variables in our typology are farming
as the main source of income, and total diversity (crop and livestock).

varImpPlot(rf_type)

<figure><img src="image_33.png" /></figure>
4.4.2 Summary table

To further interprete the typology, we produce a table from variables used in the
delineation of the typology, displaying means and standard deviations for con-
tinuous ones and % for binary ones for the different types, using tbl_summary()
from the R package gtsummary (see previous chapter.

add_p() adds a column of statistical test, and add_overall() values for the whole
population.

data[, -c(1:3)] %>
  tbl_summary(by = type, statistic = list(all continuous() ~ "{mean}")
  ({sd}), all categorica l() ~
  "{p}%"), type = list(hdds ~ "continuous"), digits =
  all continuous() ~ 1,
  label = list(age_hhh ~ "Age of the head of the household",
  sex_hhh ~ "Female-headed households",
  edu_hhh ~ "Education of the head of the household higher than
  primary",
  family_size ~ "Family size", eq_value ~ "Equipment value
  (USD)", cropped_area ~
  "Total cropped area (ha)", fallow ~ "Fallow land (ha)",
  prop_non cereals ~
  "Non-cereal crops (% total cropped area)", fertilizers ~
  "Total fertilizer used (kg)",
  amendment ~ "Total manure + compost used (kg)", cattle ~
  "Cattle (n)",
  small_rum ~ "Small ruminants (n)", poultry ~ "Poultry (n)",
  bee_hives ~
  "Bee keeping", ind garden ~ "Owning an individual
garden", com garden ~
  "Having access to a communal garden", main food source ~
  "Own production as main source of food",
  main income source ~ "Farming as main source of income",
  food security ~
  "Proportion of the year being food secured", hdds ~ "24H
  household dietary diversity score (0-12)",
  community_seed_banks ~ "Using a community seed bank",
  small grains ~
  "Using small grains", crop rotation ~ "Using crop
  rotation", intercropping ~
  "Using intercropping", cover crops ~ "Using cover crops",
  mulching ~
  "Using mulching", integrated pest management ~ "Using
  integrated pest management",
  compost manure ~ "Using compost and manure",
  tree planting on farm ~
  "Planting trees on-farm",
  tree retention natural regeneration on farm ~
  "Retaining naturally regenerated trees on-farm",
  homemade feeds ~
```
"Using homemade animal feed", fodder_production ~
  "Producing fodder",
    fodder_preservation ~ "Preserving fodder", survivalaten ~
  "Using survival feeding",
    tot_div ~ "Total farm diversity (n species)", cereal produ ~
  "Total cereal production (kg/yr)",
    offtake ~ "Total livestock offtake (TLU/yr)")) %>
add_p(test = list(allcontinuous() ~ "kruskal.test",
  all_categorical() ~ "chisq.test")) %>
add_overall() %>
as_kable(format = "latex", booktabs = TRUE, longtable = TRUE) %>
kable_styling(latex_options = c("scale_down", "repeat_header"),
  font_size = 6) %>
landscape()
```
<table><thead><tr><th>**Characteristic**</th><th>**Overall** N = 257</th><th>**1** N = 63</th><th>**2** N = 69</th><th>**3** N = 60</th><th>**4** N = 65</th><th>**p-value**</th></tr></thead><tbody><tr><td>Age of the head of the household</td><td>51.6 (15.8)</td><td>50.8 (13.3)</td><td>50.3 (16.5)</td><td>57.0 (17.5)</td><td>48.8 (14.7)</td><td>0.035</td></tr><tr><td>Female-headed households</td><td>24%</td><td>17%</td><td>35%</td><td>27%</td><td>15%</td><td>0.032</td></tr><tr><td>Education of the head of the household higher than primary</td><td>38%</td><td>48%</td><td>39%</td><td>5.0%</td><td>58%</td><td>&lt;0.001</td></tr><tr><td>Family size</td><td>5.4 (2.5)</td><td>5.9 (2.5)</td><td>4.5 (2.1)</td><td>4.9 (2.1)</td><td>6.4 (3.0)</td><td>&lt;0.001</td></tr><tr><td>Total cropped area (ha)</td><td>1.7 (1.4)</td><td>1.5 (0.9)</td><td>0.9 (0.9)</td><td>1.9 (1.4)</td><td>2.7 (1.7)</td><td>&lt;0.001</td></tr><tr><td>Fallow land (ha)</td><td>0.7 (1.4)</td><td>1.0 (1.1)</td><td>0.3 (0.7)</td><td>0.7 (1.4)</td><td>1.0 (1.9)</td><td>&lt;0.001</td></tr><tr><td>Non-cereal crops (% total cropped area)</td><td>29.3 (29.1)</td><td>30.6 (23.6)</td><td>15.7 (27.8)</td><td>31.4 (32.0)</td><td>40.5 (27.3)</td><td>&lt;0.001</td></tr><tr><td>Total fertilizer used (kg)</td><td>85.6 (115.9)</td><td>77.9 (69.3)</td><td>32.4 (56.9)</td><td>59.7 (72.5)</td><td>173.7 (170.7)</td><td>&lt;0.001</td></tr><tr><td>Total manure + compost used (kg)</td><td>258.0 (543.0)</td><td>328.7 (384.9)</td><td>29.5 (71.5)</td><td>148.5 (454.2)</td><td>533.2 (831.1)</td><td>&lt;0.001</td></tr><tr><td>Cattle (n)</td><td>3.6 (6.7)</td><td>2.7 (3.6)</td><td>0.6 (2.2)</td><td>2.0 (3.3)</td><td>9.2 (10.4)</td><td>&lt;0.001</td></tr><tr><td>Small ruminants (n)</td><td>8.5 (11.2)</td><td>9.9 (12.1)</td><td>2.9 (3.9)</td><td>4.1 (4.7)</td><td>17.0 (14.1)</td><td>&lt;0.001</td></tr><tr><td>Poultry (n)</td><td>9.4 (11.3)</td><td>10.4 (9.7)</td><td>4.0 (5.2)</td><td>6.8 (11.1)</td><td>16.6 (13.5)</td><td>&lt;0.001</td></tr><tr><td>Bee keeping</td><td>11%</td><td>11%</td><td>2.9%</td><td>17%</td><td>12%</td><td>0.075</td></tr><tr><td>Owning an individual garden</td><td>47%</td><td>95%</td><td>23%</td><td>13%</td><td>55%</td><td>&lt;0.001</td></tr><tr><td>Having access to a communal garden</td><td>1.2%</td><td>1.6%</td><td>0%</td><td>1.7%</td><td>1.5%</td><td>0.8</td></tr><tr><td>Own production as main source of food</td><td>75%</td><td>62%</td><td>57%</td><td>87%</td><td>95%</td><td>&lt;0.001</td></tr><tr><td>Farming as main source of income</td><td>51%</td><td>13%</td><td>4.3%</td><td>100%</td><td>91%</td><td>&lt;0.001</td></tr><tr><td>Proportion of the year being food secured</td><td>0.5 (0.3)</td><td>0.6 (0.2)</td><td>0.4 (0.3)</td><td>0.5 (0.3)</td><td>0.7 (0.2)</td><td>&lt;0.001</td></tr><tr><td>24H household dietary diversity score (0-12)</td><td>3.1 (1.0)</td><td>3.3 (0.9)</td><td>2.8 (0.8)</td><td>2.7 (1.0)</td><td>3.6 (1.2)</td><td>&lt;0.001</td></tr><tr><td>Using a community seed bank</td><td>12%</td><td>6.3%</td><td>5.8%</td><td>23%</td><td>12%</td><td>0.007</td></tr><tr><td>Using small grains</td><td>86%</td><td>89%</td><td>72%</td><td>93%</td><td>92%</td><td>0.001</td></tr><tr><td>Using crop rotation</td><td>50%</td><td>62%</td><td>30%</td><td>37%</td><td>72%</td><td>&lt;0.001</td></tr><tr><td>Using intercropping</td><td>16%</td><td>27%</td><td>13%</td><td>0%</td><td>25%</td><td>&lt;0.001</td></tr><tr><td>Using cover crops</td><td>8.2%</td><td>6.3%</td><td>5.8%</td><td>13%</td><td>7.7%</td><td>0.4</td></tr><tr><td>Using mulching</td><td>9.7%</td><td>13%</td><td>1.4%</td><td>8.3%</td><td>17%</td><td>0.019</td></tr><tr><td>Using integrated pest management</td><td>49%</td><td>57%</td><td>25%</td><td>47%</td><td>71%</td><td>&lt;0.001</td></tr><tr><td>Using compost and manure</td><td>16%</td><td>16%</td><td>2.9%</td><td>17%</td><td>28%</td><td>0.001</td></tr><tr><td>Planting trees on-farm</td><td>14%</td><td>21%</td><td>5.8%</td><td>13%</td><td>15%</td><td>0.094</td></tr><tr><td>Retaining naturally regenerated trees on-farm</td><td>36%</td><td>43%</td><td>32%</td><td>23%</td><td>45%</td><td>0.044</td></tr><tr><td>Using homemade animal feed</td><td>25%</td><td>30%</td><td>7.2%</td><td>17%</td><td>45%</td><td>&lt;0.001</td></tr><tr><td>Producing fodder</td><td>3.9%</td><td>6.3%</td><td>0%</td><td>3.3%</td><td>6.2%</td><td>0.2</td></tr><tr><td>Preserving fodder</td><td>3.1%</td><td>3.2%</td><td>0%</td><td>0%</td><td>9.2%</td><td>0.007</td></tr><tr><td>Using survival feeding</td><td>3.9%</td><td>4.8%</td><td>0%</td><td>1.7%</td><td>9.2%</td><td>0.034</td></tr><tr><td>Total farm diversity (n species)</td><td>6.7 (3.5)</td><td>9.1 (2.1)</td><td>3.8 (1.8)</td><td>4.8 (2.4)</td><td>9.1 (3.6)</td><td>&lt;0.001</td></tr><tr><td>Total cereal production (kg/yr)</td><td>839.4 (923.3)</td><td>734.2 (790.5)</td><td>410.5 (454.5)</td><td>840.5 (773.1)</td><td>1,395.8 (1,228.9)</td><td>&lt;0.001</td></tr><tr><td>Total livestock offtake (TLU/yr)</td><td>0.3 (0.5)</td><td>0.2 (0.4)</td><td>0.1 (0.1)</td><td>0.1 (0.3)</td><td>0.6 (0.8)</td><td>&lt;0.001</td></tr></tbody></table>
(continued)

<table><thead><tr><th>**Characteristic**</th><th>**Overall** N = 257</th><th>**1** N = 63</th><th>**2** N = 69</th><th>**3** N = 60</th><th>**4** N = 65</th><th>**p-value**</th></tr></thead><tbody><tr><td>Equipment value (USD)</td><td>231.5 (296.5)</td><td>207.9 (250.5)</td><td>58.0 (91.8)</td><td>202.5 (275.0)</td><td>465.4 (350.4)</td><td>&lt;0.001</td></tr></tbody></table>
### 4.4.3 Most representative farm of each type

Finally, we pursue our interpretation of farm types by analyzing the most representative farm of each class, as the farm closest to the medoid of its farm type.

For that, we first compute distances on PCO variables.

```
dist Mat = dist(pco_var)
```

We then extract medoid (e.g., the most central observation) for each cluster/farm type, and we display IDs of these medoids (row numbers from 'data')

```
medoids = sapply unique grpPCO, function cl)
  which min colSums as matrix dist Mat[grpPCO == cl, grpPCO == cl])
medoids
```

```
243 17 207 156
55 13 42 38
```

We then retrieve the full farm data for these medoids/most representative farms.

```
most_representative farms = data medoids, ]
most_representative farms
```

<table><thead><tr><th>ID</th><th>ward</th><th>district</th><th>age_hhh</th><th>sex_hhh</th><th>edu_hhh</th><th>family_size</th><th>cropped_area</th></tr></thead><tbody><tr><td>55</td><td>363</td><td>Ward 3</td><td>9 12 17</td><td>Mbire</td><td>57</td><td>0</td><td>0</td><td>5</td><td>1.0</td></tr><tr><td>13</td><td>318</td><td>Ward 2</td><td></td><td>Mbire</td><td>69</td><td>0</td><td>0</td><td>8</td><td>0.5</td></tr><tr><td>42</td><td>349</td><td>Ward 3</td><td>9 12 17</td><td>Mbire</td><td>39</td><td>0</td><td>1</td><td>4</td><td>1.0</td></tr><tr><td>38</td><td>345</td><td>Ward 3</td><td>9 12 17</td><td>Mbire</td><td>60</td><td>0</td><td>0</td><td>7</td><td>0.0</td></tr><tr><td></td><td></td><td></td><td></td><td colspan="5">fallow prop_non_cereals fertilizers amendment cattle small_rum poultry</td></tr><tr><td>55</td><td>2.0</td><td></td><td></td><td>0.00000</td><td>50</td><td>600</td><td>0</td><td>5</td><td>0</td></tr><tr><td>13</td><td>0.0</td><td></td><td></td><td>0.00000</td><td>5</td><td>0</td><td>0</td><td>0</td><td>20</td></tr><tr><td>42</td><td>0.8</td><td></td><td></td><td>19.98002</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>38</td><td>1.2</td><td></td><td></td><td>0.00000</td><td>10</td><td>250</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td></td><td></td><td></td><td colspan="5">bee_hives ind_garden com_garden main_food_source main_income_source</td></tr><tr><td>55</td><td></td><td>1</td><td></td><td>1</td><td>0</td><td></td><td>0</td><td></td><td>0</td></tr><tr><td>13</td><td></td><td>1</td><td></td><td>0</td><td>0</td><td></td><td>0</td><td></td><td>1</td></tr><tr><td>42</td><td></td><td>0</td><td></td><td>0</td><td>0</td><td></td><td>1</td><td></td><td>0</td></tr><tr><td>38</td><td></td><td>0</td><td></td><td>1</td><td>0</td><td></td><td>0</td><td></td><td>0</td></tr><tr><td></td><td></td><td></td><td></td><td colspan="5">food_security hdds community_seed_banks small_grains crop_rotateion</td></tr><tr><td>55</td><td></td><td>0.7500000</td><td>2</td><td></td><td>1</td><td></td><td>1</td><td></td><td>1</td></tr><tr><td>13</td><td></td><td>0.3333333</td><td>3</td><td></td><td>0</td><td></td><td>1</td><td></td><td>1</td></tr></tbody></table>
<table><tbody><tr><td>42</td><td>0.5000000</td><td>3</td><td>0</td><td>1</td><td>1</td></tr><tr><td>38</td><td>0.0000000</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td colspan="6">intercropping cover_crops mulching integrated_pest_management compost_manure</td></tr><tr><td>55</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>13</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>42</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>38</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td colspan="6">tree_planting_on_farm tree_retention_natural_regeneration_on_farm</td></tr><tr><td>55</td><td>1</td><td></td><td></td><td>0</td><td></td></tr><tr><td>13</td><td>0</td><td></td><td></td><td>0</td><td></td></tr><tr><td>42</td><td>0</td><td></td><td></td><td>1</td><td></td></tr><tr><td>38</td><td>0</td><td></td><td></td><td>0</td><td></td></tr><tr><td colspan="6">homemade SOLID FEDERATION</td></tr><tr><td>55</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td></tr><tr><td>13</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>42</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>38</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td colspan="6">tot_div cereal SOLID FEDERATION</td></tr><tr><td>55</td><td>10</td><td>257</td><td>0.20</td><td>150</td><td>1</td></tr><tr><td>13</td><td>4</td><td>170</td><td>0.02</td><td>0</td><td>3</td></tr><tr><td>42</td><td>5</td><td>550</td><td>0.00</td><td>50</td><td>2</td></tr><tr><td>38</td><td>3</td><td>0</td><td>0.00</td><td>100</td><td>2</td></tr></tbody></table>

4.5 Selection of stratified samples

Typologies can be used to select a sample of farms representative of the overall population to e.g., host on-farm experiments.

Below, we select 30 representative farms using farm type as single stratum.

sample1 = as.data.frame(stratified(data, "type", 30/nrow(data)))
table(sample1 type)

1 2 3 4
7 8 7 8

We can also use more than one stratum: here we select 30 representative farms using ward (an administrative unit) and farm type as strata.

sample2 = as.data.frame(stratified(data, c("ward", "type"),
  30/nrow(data)))
table(sample2 type, sample2 ward)
<table><thead><tr><th></th><th>Ward 2</th><th>Ward 3</th><th>Ward 9</th><th>Ward 12</th><th>Ward 17</th></tr></thead><tbody><tr><th>1</th><td>4</td><td></td><td></td><td></td><td>3</td></tr><tr><th>2</th><td>5</td><td></td><td></td><td></td><td>3</td></tr><tr><th>3</th><td>2</td><td></td><td></td><td></td><td>5</td></tr><tr><th>4</th><td>3</td><td></td><td></td><td></td><td>4</td></tr></tbody></table>
Chapter 5

Managing on-farm
experiments

In this chapter, we will go through the steps needed to collect data from on-farm experiments - where farms/trials and plots/treatments are identified by QR codes to reduce misidentification errors - using KoboToolBox for data collection and R for data handling and analysis.

The use of QR codes provides an efficient and reliable system for identify-
ing farms/trials and plots/treatments during data collection in on-farm experi-
ments.

Each QR code encodes a unique identifier that corresponds to a farmer/experimental site or a specific plot/treatment. During field monitoring, these codes can be scanned directly in KoboToolbox, automatically filling in the correct ID in the data collection form. This approach minimizes transcription and labeling errors, ensures consistent tracking of repeated measurements, and facilitates seamless integration of field data with analytical workflows in R.

In the CGIAR Agroecology initiative in Zimbabwe, each farmer hosting a trial received an ID card with a unique QR code. Before any data collection on a given farm/trial, the card was scanned to confirm the correct trial site and verify that the farmer was present. Each plot within a trial also had its own signboard displaying a QR code linked to its treatment, ensuring every observation was associated with the correct experimental condition.

5.1 Generating QR codes

We will now demonstrate how to generate unique QR codes for both farms/trials
and plots/treatments using the R package **baRcodeR** (Wu and Colautti, 2022).
<figure><img src="image_34.png" /><figcaption>Figure 5.1: ID card with unique QR code received by each farmer hosting a trial in the CGIAR Agroecology initiative in Zimbabwe (left) and signage with a unique QR code for each plot/treatment in the corresponding trial</figcaption></figure>

```r
rm(list = ls())

if (!require("baRcodeR")) install.packages("baRcodeR")

library(baRcodeR)
```

For demonstration and data protection compliance, we will use fictional farmer names that do not correspond to any real individuals. This ensures that no personal information is used while illustrating the workflow.

```r
set.seed(123)

given = c("Tendai", "Fungai", "Tawanda", "Ropafadzo", "Tafadzwa",
  "Nyasha", "Chiedza",
  "Anesu", "Munyaradzi", "Chenai", "Kudzai", "Tariro", "Simba",
  "Rutendo", "Kudzanai",
  "Vimbai", "Tanaka", "Blessing", "Shamiso", "Rudo", "Kuda", "Tinashe",
  "Chipo",
  "Prosper", "Sekai", "Loveness", "Godfrey", "Mercy", "Brenna",
  "Nokutenda", "Sibusiso",
  "Vusumuzi", "Eliazar", "Palesa", "Tafara", "Mandla", "Zvikomborero",
  "Tanyaradzwa")

surname = c("Mupfumi", "Chikwepa", "Mazarura", "Nyamadzawo", "Chikora",
  "Makarau",
  "Gwekwerere", "Mandizha", "Sibanda", "Gonye", "Mtsvera", "Ndlovu",
  "Moyo", "Ncube",
  "Mapfumo", "Chidavaenzi", "Mutsvangwa", "Dube", "Mashingaidze",
  "Nyamapfene",
  "Gunda", "Mungano", "Chikowore", "Matambanadzo", "Zimuto", "Mabhunu",
  "Chinyoka",
```
"Chimanikire", "Chiriseri", "Gwanyanya", "Madziva", "Mapiravana",
→ "Gwatidzo",
"Chikore", "Mashingairi", "Chirwa", "Chinyamurindi")

n = 25

fake_names = paste(sample(given, n, replace = TRUE), sample(surname, n,
→ replace = TRUE))

fake_names

[1] "Sibusiso Mapiravana"
[4] "Tawanda Gonye"
[7] "Sekai Gwekwerere"
[10] "Tafadzwa Zimuto"
[13] "Munyaradzi Chikora"
[16] "Anesu Moyo"
[19] "Munyaradzi Chinyoka"
[22] "Rutendo Mapfumo"
[25] "Kudzanai Chidavaenzi"

"Kudzanai Gwekwerere"
"Zvikomborero Chikowore"
"Loveness Chinyoka"
"Godfrey Chikore"
"Brenna Mandizha"
"Loveness Dube"
"Shamiso Zimuto"
"Tanaka Mabhunu"

"Rutendo Sibanda"
"Rutendo Chinyoka"
"Godfrey Mapiravana"
"Mercy Chiriseri"
"Tafara Ndlovu"
"Chiedza Gwatidzo"
"Mandla Gunda"
"Tariro Madziva"

5.1.1 Generating QR codes for farmers/trials

We now generate a first set of QR codes corresponding to farmers or trials. Each
QR code encodes one of the fictional farmer names.

The create_PDF() function from baRcodeR generates a printable PDF con-
taining all QR codes and their labels, which can later be laminated or used on
ID cards in the field (as represented in Figure 5.1 above).

create_PDF(Labels = fake_names, numcol = 2, numrow = 4, Fsz = 8, name =
"farms_qrcodes")

knitr::include(
<figure><img src="image_35.png" /></figure>

Sibusiso Mapira
vana

<figure><img src="image_36.png" /></figure>

Kudzanai Gwekwe
rere

<figure><img src="image_37.png" /></figure>

Rutendo Sibanda

<figure><img src="image_38.png" /></figure>

Tawanda Gonye

<figure><img src="image_39.png" /></figure>

Zvikomborero Ch
ikowore

<figure><img src="image_40.png" /></figure>

Rutendo Chinyok
a

<figure><img src="image_41.png" /></figure>

Sekai Gwekwerer
e

<figure><img src="image_42.png" /></figure>

Loveness Chinyo
ka

Once the PDF is created, each farmer has an identification label that can be scanned during data collection. Scanning the QR code in Kobo Toolbox automatically retrieves the farmer's unique ID, ensuring that all data collected belong to the correct individual/site.
5.1.2 Generating QR codes for plots/treatments

We will now generate QR codes for every combination of farmer × treatment. Each farm typically hosts multiple treatments or plots. In this example, there are three treatments per trial: ‘conventional practice’ (CONV), ‘conservation agriculture’ (CA), and ‘push-pull’ (PPULL). The *expand.grid()* function creates all combinations, which are then concatenated into unique identifiers such as: ‘Tendai Mupfumi-CONV’, ‘Tendai Mupfumi-CA’, ‘Tendai Mupfumi-PPULL’, etc. Each of these labels corresponds to a unique plot within a farmer’s trial.

l = expand.grid(FR = fake_names, TR = c("CONV", "CA", "PPULL"))
l = paste(l$FR, l$TR, sep = "-")

create_PDF(Labels = 1, numcol = 2, numrow = 4, Fsz = 8, name =
"plots_qrcodes")

knitr::include("plots_qrcodes.pdf")
<figure><img src="image_43.png" /></figure>

Sibusiso Mapira
vana-CONV

<figure><img src="image_44.png" /></figure>

Kudzanai Gwekwe
rere-CONV

<figure><img src="image_45.png" /></figure>

Rutendo Sibanda
-CONV

<figure><img src="image_46.png" /></figure>

Tawanda Gonye-C
ONV

<figure><img src="image_47.png" /></figure>

Zvikomborero Ch
ikowore-CONV

<figure><img src="image_48.png" /></figure>

Rutendo Chinyok
a-CONV

<figure><img src="image_49.png" /></figure>

Sekai Gwekwerer
e-CONV

<figure><img src="image_50.png" /></figure>

Loveness Chinyo
ka-CONV

The resulting PDF contains QR codes for all farmer × treatment combinations, which can be printed on plot signboards (as represented in Figure 5.1 above). During fieldwork, scanning the plot's QR code ensures that measurements (e.g., crop establishment, assessment of pest damage, estimation of yield) are automatically attributed to the correct farm and treatment.
## 5.2 Integrating QR codes in KoboToolBox surveys

QR codes can be easily integrated into KoboToolbox surveys to automatically record the identity of farms/trials and plots/treatments during data collection. This ensures that each observation is correctly linked to its corresponding farmer/trial and plot/treatment without manual typing or risk of mislabeling.

As an example, consider the survey questionnaire used to monitor on-farm experiments under the CGIAR Agroecology Initiative in Zimbabwe. The corresponding XLSForm (used to build the questionnaire in KoboToolbox) is available here: https://github.com/FBaudron/Supporting-codesign/blob/main/05-managing-experiments/EXPERIMENT%20SURVEY%20TOOL.xlsx. A PDF version of the questionnaire is also provided: https://github.com/FBaudron/Supporting-codesign/blob/main/05-managing-experiments/EXPERIMENT%20SURVEY%20TOOL.pdf. You can adapt this survey questionnaire to your own use.

In this questionnaire, several variables are of type 'barcode': 'farmID', 'plotID_1', 'plotID_2', and 'plotID_3' (each trial consists of three treatments arranged in three separate plots). The three treatments are: 'conventional practice' (CONV), 'conservation agriculture' (CA), and 'push-pull' (PPULL).

Each of these variables corresponds to a question where the enumerator scans a QR code instead of typing an ID. The farmID is obtained by scanning the farmer's ID card, while each plotID corresponds to the QR code placed on a specific plot signboard (Figure 5.1).

In the Kobo form, the questions for the three plots are grouped between 'begin_repeated' and 'end_repeated' statements. This looping structure allows enumerators to enter repeated measurements (e.g., pest assessment in several points in each plot, typically 5, following a W-shaped sampling pattern) efficiently and consistently.

### 5.2.1 Data from KoboToolBox surveys that include QR codes

Once data have been collected in KoboToolbox, they can be exported in Excel format for analysis in R. To illustrate this, you can download a sample dataset generated from the survey described above here: https://github.com/FBaudron/Supporting-codesign/blob/main/05-managing-experiments/EXPERIMENT%20SURVEY%20DATA.xlsx.

This dataset was downloaded directly from KoboToolbox as an Excel file. It contain pest damage data collected on sorghum over one season in three different plots (corresponding to the different treatments described above). The
file is made of four sheets: one main sheet with general information (e.g., farm ID, enumerator name, date), and three treatment sheets, each corresponding to pest assessment data in one of the three plots recorded through the ‘begin_re repeat’ / ‘end_re repeat’ structure in the questionnaire.

A standalone R script containing the full workflow below is avail-
able at the following link: https://github.com/FBaudron/Supporting-
codesign/blob/main/05-managing-experiments/managing_experiments.R.

Before importing the data in R, we first clear the workspace and ensure that
the necessary packages are installed and loaded.

rm(list = ls())

if (!require("openxlsx")) install.packages("openxlsx")
if (!require("tidyr")) install.packages("tidyr")
if (!require("lmerTest")) install.packages("lmerTest")
if (!require("performance")) install.packages("performance")
if (!require("ggeffects")) install.packages("ggeffects")
if (!require("sjPlot")) install.packages("sjPlot")

library(openxlsx)
library(tidyr)
library(lmerTest)
library(performance)
library(ggeffects)
library(sjPlot)

data = read.xlsx("EXPERIMENT SURVEY DATA.xlsx", sheet = 1)

sc1 = read.xlsx("EXPERIMENT SURVEY DATA.xlsx", sheet = 2)
sc2 = read.xlsx("EXPERIMENT SURVEY DATA.xlsx", sheet = 3)
sc3 = read.xlsx("EXPERIMENT SURVEY DATA.xlsx", sheet = 4)

Each plot’s QR code (scanned in the field) produces a text string such as ‘Tendai Mupfumi-CA’ or ‘Rudo Makarau-CONV’. These strings combine the farmer’s name and the treatment applied to that plot.

To make the data easier to analyze, we split each plot ID column into two new
columns: farmer - extracted from the first part of the QR code label - and
a renamed column for the treatment in that plot (e.g., ‘1st plot’, ‘2nd plot’).
Since the farmer’s name is already recorded via the farmID QR code at the
beginning of the survey, we later drop the redundant farmer column from each
plot section. This leaves a clean dataset with one row per trial, including the
farmer ID and treatments for the three plots. We also retain the submission ID
(‘_index’) to be able to later merge this data with pest assessment data.
data = separate(data = data, col = "Scan.the.QR.code.of.the.1st.plot",
  into = c("Farmer",
  "1st plot"), sep = "-")

data = data[, c(3, 17, 5:7)]

data = separate(data = data, col = "Scan.the.QR.code.of.the.2nd.plot",
  into = c("Farmer",
  "2nd plot"), sep = "-")

data = data[, c(1:3, 5:6)]

data = separate(data = data, col = "Scan.the.QR.code.of.the.3rd.plot",
  into = c("Farmer",
  "3rd plot"), sep = "-")

data = data[, c(1:4, 6)]

We split the main dataset into three separate datasets, one for each treatment. Each of these treatment-specific datasets contains three key fields: the farmer name, the treatment name, and the submission ID. Submission ID will allow merging with the repeated pest assessment data collected in KoboToolbox.

For each treatment dataframe, we reorder the columns to place the submission ID first. We also remove columns that are not needed for analysis (e.g., submission date), keeping only the fields required for analysis.

data1 = data[, c(1:3)]
data2 = data[, c(1:2, 4)]
data3 = data[, c(1:2, 5)]

sc1 = sc1[, c(11, 1:8)]
sc2 = sc2[, c(11, 1:8)]
sc3 = sc3[, c(11, 1:8)]

We merge each treatment dataset with its corresponding pest assessment data using the submission ID. In KoboToolbox exports, `_index` identifies the main submission, while `_parent_index` links repeated sections. This merge ensures that each plot's pest assessment data aligns correctly with the corresponding farmer and treatment. After merging, we rename the plot columns to 'Treatment'.

We then combine the six treatment datasets into a single unified dataset. The submission ID is removed since it is no longer needed after merging. Selected columns are renamed to clearly reflect their contents, including the 'Farmer' column (captured from the QR code scanned at the start of the form) and various pest damage scores.
```r
data1 = merge(data1, sc1, by.x = "_index", by.y = "_parent_index", all.y
             % TRUE)
data2 = merge(data2, sc2, by.x = "_index", by.y = "_parent_index", all.y
             % TRUE)
data3 = merge(data3, sc3, by.x = "_index", by.y = "_parent_index", all.y
             % TRUE)

names(data1)[3] = "Treatment"
names(data2)[3] = "Treatment"
names(data3)[3] = "Treatment"

data = rbind(data1, data2, data3)

data = data[, -c(1)]

names(data)[c(1, 4:10)] = c("Farmer", "armmoured_cricket_leaf damage",
                           "armmoured_cricket damage score",
                           "fallArmyworm_leaf damage", "fallArmyworm Frass whorl",
                           "fallArmyworm damage score",
                           "stalk_borer_leaf damage", "stalk_borer damage score")
```

### 5.2.2 Running basic statistical models on the cleaned data

Using the cleaned and merged dataset, we fit a mixed-effects model to analyze the severity of fall armyworm damage.

We first set the levels of the variable Treatment to ensure that 'CONV' (conventional practice) is used as the reference category, since R orders factor levels alphabetically by default. We then fit a generalized linear mixed-effects model to explain the variability in fall armyworm severity scores. Treatment and crop growth stage (Vstage) are included as fixed effects, while the farm (Farmer) is included as a random effect to account for variability among farms. Because the severity score is discrete (1-9), we use a Poisson distribution.

We evaluate model performance using the function *check_model()* of the R package *performance* (Lüdecke et al., 2021) to ensure the assumptions of the mixed-effects model are reasonably met.

```r
data$Treatment = factor(data$Treatment, levels = c("CONV", "CA",
                                                  "PPULL"))

mod_faw_severity = glmer(fallArmyworm damage score ~ Treatment + Vstage
                           % + (1 | Farmer),
                           data = data, family = poisson())

check_plot = check_model(mod_faw_severity)

print(check plot)
```
<figure><img src="image_51.png" /></figure>

We examine the model results using both ANOVA and the model summary to assess the significance of fixed effects and to understand the estimated coefficients.

```
anova(mod_faw_severity)
```

Analysis of Variance Table

npar Sum Sq Mean Sq F value

Treatment 2 0.0169 0.0085 0.0085

Vstage 1 5.4042 5.4042 5.4042

```
summary(mod_faw_severity)
```

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]

Family: poisson ( log )

Formula: fallArmyworm_damag score ~ Treatment + Vstage + (1 | Farmer)

Data: data
AIC BIC logLik -2*log(L) df.resid
424.5 439.3 -207.2 414.5 138

Scaled residuals:

Min 1Q Median 3Q Max
-1.01490 -0.33144 -0.04122 0.28653 3.09858

Random effects:

Groups Name Variance Std.Dev.

Farmer (Intercept) 0.07013 0.2648

Number of obs: 143, groups: Farmer, 11

Fixed effects:

Estimate Std. Error z value Pr(>z|)
(Intercept) 0.36511 0.19565 1.866 0.0620 .
TreatmentCA 0.01969 0.14880 0.132 0.8947
TreatmentPPULL 0.01554 0.14022 0.111 0.9118
Vstage 0.09504 0.04069 2.336 0.0195 *

---

Signif. codes: 0 '**' 0.001 '**' 0.01 '**' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:

(Intr) TrtmCA TPPULL

TreatmentCA -0.372

TrtmntPPULL -0.378 0.492

Vstage -0.740 0.001 0.005

We calculate predicted mean severity scores for each treatment level based on the fitted model. This allows us to interpret the expected damage under different management practices.

predict_response(mod_faw_severity, terms = "Treatment")

# Predicted counts of fallArmyworm_damag score

Treatment | Predicted | 95% CI
---
CONV | 1.98 | 1.53, 2.56
CA | 2.02 | 1.55, 2.63
PPULL | 2.01 | 1.56, 2.59

Adjusted for:

* Vstage = 3.35

* Farmer = 0 (population-level)
Finally, we visualize the predicted mean severity scores for each treatment, which helps to compare the expected pest damage across different treatments.

```r
plot_model(mod_faw_severity, type = "pred", pred.type = "fe", terms =
  ~ Treatment)
```

<figure><img src="image_52.png" /></figure>
Chapter 6

Evaluating on-farm
experiments

In this chapter, we focus on the participatory evaluation of on-farm experiments, involving farmers and other stakeholders (e.g., extension agents, NGOs, local authorities) directly in the interpretation of the results. Building on a dataset covering two seasons (2022/23 and 2023/24) of on-farm trials on pest damage, crop yield, and farmers’ perceptions, we show how results can be explored in-teractively during field days, feedback meetings, or local workshops using Shiny applications.

Shiny apps are interactive web applications built directly in R. They allow users to explore data dynamically — e.g., by selecting seasons or indicators — without requiring any programming knowledge. Because graphs and summaries update automatically when the user changes a selection, Shiny provides an intuitive interface for collective interpretation of results. These digital tools also enable farmers to see how their own fields contribute to wider patterns, creating a space where observations, interpretations, and preferences can be discussed openly. In this way, the Shiny app acts as a bridge between quantitative measurements and farmers’ qualitative, experiential knowledge, helping transform experimental results into practical insights for decision-making.

Download the dataset used throughout this chapter here: https://github.com/FBaudron/Supporting-codesign/blob/main/06-evaluating-experiments/DATA%20FOR%20SHINY.xlsx.

The data were collected on sorghum over two seasons and include pest incidence and severity, crop yield, and farmers’ overall evaluation of the three experimental treatments: conventional practice (CONV), conservation agriculture (CA), and push-pull (PPULL). The next code chunk loads the raw dataset, prepares factor variables, and separates the two seasons for later comparisons.
A standalone R script containing the full workflow below is available at the following link: https://github.com/FBaudron/Supporting-codesign/blob/main/06-evaluating-experiments/visualization_of_treatments_performance.R.

## 6.1 Visualization of treatments' performance

Before introducing the interactive application, we build the graphs that will later be embedded in the Shiny interface, using data from the 2022/23 season (the app will allow the user to switch between this season and 2023/24).

We begin by clearing the workspace, loading the required packages, importing the dataset, and splitting it into separate objects for the two seasons.

```r
rm(list = ls())

if (!require("openxlsx")) install.packages("openxlsx")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("ggdist")) install.packages("ggdist")
if (!require("ggthemes")) install.packages("ggthemes")
if (!require("egg")) install.packages("egg")
if (!require("ggradar")) install.packages("ggradar")

library(openxlsx)
library(ggplot2)
library(ggdist)
library(ggthemes)
library(egg)
library(ggradar)

data = read.xlxs("DATA FOR SHINY.xlsx", sheet = 1)

data$Treatment = factor(data$Treatment, levels = c("CONV", "CA",
"PPULL"))

data_2022_23 = subset(data, Season == "2022/23")
data_2023_24 = subset(data, Season == "2023/24")
```

### 6.1.1 Pest pressure in different treatments

We first explore pest pressure by visualising incidence and severity for three major pests: armoured cricket, fall armyworm, and stemborer. We combine half-eye distributions, boxplots, and jittered points to convey both the central tendencies and underlying variability across farmers' fields. We use the R packages `ggplot2` (Wickham, 2016) and `ggdist` (Kay, 2025).
```r
cri = data_2022_23[, c(1:3, 4:5)]
names(cri)[4] = "Incidence"
names(cri)[5] = "Severity"
```

```r
faw = data_2022_23[, c(1:3, 6:7)]
names(faw)[4] = "Incidence"
names(faw)[5] = "Severity"
```

```r
stm = data_2022_23[, c(1:3, 8:9)]
names(stm)[4] = "Incidence"
names(stm)[5] = "Severity"
```

```r
p1 = ggplot(cri, aes(x = Treatment, y = Incidence, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
  width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
  "#DD513AFF", "#932667FF")) +
  ylab("") + xlab("") + themeFew() + ggtitle("Armored cricket -
  Incidence (%)") +
  theme(legendDash = 0, plot.title = element_text(size = 16,
  face = "bold",
  margin = ggplot2::margin(0, 0, 10, 0)), axis.text =
  element_text(size = 10),
  axis.title = element盲点(), legend.title = element盲点()) +
  coordflip()
```

```r
p2 = ggplot(faw, aes(x = Treatment, y = Incidence, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
  width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
  "#DD513AFF", "#932667FF")) +
  ylab("") + xlab("") + themeFew() + ggtitle("Fall armyworm -
  Incidence (%)") +
  theme(legendDash = 0, plot.title = element_text(size = 16,
  face = "bold",
  margin = ggplot2::margin(0, 0, 10, 0)), axis.text =
  element_text(size = 10),
  axis.title = element盲点(), legend.title = element盲点()) +
  coordflip()
```
```r
p3 = ggplot(stm, aes(x = Treatment, y = Incidence, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point.shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
  width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
  "#DD513AFF", "#932667FF")) +
  ylab("") + xlab("") + theme_few() + ggtitle("Stemborer - Incidence
(%)") + theme(legendDash = "none",
  plot.title = element_text(size = 16, face = "bold", margin =
  ggplot2::margin(0,
  0, 10, 0)), axis.text = element_text(size = 10), axis.title =
  elementBlanks(),
  legend.title = elementBlanks()) + coord flip()
```

```r
p4 = ggplot(cri, aes(x = Treatment, y = Severity, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point.shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
  width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
  "#DD513AFF", "#932667FF")) +
  ylab("") + xlab("") + theme_few() + ggtitle("Armored cricket -
Severity (1-9)") +
  theme(legendDash = "none", plot.title = element_text(size = 16,
  face = "bold",
  margin = ggplot2::margin(0, 0, 10, 0)), axis.text =
  elementText(size = 10),
  axis.title = elementBlank(), legend.title = elementBlank()) +
  coord flip()
```

```r
p5 = ggplot(faw, aes(x = Treatment, y = Severity, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point.shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
  width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
  "#DD513AFF", "#932667FF")) +
  ylab("") + xlab("") + theme_few() + ggtitle("Fall armyworm - Severity
(1-9)") +
```
```r
theme(legendITA position = "none", plot.title = element_text(size = 16,
  face = "bold",
  margin = ggplot2::margin(0, 0, 10, 0)), axis.text =
  element_text(size = 10),
  axis.title = element CAT blank(), legend.title = element CAT blank()) +
  coordFlip()

p6 = ggplot(stm, aes(x = Treatment, y = Severity, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
  width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
  "#DD513AFF", "#932667FF")) +
  ylab("") + xlab("") + themeFew() + ggtitle("Stemborer - Severity
(1-9)") + theme(legendner position = "none",
  plot.title = element_text(size = 16, face = "bold", margin =
  ggplot2::margin(0,
  0, 10, 0)), axis.text = element_text(size = 10), axis.title =
  element CAT blank(),
  legend.title = element CAT blank()) + coordFlip()

ggarrange(p1, p4, p2, p5, p3, p6, ncol = 2, nrow = 3, widths = c(1, 1),
  heights = c(1,
  1, 1))
```
<figure><img src="image_53.png" /></figure>

### 6.1.2 Productivity of different treatments

In the next set of graphs, we show grain and biomass yields for each treatment. As for the previous panel of graphs, We use half-eye distributions, boxplots and jittered points.

```r
yldcrl = data_2022_23[, c(1:3, 10:11)]
yldcrl = unique(yldcrl)

p7 = ggplot(yldcrl, aes(x = Treatment, y = YLD, fill = Treatment)) +
  ggdist::stat_halfeye(adjust = 0.6,
  width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
  alpha = 0.5,
  color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
  alpha = 0.5) +
  geom_point.shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1,
```
```r
width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
"#DD513AFF", "#932667FF")) +
ylab("") + xlab("") + themeFew() + ggTitle("Cereal yield - Grain
(t/ha)") +
theme(legendDash = "none", plot.title = element_text(size = 16,
face = "bold",
margin = ggplot2::margin(0, 0, 10, 0)), axis.text =
element_text(size = 10),
axis.title = elementBlanks(), legend.title = elementBlanks()) +
coord flip()

p8 = ggplot(yldcr1, aes(x = Treatment, y = BIO, fill = Treatment)) +
 ggdist::statHalfeye(adjust = 0.6,
width = 0.4, .width = 0, justification = -0.4, pointColour = NA,
alpha = 0.5,
color = "black") + geom_boxplot(width = 0.2, outlier.shape = NA,
alpha = 0.5) +
geom_point shape = 21, size = 2, alpha = 0.5, position =
position Jitter(seed = 1,
width = 0.1)) + scale_fill_manual(values = c("#FCA50AFF",
"#DD513AFF", "#932667FF")) +
ylab("") + xlab("") + themeFew() + ggTitle("Cereal yield - Biomass
(t/ha)") +
theme(legendDash = "none", plot.title = element_text(size = 16,
face = "bold",
margin = ggplot2::margin(0, 0, 10, 0)), axis.text =
element_text(size = 10),
axis.title = elementBlanks(), legend.title = elementBlanks()) +
coord flip()

ggarrange(p7, p8, ncol = 2, nrow = 1, widths = c(1, 1), heights = c(1))
```
<figure><img src="image_54.png" /></figure>

### 6.1.3 Rating of treatments by different demographic groups

In the last set of plot, we compare how different demographic groups — older females, older males, young females, and young males — evaluated the three treatments using radar plots built with the R package `ggradar` (Bion, 2023).

The plots below display the composite scores for each group, allowing a visualization of convergences and divergences in perceptions.

```
of = data_2022_23[, c(1:2, 12:18)]
of$Group = rep("Older females", nrow(of))
of = unique(of)
om = data_2022_23[, c(1:2, 19:25)]
om$Group = rep("Older males", nrow(om))
om = unique(om)
yf = data_2022_23[, c(1:2, 26:32)]
yf$Group = rep("Young females", nrow(yf))
yf = unique(yf)
ym = data_2022_23[, c(1:2, 33:39)]
ym$Group = rep("Young males", nrow(ym))
ym = unique(ym)
names(of)[c(3:9)] = c("Grain", "Biomass", "Inputs", "Labour", "Pests",
  "Drought",
  ...
```
```
"Soil")

names(om)[c(3:9)] = c("Grain", "Biomass", "Inputs", "Labour", "Pests",
  én "Drought",
  "Soil")
```

```
names(yf)[c(3:9)] = c("Grain", "Biomass", "Inputs", "Labour", "Pests",
  én "Drought",
  "Soil")
```

```
names(ym)[c(3:9)] = c("Grain", "Biomass", "Inputs", "Labour", "Pests",
  én "Drought",
  "Soil")
```

```
rat = rbind(of, om, yf, ym)
```

```
conv = subset(rat, Treatment == "CONV")
```

```
ca = subset(rat, Treatment == "CA")
```

```
ppull = subset(rat, Treatment == "PPULL")
```

```
p9 = ggradar(conv[, c(10, 3:9)], plot.title = "Conventional Practice",
  én values.radar = c("0",
  "5", "10"), grid.min = 0, grid.mid = 5, grid.max = 10,
```

```
  én group line width = 1,
  én group.point.size = 3, group.colours = c("#FCA50AFF", "#DD513AFF",
  "#932667FF",
```

```
  "#420A68FF"), fill = TRUE, fill.alpha = 0.2,
```

```
  én background.circle.colour = "white",
  gridline.mid.colour = "grey", legendDash = "none") +
```

```
  én theme plot title = element_text(size = 18,
  face = "bold", hjust = 0.5))
```

```
p10 = ggradar(ca[, c(10, 3:9)], plot.title = "Conservation Agriculture",
  én values.radar = c("0",
  "5", "10"), grid.min = 0, grid.mid = 5, grid.max = 10,
```

```
  én group line width = 1,
  én group.point.size = 3, group.colours = c("#FCA50AFF", "#DD513AFF",
  "#932667FF",
```

```
  "#420A68FF"), fill = TRUE, fill.alpha = 0.2,
```

```
  én background.circle.colour = "white",
  gridline.mid.colour = "grey", legendDash = "bottom") +
```

```
  én theme plot title = element_text(size = 18,
  face = "bold", hjust = 0.5))
```

```
p11 = ggradar(ppull[, c(10, 3:9)], plot.title = "Push-Pull", values.radar
  én c("0",
  "5", "10"), grid.min = 0, grid.mid = 5, grid.max = 10,
```

```
  én group line width = 1,
  én group.point.size = 3, group.colours = c("#FCA50AFF", "#DD513AFF",
  "#932667FF",
```
```r
"#420A68FF"), fill = TRUE, fill.alpha = 0.2,
  background.circle.colour = "white",
  gridline.mid.colour = "grey", legendposition = "none") +
  theme plot.title = element_text(size = 18,
  face = "bold", hjust = 0.5))

ggarrange(p9, p10, p11, ncol = 3, nrow = 1, widths = c(1, 1, 1), heights
  = c(1))
```

<figure><img src="image_55.png" /></figure>

## 6.2 Shiny app

### 6.2.1 Description of a Shiny app

A Shiny application is built around two core components: the user interface and the server logic. The user interface (ui) defines what participants see on the screen: dropdown menus (e.g., to choose a season), tabs to navigate between indicators, and the spaces where plots or tables will appear. It provides the visual structure of the app and determines how users interact with the data.

The server component is where the analytical work takes place. It receives the user's selections, filters or processes the data accordingly, performs calculations or summarisation, and generates the updated plots and tables that appear on the screen.

Several R packages support this process. The core package is **shiny** (Chang et al., 2024), which provides the framework for building reactive web applications. Two additional packages improve the user experience of the app. **bslib** (Sievert et al., 2025) allows customisation of the app's appearance using Bootstrap themes, making the layout cleaner and compatible across devices. **shinycssloaders** (Attali and Sali, 2024) adds animated loading indicators ('spinners') to plots so that users immediately see when a graph is being recalculated, which is helpful when working with larger datasets.
What makes Shiny powerful is its use of reactive programming. Whenever a user changes a selection—choosing a different indicator or switching seasons—the app automatically re-runs only the parts of the code affected by that change and updates the outputs instantly. This creates the experience of a live, responsive system, even though all processing occurs within R. Thanks to this structure, Shiny apps can serve as accessible interfaces to complex datasets, allowing farmers and other stakeholders to explore results interactively and contribute to the interpretation.

### 6.2.2 Example of a Shiny app

The Shiny script below reproduces the interactive application used in this chapter. It allows users to visualize the graphs developed above for a selected season and to explore three domains: pest damage, crop yield, and farmers’ perceptions.

A standalone R script for the app is available at the following link: https://github.com/FBaudron/Supporting-codesign/blob/main/06-evaluating-experiments/app.R.

```r
library(openxlsx)
library(ggplot2)
library(ggdist)
library(ggthemes)
library(egg)
library(ggradar)
library(shiny)
library bslib)
library(shinycssloaders)

data = read.xlsx("DATA FOR SHINY.xlsx", sheet = 1)

data$Treatment = factor(data$Treatment, levels = c("CONV", "CA",
"PPULL"))

data_2022_23 = subset(data, Season == "2022/23")
data_2023_24 = subset(data, Season == "2023/24")

# user interface

ui = fluidPage(
  theme = bs_theme(version = 5, bootswatch = "flatly"),
  titlePanel(
    h2("AEI-Zimbabwe On-Farm Experiments Evaluation App (Mbire)",
    style = "font-size: 24px; font-weight: bold;")
  ),
```
```r
 sidebarLayout(
   sidebarPanel(
     selectInput(inputId = "dataset",
                 label = "Choose season:",
                 choices = c("2022/23", "2023/24"))
   ),
   mainPanel(
     tabsetPanel(type = "tabs",
                 tabPanel("Pest damage",
                   withSpinner(plotOutput("plot1")))
                 tabPanel("Crop yield",
                   withSpinner(p
```
```r
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("#FCA50AFF", "#DD513AFF",
  "#932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Armored cricket - Incidence (%)") +
  theme(legendDash = 2, position = "none", 
        plotTitle = element_text(size = 16, face = "bold", margin =
        ggplot2::margin(0, 0, 10, 0)),
  axisText = element_text(size = 10),
  axisTitle = elementalter(),
  legendTitle = elementalter()) +
  coordFlip()

p2 = ggplot(faw, aes(x = Treatment, y = Incidence, fill = Treatment))
+ 
  ggdist::stat_halframe(adjust = 0.6, width = 0.4, width = 0,
  justification = -.4, pointShape = NA, alpha = 0.5, color = "black")
+ 
  geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("#FCA50AFF", "#DD513AFF",
  "#932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Fall armyworm - Incidence (%)") +
  theme(legendDash = 2, position = "none", 
        plotTitle = element_text(size = 16, face = "bold", margin =
        ggplot2::margin(0, 0, 10, 0)),
  axisText = element_text(size = 10),
  axisTitle = elementalter(),
  legendTitle = elementalter()) +
  coordFlip()

p3 = ggplot(stm, aes(x = Treatment, y = Incidence, fill = Treatment))
+ 
  ggdist::stat_halframe(adjust = 0.6, width = 0.4, width = 0,
  justification = -.4, pointShape = NA, alpha = 0.5, color = "black")
+ 
  geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("#FCA50AFF", "#DD513AFF",
  "#932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Stemborer - Incidence (%)") +
```
```r
theme(legendposition = "none",
      plot.title = element_text(size = 16, face = "bold", margin =
  ggplot2::margin(0, 0, 10, 0)),
      axis.text = element_text(size = 10),
      axis.title = elementBlanks,
      legend.title = elementBlanks()) +
      coordflip()

p4 = ggplot(cri, aes(x = Treatment, y = Severity, fill = Treatment))
+ ggdist::statHalfeye(adjust = 0.6, width = 0.4, .width = 0,
justification = -.4, pointColour = NA, alpha = 0.5, color = "black")
+ geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position Jitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("FCA50AFF", "DD513AFF",
"932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Armored cricket - Severity (1-9)") +
  theme(legend position = "none",
      plot.title = element_text(size = 16, face = "bold", margin =
  ggplot2::margin(0, 0, 10, 0)),
      axis.text = element_text(size = 10),
      axis.title = elementBlanks,
      legend.title = elementBlanks()) +
      coordflip()

p5 = ggplot(faw, aes(x = Treatment, y = Severity, fill = Treatment))
+ ggdist::statHalfeye(adjust = 0.6, width = 0.4, .width = 0,
justification = -.4, pointColour = NA, alpha = 0.5, color = "black")
+ geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position Jitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("FCA50AFF", "DD513AFF",
"932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Fall armyworm - Severity (1-9)") +
  theme(legend position = "none",
      plot.title = element_text(size = 16, face = "bold", margin =
  ggplot2::margin(0, 0, 10, 0)),
      axis.text = element_text(size = 10),
      axis.title = elementBlanks,
      legend.title = elementBlanks()) +
      coordflip()
```
```r
p6 = ggplot(stm, aes(x = Treatment, y = Severity, fill = Treatment))

+ ggdist::stat_halfeye(adjust = 0.6, width = 0.4, .width = 0,
justification = -.4, pointShape = NA, alpha = 0.5, color = "black")

+ geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  position_jitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("#FCA50AFF", "#DD513AFF",
"#932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Stemborer - Severity (1-9)") +
  theme(legendDash = 1, width = 1, alpha = 0.5) +
  plot.title = element_text(size = 16, face = "bold", margin =
  ggplot2::margin(0, 0, 10, 0),
  axis.text = element_text(size = 10),
  axis.title = elementekte,
  legend.title = elementblanks)) +
  coordflip()

ggarrange(p1, p4, p2, p5, p3, p6,
  ncol = 2, nrow = 3, widths = c(1, 1), heights = c(1,
  1, 1))

}, height = 480)
```

```r
output$plot2 = renderPlot({
  data = datasetInput()

  yldcrl = data[, c(1:3, 10:11)]
  yldcrl = unique(yldcrl)

  p7 = ggplot(yldcrl, aes(x = Treatment, y = YLD, fill = Treatment)) +
    ggdist::stat_halfeye(adjust = 0.6, width = 0.4, .width = 0,
    justification = -.4, pointShape = NA, alpha = 0.5, color = "black")

  + geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
    geom_point shape = 21, size = 2, alpha = 0.5, position =
    position_jitter(seed = 1, width = .1)) +
    scale_fill_manual(values = c("#FCA50AFF", "#DD513AFF",
    "#932667FF")) +
    ylab("") + xlab("") +
    themeFew() +
    ggtitle("Cereal yield - Grain (t/ha)") +
    theme(legendDash = 1, width = 1, alpha = 0.5) +
    plot.title = element_text(size = 16, face = "bold", margin =
    ggplot2::margin(0, 0, 10, 0)),
    axis.text = element_text(size = 10),
```
```r
axis.title = elementBlanks()
legend.title = elementBlank() +
coordFlip()

p8 = ggplot(yldcr1, aes(x = Treatment, y = BIO, fill = Treatment)) +
  ggdist::statHalfeye(adjust = 0.6, width = 0.4, .width = 0,
  justification = -.4, pointColour = NA, alpha = 0.5, color = "black")

+
  geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.5) +
  geom_point shape = 21, size = 2, alpha = 0.5, position =
  positionJitter(seed = 1, width = .1)) +
  scale_fill_manual(values = c("FCA50AFF", "DD513AFF",
  "932667FF")) +
  ylab("") + xlab("") +
  themeFew() +
  ggtitle("Cereal yield - Biomass (t/ha)") +
  theme(legendDash = "none",
  plotTitle = elementText(size = 16, face = "bold", margin =
  ggplot2::margin(0, 0, 10, 0)),
  axisText = elementText(size = 10),
  axisTitle = elementBlank(),
  legendTitle = elementBlank()) +
coordFlip()

ggarrange(p7, p8,
  ncol = 2, nrow = 1, widths = c(1, 1), heights =
  c(1))

}, height = 480)
```

```r
output$plot3 = renderPlot(
  data = datasetInput()

  of = data[, c(1:2, 12:18)]
  of$Group = rep("Older females", nrow(of))
  of = unique(of)

  om = data[, c(1:2, 19:25)]
  om$Group = rep("Older males", nrow(om))
  om = unique(om)

  yf = data[, c(1:2, 26:32)]
  yf$Group = rep("Young females", nrow(yf))
  yf = unique(yf)

  ym = data[, c(1:2, 33:39)]
  ym$Group = rep("Young males", nrow(ym))
  ym = unique(ym)
```
```r
legendITA = "bottom") +
  theme plot = element_text (size = 18, face = "bold", hjust =
    0.5))

p11 = ggradar(
  ppull[, c(10, 3:9)],
  plot.title = "Push-Pull",
  values.radar = c("0", "5", "10"),
  grid.min = 0, grid.mid = 5, grid.max = 10,
  group line. width = 1,
  group. point. size = 3,
  group. colours = c("#FCA50AFF", "#DD513AFF", "#932667FF",
    "#420A68FF"),
  fill = TRUE, fill. alpha = 0.2,
  background. circle. colour = "white",
  gridline. mid. colour = "grey",
  legend. position = "none") +
  theme (plot. title = element_text (size = 18, face = "bold", hjust =
    0.5))

ggarrange(p9, p10, p11,
  ncol = 3, nrow = 1, widths = c(1, 1, 1),
  heights = c(1))

}, height = 480)
}

# Run the app

shinyApp (ui, server)
```

The deployed app is accessible here: https://fredericbaudron.shinyapps.io/onFarm_experiments_evalua tion_app/.
# Chapter 7

## Scaling

In this chapter, we assess the transferability potential of technologies generated in a specific location — the section of Mbire district where the CGIAR Agroecology Initiative was implemented. Areas sharing similar farming system, environmental, and socio-economic characteristics with the project site are assumed to have high potential for technology transfer, following the approach of Václavík et al. (2016).

We use elevation, soil texture, mean temperature, precipitation, population density, and travel time to the nearest city of at least 50,000 inhabitants as predictors. All data are downloaded via the R package **geodata** (Hijmans et al., 2024).

A standalone R script containing the full workflow is available at the following link: https://github.com/FBaudron/Supporting-codesign/blob/main/07-scaling/scaling.R.

### 7.1 Data preparation

We begin by clearing the workspace, loading the required packages, and importing the baseline dataset, which will be used to locate the project site ('project archetype').

```r
rm(list = ls())

if (!require("openxlsx")) install.packages("openxlsx")
if (!require("terra")) install.packages("terra")
if (!require("geodata")) install.packages("geodata")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("ggspatial")) install.packages("ggspatial")
```
```r
if (!require("tidyterra")) install.packages("tidyterra")

library(openxlsex)
library(terra)
library(geodata)
library(ggplot2)
library(ggspatial)
library(tidyterra)

data = read.xlsx("BASELINE DATA.xlsx")
```

## 7.1.1 Administrative shapefiles

We download shapefiles for Zimbabwe at country level (level 0) and province level (level 1) using the *gadm()* function from **geodata**.

```r
zwe0 = gadm(country = "ZWE", level = 0, path = "Geodata")
plot(zwe0)
```

<figure><img src="image_56.png" /></figure>

```r
zwe1 = gadm(country = "ZWE", level = 1, path = "Geodata")
plot(zwe1)
```
<figure><img src="image_57.png" /></figure>

### 7.1.2 Cropland mask

We download a raster of probability of cropland occurrence and convert it to a presence/absence raster. We crop it to the bounding box of Zimbabwe to reduce extent (using the function `crop()` of `terra`) and mask it so that only cells within the national boundary are retained (using the function `mask()` of `terra`).

This raster represents our region of interest and will be used as a mask for all other rasters to exclude non-cropped areas.

```r
cropland = geodata::cropland(source = "WorldCover", year = 2019, path =
  "Geodata")
cropland = ifel(cropland$cropland > 0, 1, NA)
```

```
cropland = crop(cropland, zwe0)
cropland = mask(cropland, zwe0)

plot(zwe0)
plot(cropland, col = "forestgreen", add = T, legend = F)
plot(zwe1, add = T)
```

```r
cropland = crop(cropland, zwe0)
cropland = mask(cropland, zwe0)

plot(zwe0)
plot(cropland, col = "forestgreen", add = T, legend = F)
plot(zwe1, add = T)
```
<figure><img src="image_58.png" /></figure>

7.1.3 Rasters of predictors

We download and process rasters for all predictors, cropping and masking them
to the cropland raster. We start with the raster of elevation.

elevation = geodata::elevation_30s(country = "ZWE", path = "Geodata")
elevation = crop(elevation, cropland)
elevation = mask(elevation, cropland)
names(elevation) = "elevation"

plot(zwe0)
plot(elevation, add = T, legend = T)
plot(zwe1, add = T)
<figure><img src="image_59.png" /></figure>

We compute a raster of weighted average of sand content across soil layers (0–5 cm, 5–15 cm, 15–30 cm), and crop and mask it to the cropland raster.

```r
sand_5 = geodata::soil_af(var = "sand", depth = 5, path = "Geodata")
sand_15 = geodata::soil_af(var = "sand", depth = 15, path = "Geodata")
sand_30 = geodata::soil_af(var = "sand", depth = 30, path = "Geodata")
sand = (sand_5 * 5 + sand_15 * 10 + sand_30 * 15)/(5 + 10 + 15)
sand = terra::crop(sand, cropland)
sand = terra::mask(sand, cropland)
names(sand) = "sand"

plot(zwe0)
plot(sand, add = T, legend = T)
plot(zwe1, add = T)
```
<figure><img src="image_60.png" /></figure>

We compute a raster of mean annual temperatures, and crop and mask it to the cropland raster.

```r
tavg = worldclim_country(country = "ZWE", var = "tavg", path = "Geodata",
                          res = 0.5)
tavg = mean(tavg)
tavg = terra::crop(tavg, cropland)
tavg = terra::mask(tavg, cropland)
names(tavg) = "tavg"

plot(zwe0)
plot(tavg, add = T, legend = T)
plot(zwe1, add = T)
```
<figure><img src="image_61.png" /></figure>

We compute a raster of annual precipitation, and crop and mask it to the
cropland raster.

prec = worldclim_country(country = "ZWE", var = "prec", path = "Geodata",
  res = 0.5)
prec = sum(prec)
prec = terra::crop(prec, cropland)
prec = terra::mask(prec, cropland)
names(prec) = "prec"

plot(zwe0)
plot(prec, add = T, legend = T)
plot(zwe1, add = T)
<figure><img src="image_62.png" /></figure>

We resample the population raster to match the cropland raster resolution and extent for consistent spatial alignment, and crop and mask it to the cropland raster.

pop = geodata::population(year = 2020, path = "Geodata")

pop = resample(pop, cropland)

pop = crop(pop, cropland)

pop = mask(pop, cropland)

names(pop) = "pop"

plot(zwe0)

plot(pop, add = T, legend = T)

plot(zwe1, add = T)
<figure><img src="image_63.png" /></figure>

Finally, we download a raster of travel time to a city of at least 50,000 inhabitants. We crop and mask it to the cropland raster.

```r
travel_time = travel_time(to = "city", size = 6, up = TRUE, path =
  ↔ "Geodata")
travel_time = crop(travel_time, cropland)
travel_time = mask(travel_time, cropland)
names(travel_time) = "travel_time"

plot(zwe0)
plot(travel_time, add = T, legend = T)
plot(zwe1, add = T)
```

<figure><img src="image_64.png" /></figure>
### 7.1.4 Preprocessing

We merge all rasters into a single stack and normalize each variable to zero mean and unit variance.

```
all = c(elevation, sand, tavg, prec, pop, travel_time)

means = as.numeric(global(all, "mean", na.rm = TRUE)[, 1])
sds = as.numeric(global(all, "sd", na.rm = TRUE)[, 1])
all_z = (all - means)/sds
names(all_z) = paste0 names(all), "_n")
```

We convert the baseline data to a spatial object, extract raster values for each surveyed farm, and compute the mean for the project site.

```
data = vect(data, geom =
  c("_Record.your.current.location_longitude",
    "_Record.your.current.location_lat),
  crs = "EPSG:4326")

vals = terra::extract(all_z, data)

mean_site = colMeans vals[, -1], na.rm = TRUE)

mean_site DF = data.frame variable = names(all_z), mean_value =
  mean_site)

mean_site DF
```

<table><thead><tr><th></th><th>variable</th><th>mean_value</th></tr></thead><tbody><tr><td>elevation_n</td><td>elevation_n</td><td>-2.0614019</td></tr><tr><td>sand_n</td><td>sand_n</td><td>0.3985115</td></tr><tr><td>tavg_n</td><td>tavg_n</td><td>1.7957409</td></tr><tr><td>prec_n</td><td>prec_n</td><td>0.2687717</td></tr><tr><td>pop_n</td><td>pop_n</td><td>-0.1931808</td></tr><tr><td>travel_time_n</td><td>travel_time_n</td><td>0.1083564</td></tr></tbody></table>

## 7.2 Transferability potential

### 7.2.1 Distance to project prototype

We compute the Euclidean distance to the project prototype for each cell, following the method of Václavík et al. (2016).

We plot this raster with the centroid of the project site.
```r
dist_to_site = sqrt(sum((all_z - mean_site)^2))

centroid = colMeans(crds(dataVect))

plot(dist_to_site, main = "Euclidean distance to project prototype")
plot(zwei, add = T)
points(centroid[x"], centroid[y"], pch = 21, bg = "orangered", cex = 2)
```

<figure><img src="image_65.png" /></figure>

### 7.2.2 Transferability potential

We derive transferability potential, bounded between 0 and 1, and plot the
resulting raster with the centroid of the project site.

```r
tp_site = 1/(1 + dist_to_site)
names.tp_site = "tp"

plot.tp_site, main = "Transferability potential")
plot(zwei, add = T)
points(centroid[x"], centroid[y"], pch = 21, bg = "orangered", cex = 2)
```
<figure><img src="image_66.png" /></figure>

To facilitate interpretability, we discretize the raster into four classes (low, moderate, high, very high transferability) using quartiles.

```r
tp_site_q = classify.tp_site$tp, c(quantile(values.tp_site$tp, na.rm =
  T))[1], quantile(values.tp_site$tp,
  na.rm = T))[2], quantile(values.tp_site$tp, na.rm = T))[3],
  quantile(values.tp_site$tp,
  na.rm = T))[4], quantile(values.tp_site$tp, na.rm = T))[5]))
```

```r
plot.tp_site_q, type = "classes", levels = c("Low", "Moderate", "High",
  "Very high"),
  main = "Transferability potential")
plot(zwe1, add = T)
points(centroid["x"], centroid["y"], pch = 21, bg = "orangered", cex = 2)
```
<figure><img src="image_67.png" /></figure>

### 7.2.3 Publication-ready map

Below, we propose a more publication-ready version of the map using the R packages `ggplot2` (Wickham, 2016), `tidyterra` (Hernangómez, 2023), and `ggspatial` (Dunnington, 2023).

We convert the raster (and the centroid) to a dataframe, assign ordered factor labels, and plot the polished map.

```r
tp_site_q DF = as.data.frame.tp_site_q, xy = TRUE)

names.tp_site_q DF)[3] = "class_id"

tp_site_q DF$class_id_num = as.numeric.tp_site_q DF$class_id)

tp_site_q DF$class_label = factor.tp_site_q DF$class_id_num, levels =
  1:4, labels = c("Low",
  "Moderate", "High", "Very high"), ordered = TRUE)

centroid DF = data.frame(x = centroid[x], y = centroid[y])

ggplot() + themeBW() + geom_raster(data = na.omit.tp_site_q DF), aes(x =
  x, y = y,
  fill = class_label)) + geom_spatvector(data = zwe1, fill = NA,
  linewidth = 1,
  color = "black") + geom_point(data = centroid DF, aes(x = x, y = y),
  color = "black",
  fill = "tomato", shape = 21, size = 5, stroke = 2) +
  scale_fill_manual(values = c(Low = "#420A68FF",
  Moderate = "#932667FF", High = "#DD513AFF", `Very high` =
  "#FCA50AFF")) + labs(fill = "Transferability") +
```
```r
theme plot.title = element CAT(), axis title = element CAT(),
  axis.text.y = element_text(angle = 90,
    hjust = 0.5, size = 12), axis.text.x = element_text(size = 12),
  legend.title = element_text(size = 14),
  legend.text = element_text(size = 12), legend的位置 = "right",
  plot-margin = ggplot2::margin(10,
    10, 10, 10)) + annotation_scale(location = "bl", width_hint =
    0.5, height = unit(0.25,
      "cm"), text_cex = 1, pad_x = unit(0.15, "cm"), pad_y = unit(0.15,
      "cm")) + annotation_northArrow(location = "tl",
      which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"),
      pad_x = unit(0.15,
        "cm"), pad_y = unit(0.15, "cm"), style =
      northArrow_fancyOrienteering(text_size = 10))
```
<figure><img src="image_68.png" /></figure>

### 7.2.4 Area and population by transferability domain

```
area_tbl = zonal(cellSize.tp_site_q, unit = "km"), tp_site_q, "sum",
  na.rm = TRUE)
colnames(area_tbl) = c("Class", "Area_km2")
pop_tbl = zonal(pop, tp_site_q, fun = "sum", na.rm = TRUE)
```
colnames(pop_tbl) = c("Class", "Population")

result = merge(area_tbl, pop_tbl, by = "Class")

result$Class = c("Low", "Moderate", "High", "Very high")

result

<table><thead><tr><th></th><th>Class</th><th>Area_km2</th><th>Population</th></tr></thead><tbody><tr><td>1</td><td>Low</td><td>64435.69</td><td>7306546</td></tr><tr><td>2</td><td>Moderate</td><td>64312.90</td><td>2525554</td></tr><tr><td>3</td><td>High</td><td>64213.90</td><td>2242085</td></tr><tr><td>4</td><td>Very high</td><td>64444.78</td><td>2435429</td></tr></tbody></table>
# Chapter 8

# Conclusions

This book demonstrates that data can enrich every step of the codesign process, from the diagnostic stage, to the testing and validation stage, and the evaluation stage. Although related scripts have been presented sequentially, the process should be seen as iterative, rather than linear, as described for example in the DEED (Describe, Explain, Explore, Design) cycle (Giller et al., 2008; Falconnier et al., 2017).

Negotiation is at the center of this process. The role of quantitative analysis (e.g., methods presented in this book) is to make on-going negotiations more explicit, creative and equitable. This implies a nonpartisan posture from researchers towards particular technologies and implies that the main function of living labs may be to serve as negotiation platforms (Moruzzo et al., 2024).

We hope this book can help scientists and practitioners make more use of data and data analysis in their codesign work. Additional analytics/chapters could be added in a next edition. For example, the selection of project site(s) could be guided by the delineation of spatial domains (Muthoni et al., 2017). The selection of innovations to be trialed could be informed by the identification of determinants of crop and livestock productivity (through e.g., stochastic frontier analysis) (Baudron et al., 2024) and/or the analysis of resource flows (nutrients, energy, labour, etc) within representative farms (Chavarría et al., 2018). These innovations could also include practices used by ‘positive deviants’ (Mutsamba-Magwaza et al., 2025).
Acknowledgments

This study was supported by the CGIAR initiative 'Transformational Agroecology across Food, Land, and Water Systems' and the CGIAR 'Multifunctional Landscapes' program
Bibliography

Allaire, J., Xie, Y., Dervieux, C., McPherson, J., Luraschi, J., Ushey, K., Atkins, A., Wickham, H., Cheng, J., Chang, W., and Iannone, R. (2024). *rmarkdown: Dynamic Documents for R*. R package version 2.29.

Alvarez, S., Timler, C. J., Michalscheck, M., Paas, W., Descheemaeker, K., Tittonell, P., Andersson, J. A., and Groot, J. C. (2018). Capturing farm diversity with hypothesis-based typologies: An innovative methodological framework for farming system typology development. *PLoS ONE*, 13.

Asten, P. J. V., Kaaria, S., Fermont, A. M., and Delve, R. J. (2009). Challenges and lessons when using farmer knowledge in agricultural research and development projects in Africa. *Experimental Agriculture*, 45:1–14.

Attali, D. and Sali, A. (2024). *shinycssloaders: Add Loading Animations to a 'shiny' Output While It's Recalculating*. R package version 1.1.0.

Baudron, F. (2024). *Supporting biodiversity-smart smallholder farming in Africa through systems agronomy*.

Baudron, F., Corbeels, M., Andersson, J. A., Sibanda, M., and Giller, K. E. (2011). Delineating the drivers of waning wildlife habitat: The predominance of cotton farming on the fringe of protected areas in the mid-zambezi valley, Zimbabwe. *Biological Conservation*, 144:1481–1493.

Baudron, F., Guerrini, L., Chimimba, E., Chimusimbe, E., and Giller, K. E. (2022). Commodity crops in biodiversity-rich production landscapes: Friends or foes? The example of cotton in the mid zambezi valley, Zimbabwe. *Biological Conservation*, 267.

Baudron, F., Tui, S. H. K., Silva, J. V., Chakoma, I., Matangi, D., Nyagumbo, I., and Dube, S. (2024). Tailoring interventions through a combination of statistical typology and frontier analysis: a study of mixed crop-livestock farms in semi-arid Zimbabwe. *Experimental Agriculture*, 60.

Bion, R. (2023). *ggradar: Create radar charts using ggplot2*. R package version 0.2.
Chang, W., Cheng, J., Allaire, J., Sievert, C., Schloerke, B., Xie, Y., Allen, J., McPherson, J., Dipert, A., and Borges, B. (2024). *shiny: Web Application Framework for R*. R package version 1.10.0.

Charrad, M., Ghazzali, N., Boiteau, V., and Niknafs, A. (2014). NbClust: An R package for determining the relevant number of clusters in a data set. *Journal of Statistical Software*, 61(6):1–36.

Chavarría, J. Y. D., Baudron, F., and Sunderland, T. (2018). Retaining forests within agricultural landscapes as a pathway to sustainable intensification: Evidence from southern ethiopia. *Agriculture, Ecosystems and Environment*, 263:41–52.

Coid, C., Gaidet-Drapier, N., Moya, C., Poilecot, P., Poulet, D., Renaud, P. C., Ricard, X., and Takawira, S. (2001). *Les hommes et les animaux dans la moyenne vallée du Zambèze, Zimbabwe*. Cirad edition.

Dunnington, D. (2023). *ggspatial: Spatial Data Framework for ggplot2*. R package version 1.1.9.

Falconnier, G. N., Descheemaeker, K., Mourik, T. A. V., Adam, M., Sogoba, B., and Giller, K. E. (2017). Co-learning cycles to support the design of innovative farm systems in southern mali. *European Journal of Agronomy*, 89:61–74.

Firke, S. (2024). *janitor: Simple Tools for Examining and Cleaning Dirty Data*. R package version 2.2.1.

Gaidet, N., Fritz, H., and Nyahuma, C. (2003). A participatory counting method to monitor populations of large mammals in non-protected areas: A case study of bicycle counts in the zambezi valley, zimbabwe. *Biodiversity and Conservation*, 12:1571–1585.

Galili, T. (2015). *dendextend: an r package for visualizing, adjusting, and comparing trees of hierarchical clustering*. *Bioinformatics*.

Giller, K. E., Leeuwis, C., Andersson, J. A., Andriesse, W., Brouwer, A., Frost, P. G. H., Hebinck, P., Heitk??nig, I., van Ittersum, M. K., Koning, N., Ruben, R., Slingerland, M., Udo, H., Veldkamp, T. A., de Vijver, C. V., van Wijk, M. T., and Windmeijer, P. N. (2008). Competing claims on natural resources: What role for science? *Ecology and Society*, 13.

Hammond, J., Rosenblum, N., Breseman, D., Gorman, L., Manners, R., van Wijk, M. T., Sibomana, M., Remans, R., Vanlauwe, B., and Schut, M. (2020). Towards actionable farm typologies: Scaling adoption of agricultural inputs in rwanda. *Agricultural Systems*, 183:102857.

Hassall, K. L., Baudron, F., MacLaren, C., Cairns, J. E., Ndhlela, T., McGrath, S. P., Nyagumbo, I., and Haefele, S. M. (2023). Construction of a generalised farm typology to aid selection, targeting and scaling of onfarm research. *Computers and Electronics in Agriculture*, 212.
Hernangómez, D. (2023). Using the tidyverse with terra objects: the tidyterra package. *Journal of Open Source Software*, 8(91):5751.

Hijmans, R. J. (2024). *terra: Spatial Data Analysis*. R package version 1.7-78.

Hijmans, R. J., Barbosa, M., Ghosh, A., and Mandel, A. (2024). *geodata: Download Geographic Data*. R package version 0.6-2.

Hossain, M., Leminen, S., and Westerlund, M. (2019). A systematic review of living lab literature. *Journal of Cleaner Production*, 213:976–988.

Houérou, H. N. L. and Hoste, C. H. (1977). Rangeland production and annual rainfall relations in the mediterranean basin and in the african sahelosudanian zone. *Journal of Range Management*, 30:181–189.

Hulme, D. and Murphree, M. W. (1999). Communities, wildlife and the ‘new conservation’ in Africa. *Journal of International Development*, 11:277–285.

Jahnke, H. E. (1982). *Livestock Production Systems and Livestock Development in Tropical Africa*, volume 65. Kieler wiss edition.

Kay, M. (2025). *ggdist: Visualizations of Distributions and Uncertainty*. R package version 3.3.3.

Kebede, B. (2009). Community wealth ranking and household surveys: An integrative approach. *Journal of Development Studies*, 45:1731–1746.

Kennedy, G., Berardo, A., Papavero, C., Horjus, P., Ballard, T., Dop, M., Delbaere, J., and Brouwer, I. D. (2010). Proxy measures of household food consumption for food security assessment and surveillance: comparison of the household dietary diversity and food consumption scores. *Public health nutrition*, 13:2010–2018.

Liaw, A. and Wiener, M. (2002). Classification and regression by randomforest. *R News*, 2(3):18–22.

Lüdecke, D., Ben-Shachar, M. S., Patil, I., Waggoner, P., and Makowski, D. (2021). performance: An R package for assessment, comparison and testing of statistical models. *Journal of Open Source Software*, 6(60):3139.

Meynard, J. M., Dedieu, B., and Bos, A. P. B. (2012). *Re-design and co-design of farming systems. An overview of methods and practices*, pages 405–429. Springer Netherlands.

Moruzzo, R., Diaz, S. E., Granai, G., Iacovo, F. D., and Riccioli, F. (2024). Living lab as support for co-creation of value: application to agro-biodiversity contracting solutions. *Local Environment*, pages 1–14.

Muthoni, F. K., Guo, Z., Bekunda, M., Sseguya, H., Kizito, F., Baijukya, F., and Hoeschle-Zeledon, I. (2017). Sustainable recommendation domains for scaling agricultural technologies in tanzania. *Land Use Policy*, 66:34–48.
Mutsamba-Magwaza, E. F., Baudron, F., Franke, A. C., Watt, E. V. D., and Nyagumbo, I. (2025). Learning from positive deviants' practices to improve the performance of mixed crop-livestock systems in zimbabwe. *Agricultural Systems*, 228.

Oksanen, J., Simpson, G. L., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O'Hara, R., Solymos, P., Stevens, M. H. H., Szoecs, E., Wagner, H., Barbour, M., Bedward, M., Bolker, B., Borcard, D., Borman, T., Carvalho, G., Chirico, M., De Caceres, M., Durand, S., Evangelista, H. B. A., FitzJohn, R., Friendly, M., Furneaux, B., Hannigan, G., Hill, M. O., Lahti, L., Martino, C., McGlinn, D., Ouellette, M.-H., Ribeiro Cunha, E., Smith, T., Stier, A., Ter Braak, C. J., and Weedon, J. (2025). *vegan: Community Ecology Package*. R package version 2.7-1.

Sievert, C., Cheng, J., and Aden-Buie, G. (2025). *bslib: Custom 'Bootstrap' 'Sass' Themes for 'shiny' and 'r markdown'.* R package version 0.9.0.

Sjoberg, D. D., Whiting, K., Curry, M., Lavery, J. A., and Larmarange, J. (2021). Reproducible summary tables with the gtsummary package. *The R Journal*, 13:570–580.

Taylor, R. (2009). Community based natural resource management in zimbabwe: the experience of campfire. *Biodiversity and Conservation*, 18:2563–2583.

Thioulouse, J., Dray, S., Dufour, A.-B., Siberchicot, A., Jombart, T., and Pavoine, S. (2018). *Multivariate Analysis of Ecological Data with ade4*. Springer.

Tittonell, P., Muriuki, A., Shepherd, K. D., Mugendi, D., Kaizzi, K. C., Okeyo, J., Verchot, L., Coe, R., and Vanlauwe, B. (2010). The diversity of rural livelihoods and their influence on soil fertility in agricultural systems of east Africa - a typology of smallholder farms. *Agricultural Systems*, 103:83–97.

Václavík, T., Langerwisch, F., Cotter, M., Fick, J., Häuser, I., Hotes, S., Kamp, J., Settele, J., Spangenberg, J. H., and Seppelt, R. (2016). Investigating potential transferability of place-based research in land system science. *Environmental Research Letters*, 11.

Wickham, H. (2016). *ggplot2: Elegant Graphics for Data Analysis*. Springer-Verlag New York.

Wickham, H., Vaughan, D., and Girlich, M. (2024). *tidyr: Tidy Messy Data*. R package version 1.3.1.

Wu, Y. and Colautti, R. (2022). *baRcodeR: Label Creation for Tracking and Collecting Data from Biological Samples*. R package version 0.1.7.

Xie, Y. (2025). *bookdown: Authoring Books and Technical Documents with R Markdown*. R package version 0.44.1.
Zhu, H. (2024). *kableExtra: Construct Complex Table with 'kable' and Pipe Syntax*. R package version 1.4.0.