Model for Demo Farm Treatment Analysis

[Shamsudini Amidu]

Dear Family,

We recently conducted an experiment to promote environmentally sustainable input bundles (certified seeds, inoculants, Yara legume fertilizer/TSP, organic farming aid (OFA), and biochar) aimed at increasing soybean and groundnut yields. The study used 97 demonstration farms (51 groundnut and 46 soybean), each serving as a platform to promote adoption among farmers.

Each demonstration farm had six treatment plots. For both groundnut and soybean, the treatments were:
(1) Certified seed only (control),
(2) Certified seed + inoculant,
(3) Certified seed + inoculant + OFA,
(4) Certified seed + TSP/Yara legume fertilizer,
(5) Certified seed + inoculant + biochar, and
(6) Certified seed + OFA.

For every plot, we collected fresh and dry weight yield data, resulting in 582 observations across the 97 demo farms (six treatments per farm).

The main objective of the analysis is to identify which input bundle produces the highest yield for soybeans and groundnuts. Given that each demonstration farm contributes multiple observations, the data are clustered at the farm level and observations are not independent.

To account for this structure, I estimated a mixed-effects model as follows:

mixed DriedYield i.trt|| FarmerID: if Crop=="Goundnut demo farms"
testparm i.trt
margins trt

mixed DriedYield i.trt|| FarmerID: if Crop=="soybean demo farms" testparm i.trt margins trt While this approach accounts for clustering by demo farm, I would appreciate your guidance on whether this is the most appropriate model specification, including the correct treatment of fixed and random effects, and the best way to compare treatment effects to identify the most effective input bundle.

I am happy to share a sample dataset or the do-file if helpful.

* Example generated by -dataex-. For more info, type help dataex
clear
input long FarmerID str19 Crop byte(plotlayout_index trt) double(FreshYield DriedYield)
 1 "soybean demo farms"  1 1   5.7  14.6
 1 "soybean demo farms"  2 2   4.8   6.8
 1 "soybean demo farms"  3 3   6.8  14.6
 1 "soybean demo farms"  4 4   4.7   5.7
 1 "soybean demo farms"  5 5     6  14.6
 1 "soybean demo farms"  6 6  12.9   7.8
 2 "soybean demo farms"  1 1  11.5 13.25
 2 "soybean demo farms"  2 2 12.65  19.9
 2 "soybean demo farms"  3 3 12.55 10.55
 2 "soybean demo farms"  4 4 14.25 10.65
 2 "soybean demo farms"  5 5  21.1 16.35
 2 "soybean demo farms"  6 6 17.53   9.5
 3 "soybean demo farms"  1 1  16.8 14.45
 3 "soybean demo farms"  2 2  11.2  11.2
 3 "soybean demo farms"  3 3  10.6  16.8
 3 "soybean demo farms"  4 4 20.75 20.75
 3 "soybean demo farms"  5 5 14.45  10.6
 3 "soybean demo farms"  6 6 15.71 15.71
 4 "soybean demo farms"  1 1  11.3  8.47
 4 "soybean demo farms"  2 2  7.18  8.11
 4 "soybean demo farms"  3 3 10.12  6.01
 4 "soybean demo farms"  4 4  8.08 10.53
 4 "soybean demo farms"  5 5 12.62  6.46
 4 "soybean demo farms"  6 6 10.76  9.14
 5 "soybean demo farms"  1 1 13.52 12.47
 5 "soybean demo farms"  2 2 10.12 12.29
 5 "soybean demo farms"  3 3  8.22  8.22
 5 "soybean demo farms"  4 4  14.7  7.53
 5 "soybean demo farms"  5 5  9.86     8
 5 "soybean demo farms"  6 6 13.38 11.26
 6 "Goundnut demo farms" 1 1    12   6.4
 6 "Goundnut demo farms" 2 2    11   5.4
 6 "Goundnut demo farms" 3 3    14   5.4
 6 "Goundnut demo farms" 4 4    15   5.3
 6 "Goundnut demo farms" 5 5  14.8     5
 6 "Goundnut demo farms" 6 6    20   5.3
 7 "Goundnut demo farms" 1 1    12  9.73
 7 "Goundnut demo farms" 2 2    14 13.11
 7 "Goundnut demo farms" 3 3    16  15.5
 7 "Goundnut demo farms" 4 4    10 14.87
 7 "Goundnut demo farms" 5 5    18 17.39
 7 "Goundnut demo farms" 6 6    15 13.64
 8 "Goundnut demo farms" 1 1    10  13.5
 8 "Goundnut demo farms" 2 2    14  13.3
 8 "Goundnut demo farms" 3 3    18   8.7
 8 "Goundnut demo farms" 4 4    17   8.2
 8 "Goundnut demo farms" 5 5    14   9.6
 8 "Goundnut demo farms" 6 6    14   6.4
 9 "Goundnut demo farms" 1 1    14     9
 9 "Goundnut demo farms" 2 2  13.2     7
 9 "Goundnut demo farms" 3 3    13   7.5
 9 "Goundnut demo farms" 4 4    18     6
 9 "Goundnut demo farms" 5 5  15.9   7.5
 9 "Goundnut demo farms" 6 6  18.1     8
10 "Goundnut demo farms" 1 1    14    10
10 "Goundnut demo farms" 2 2    12    10
10 "Goundnut demo farms" 3 3    15   8.5
10 "Goundnut demo farms" 4 4    13   8.4
10 "Goundnut demo farms" 5 5    10  12.5
10 "Goundnut demo farms" 6 6    12   9.4
end
label values trt Demo_type
label def Demo_type 1 "Plot 1 (Certified seed only (control)", modify
label def Demo_type 2 "Plot 2 (Certified seed + inoculant )", modify
label def Demo_type 3 "Plot 3 (Certified Seed + inoculant + OFA)", modify
label def Demo_type 4 "Plot4 (Certified Seed + TSP/Yara legume fertilizer)", modify
label def Demo_type 5 "Plot5 (Certified seed + inoculant + Biochar)", modify
label def Demo_type 6 "Plot6 (Certified seed + OFA )", modify

I need your help.

Thank you very much for your support.


Kind regards,

[Carlo Lazzaro]

Shamsudini:
two comments about your post:
1) unless the dataset you shared is an excerpt of the one you're working on, 6 obs are clearly too few for any convincing inference;
2) clustering is not a default option, but should be imposed via -vce(cluster FarmerID).

[Shamsudini Amidu]

Dear Carlo Lazzaro,
Thank you very much for your feedback.
The dataset I previously shared was only a small example for illustration purposes and did not represent the full dataset. I fully agree that the number of observations in that example was too few for any meaningful inference.

I have now provided an updated dataset with 100 observations, which better reflects the data structure and should allow for more appropriate guidance on the analysis.

* Example generated by -dataex-. For more info, type help dataex
clear
input float FarmerID str19 Crop byte(plotlayout_index trt) double(FreshYield DriedYield)
4 "soybean demo farms" 1 1 5.7 14.6
4 "soybean demo farms" 2 2 4.8 6.8
4 "soybean demo farms" 3 3 6.8 14.6
4 "soybean demo farms" 4 4 4.7 5.7
4 "soybean demo farms" 5 5 6 14.6
4 "soybean demo farms" 6 6 12.9 7.8
86 "soybean demo farms" 1 1 11.5 13.25
86 "soybean demo farms" 2 2 12.65 19.9
86 "soybean demo farms" 3 3 12.55 10.55
86 "soybean demo farms" 4 4 14.25 10.65
86 "soybean demo farms" 5 5 21.1 16.35
86 "soybean demo farms" 6 6 17.53 9.5
98 "soybean demo farms" 1 1 16.8 14.45
98 "soybean demo farms" 2 2 11.2 11.2
98 "soybean demo farms" 3 3 10.6 16.8
98 "soybean demo farms" 4 4 20.75 20.75
98 "soybean demo farms" 5 5 14.45 10.6
98 "soybean demo farms" 6 6 15.71 15.71
490 "soybean demo farms" 1 1 11.3 8.47
490 "soybean demo farms" 2 2 7.18 8.11
490 "soybean demo farms" 3 3 10.12 6.01
490 "soybean demo farms" 4 4 8.08 10.53
490 "soybean demo farms" 5 5 12.62 6.46
490 "soybean demo farms" 6 6 10.76 9.14
689 "soybean demo farms" 1 1 13.52 12.47
689 "soybean demo farms" 2 2 10.12 12.29
689 "soybean demo farms" 3 3 8.22 8.22
689 "soybean demo farms" 4 4 14.7 7.53
689 "soybean demo farms" 5 5 9.86 8
689 "soybean demo farms" 6 6 13.38 11.26
975 "soybean demo farms" 1 1 4.8 4.8
975 "soybean demo farms" 2 2 5.2 5.6
975 "soybean demo farms" 3 3 5.9 5.7
975 "soybean demo farms" 4 4 7.1 7
975 "soybean demo farms" 5 5 0 5.1
975 "soybean demo farms" 6 6 5.7 0
488 "soybean demo farms" 1 1 6 7.3
488 "soybean demo farms" 2 2 6.2 5.9
488 "soybean demo farms" 3 3 7.8 0
488 "soybean demo farms" 4 4 8.2 7.5
488 "soybean demo farms" 5 5 0 8.1
488 "soybean demo farms" 6 6 7.4 6.1
497 "soybean demo farms" 1 1 6.4 0
497 "soybean demo farms" 2 2 7.1 6.3
497 "soybean demo farms" 3 3 7.3 6.5
497 "soybean demo farms" 4 4 11.5 7.1
497 "soybean demo farms" 5 5 0 11.5
497 "soybean demo farms" 6 6 6.5 7.2
929 "soybean demo farms" 1 1 6.5 6.6
929 "soybean demo farms" 2 2 6.8 6.4
929 "soybean demo farms" 3 3 7 6.9
929 "soybean demo farms" 4 4 8.9 8.7
929 "soybean demo farms" 5 5 0 6.4
929 "soybean demo farms" 6 6 6.5 0
367 "soybean demo farms" 1 1 8 0
367 "soybean demo farms" 2 2 12.5 12.2
367 "soybean demo farms" 3 3 13 8
367 "soybean demo farms" 4 4 24 23.5
367 "soybean demo farms" 5 5 0 13
367 "soybean demo farms" 6 6 9 8.9
991 "soybean demo farms" 1 1 5.8 6.8
991 "soybean demo farms" 2 2 6.4 8.2
991 "soybean demo farms" 3 3 8.3 8.7
991 "soybean demo farms" 4 4 8.9 5.6
991 "soybean demo farms" 5 5 0 0
991 "soybean demo farms" 6 6 7 6.2
18 "soybean demo farms" 1 1 5.6 5.8
18 "soybean demo farms" 2 2 5.95 5.1
18 "soybean demo farms" 3 3 5.1 5.6
18 "soybean demo farms" 4 4 5.1 5.6
18 "soybean demo farms" 5 5 0 0
18 "soybean demo farms" 6 6 5.9 5.8
234 "Goundnut demo farms" 1 1 12 6.4
234 "Goundnut demo farms" 2 2 11 5.4
234 "Goundnut demo farms" 3 3 14 5.4
234 "Goundnut demo farms" 4 4 15 5.3
234 "Goundnut demo farms" 5 5 14.8 5
234 "Goundnut demo farms" 6 6 20 5.3
451 "Goundnut demo farms" 1 1 12 9.73
451 "Goundnut demo farms" 2 2 14 13.11
451 "Goundnut demo farms" 3 3 16 15.5
451 "Goundnut demo farms" 4 4 10 14.87
451 "Goundnut demo farms" 5 5 18 17.39
451 "Goundnut demo farms" 6 6 15 13.64
609 "Goundnut demo farms" 1 1 10 13.5
609 "Goundnut demo farms" 2 2 14 13.3
609 "Goundnut demo farms" 3 3 18 8.7
609 "Goundnut demo farms" 4 4 17 8.2
609 "Goundnut demo farms" 5 5 14 9.6
609 "Goundnut demo farms" 6 6 14 6.4
853 "Goundnut demo farms" 1 1 14 9
853 "Goundnut demo farms" 2 2 13.2 7
853 "Goundnut demo farms" 3 3 13 7.5
853 "Goundnut demo farms" 4 4 18 6
853 "Goundnut demo farms" 5 5 15.9 7.5
853 "Goundnut demo farms" 6 6 18.1 8
131 "Goundnut demo farms" 1 1 14 10
131 "Goundnut demo farms" 2 2 12 10
131 "Goundnut demo farms" 3 3 15 8.5
131 "Goundnut demo farms" 4 4 13 8.4
end
label values trt Demo_type
label def Demo_type 1 "Plot 1 (Certified seed only (control)", modify
label def Demo_type 2 "Plot 2 (Certified seed + inoculant )", modify
label def Demo_type 3 "Plot 3 (Certified Seed + inoculant + OFA)", modify
label def Demo_type 4 "Plot4 (Certified Seed + TSP/Yara legume fertilizer)", modify
label def Demo_type 5 "Plot5 (Certified seed + inoculant + Biochar)", modify
label def Demo_type 6 "Plot6 (Certified seed + OFA )", modify


Thank you once again for your time and support, and I look forward to your support.

[Carlo Lazzaro]

Shamsudini:
as your -mixed- model does not have a time-series dimension, why not taking a look at example #15 -anova- entry, Stata .pdf manual?

. anova DriedYield FarmerID trt if Crop=="Goundnut demo farms", repeated( trt )

                         Number of obs =         28    R-squared     =  0.7760
                         Root MSE      =    2.02093    Adj R-squared =  0.6640

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |  254.72438          9   28.302709      6.93  0.0003
                         |
                FarmerID |   244.8799          4   61.219975     14.99  0.0000
                     trt |  9.8652396          5   1.9730479      0.48  0.7843
                         |
                Residual |  73.514594         18   4.0841441  
              -----------+----------------------------------------------------
                   Total |  328.23897         27   12.156999  


Between-subjects error term:  FarmerID
                     Levels:  5         (4 df)
     Lowest b.s.e. variable:  FarmerID

Repeated variable: trt
                                          Huynh-Feldt epsilon        =  0.4069
                                          Greenhouse-Geisser epsilon =  0.2882
                                          Box's conservative epsilon =  0.2000

                                            ------------ Prob > F ------------
                  Source |     df      F    Regular    H-F      G-G      Box
              -----------+----------------------------------------------------
                     trt |      5     0.48   0.7843   0.6382   0.5843   0.5293
                Residual |     18
              ----------------------------------------------------------------

.

[Joseph Coveney]

The main objective of the analysis is to identify which input bundle produces the highest yield for soybeans and groundnuts.

In addition to Carlo's excellent suggestion, you could also consider fitting an omnibus regression model with both crop types. Such an approach might be more efficient in that you could include all of the farms at once instead of two subsets of farms each with a smaller ability to contribute to the precision of the estimates. You would want to estimate the residual variance separately by crop, which is possible using mixed to fit the ANOVA model.

I illustrate below using the union of your two data snippets, after correcting what appears to be a duplicate farmer ID and converting the string variable for crop identifier to numeric. (I also renamed the variables to shorter names for convenience.)

version 19

clear *

quietly input long FarmerID str19 Crop byte(plotlayout_index trt) ///
    double(FreshYield DriedYield)
<redacted for brevity>

// Duplicate FarmerID fixup
generate long row = _n
summarize row if FarmerID == 10, meanonly
quietly replace FarmerID = 40 if FarmerID == 4 & row > r(max)
drop row

// Replace cumbersome variable names
rename (FarmerID plotlayout_index FreshYield DriedYield) (fid pli fld dld)
label variable fid "Farmer ID"
label variable pli "Plot-layout Index"
label variable trt "Soil Treatment"
label variable fld "Fresh Yield (100 kg/ha)"
label variable dld "Dried Yield (100 kg/ha)"
sort fid pli

// Crop types
generate byte cid = strpos(strlower(Crop), "soybean") > 0
label variable cid Crop
label define Crops 0 Peanuts 1 Soybeans
label values cid Crops
drop Crop

// Soil treatments
#delimit ;
label define Tx 
    1 None 
    2 Inoculant 
    3 "Inoc. + OFA" 
    4 "TSP + YLF" 
    5 "Inoc. + Biochar" 
    6 OFA;
#delimit cr
label values trt Tx

// Fit omnibus regression model
mixed dld i.cid##i.trt || fid: , ///
    reml dfmethod(kroger) ///
    residuals(independent, by(cid)) ///
    nolrtest nolog

// Get the denominator degrees-of-freedom for -margins- to use
contrast cid trt cid#trt, small
local df = r(df2)[1, 3]

set scheme s2color
local o plotopts(msize(medlarge) mcolor(black) mfcolor(white) connect(none)) ///
    ciopts(lcolor(black)) ylabel( , angle(horizontal) nogrid) ///
    ytitle(`: variable label dld')

margins 0.cid#i.trt, level(68.27) df(`df')
marginsplot , `o' title(Peanuts, ring(0) position(11))
quietly graph export Peanuts.png, replace

margins 1.cid#i.trt, level(68.27) df(`df')
marginsplot , `o' title(Soybeans, ring(0) position(11))
quietly graph export Soybeans.png, replace

exit

I guessed at the units for crop yields for the plots. Unless the data that you show are made up for illustration, you might want to double check those yield values: some of the zeroes seem inconsistent between fresh and dried, for example.

Complete do-file and log file are attached if you're interested further.

[Shamsudini Amidu]

First of all, Carlo and Joseph, thank you very much for your continuous support and guidance on this analysis. I truly appreciate the time and effort you have invested in reviewing the data and providing suggestions.

I would like to provide an update on my progress:

• Initially, I shared only 100 observations (for illustration) due to data sharing constraints, but the full dataset actually contains 582 observations.
• I have now dropped zero-yield plots and am working with 568 observations. This means some farmers do not have complete observations for all six plots, as at least one plot may have zero yield. In such cases, I dropped only the affected plot rather than the entire farmer.

With your permission, I attached the results and the do file for your perusal if you are interested.

However, I need some clarification on a few points:

1. Do you recommend dropping the entire farmer if at least one plot has no harvest, or is it acceptable to drop only the affected plot?
2. Should I check for normality of the dried yield (dld) variable before running the mixed-effects model?
3. I noticed that none of the treatment estimates are statistically significant. How should I interpret this in the context of the analysis?
4. Would it be advisable to transform dld into yield per hectare (instead of using kilograms) to potentially improve the model estimation?

Your guidance on these points will be highly valuable for ensuring the analysis is robust and interpretable.

Thank you again for your support

[Carlo Lazzaro]

Shamsudini:
1) it is not scientific to delete observations from your dataset, unless they are blatant examples of mistaken data entry;
2) normality is a (weak) requirement for residuals, not for variables;
3) results are what they are. Let's try to understand why the coefficients do not reach statistical significance, instead;
4) it depends on the usual metric used in your research field for that variable.

[Shamsudini Amidu]

Carlo,
Thank you very much for taking the time to review my questions and for your clear and thoughtful explanations. I truly appreciate the clarity you provided

Your guidance has helped me better understand the issues and how to approach them more rigorously. I am grateful for your support and willingness to explain these points in such a clear manner.

Many thanks once again to you, Joseph and all of you for your time and insight.



Kind regards,

Shamsudini Amidu