# Plant Competition Experiments

Competition experiments are a staple of weed science. Typically, we often want to assess the effect of weed density or duration of competition on crop yield. There is an ongoing debate about the appropriateness of using density and not for example plant cover. We will not go into this debate, but stick to density of plant, because the methods of analyzing data remain the same whether the independent variable, x, is density or plant cover.

One of the important questions is,how do we assess competition and when does it start. In principle competition starts at germination and is a question of the resources:

• Light
• Nutrients
• Water
• “Space”

The three first factors are rather easy to quantify, but the fourth one is more intangible it is dependent on the growth habits of the weeds, e.g. prostrate, erect and the ability of the crop to outgrow the weeds. Whatever the reason for competition, it often boils down to the relationship in Figure1; when will the relationship diverts from a straight line.

## Yield, Density and Yield Loss

The general plant competition and crop yield loss relationships are consider the same, a rectangular hyperbola. In Figure 1 we have a classical intra-specific competition relationship,A, and a yield loss relationship in B.

In Figure 1A competition materializes when the curve diverts from the straight line. The initial straight line means that putting a new plant into the system just increases the yield the same way as all the other individuals contribute initially. When individual plants begin compete with each other for resources, because of high density, then the curves diverts from the straight line.

In Figure 1B the percentage yield loss is based upon the yield without the presence of weeds. The competition (inter-specific competition) for resources materializes itself immediately. When the line diverts from the straight line relationship there will also be some intra-specific competition among the weeds. If there is no competition between crop and weeds at all, then the slope of the curve in Figure 1B would be zero, viz no change in yield whatever the density of weeds.

Figure 1. A, competition within the same species, often denoted intra-specific competition. B, yield loss function based on the percentage yield loss relative to the yield in weed free environment.

If the yield is a crop and the density is weeds per unit area then the the competition (inter-specific competition) materializes in exactly the same way. But now the competition begins from the very start. If there is no competition between crop and weed then the slope of the curve would be zero, viz no change in yield whatever the density of weeds.

The first example will utilize a study conducted near Lingle, Wyoming over two years. In this study, volunteer corn densities ranging from 0 to 2.4 plants/$m^2$ were planted along with dry edible beans to document the bean yield loss from increasing volunteer corn density.

read.csv("http://rstats4ag.org/data/volcorn.csv")->VolCorn
VolCorn$yr<-as.factor(VolCorn$yr)
VolCorn$reps<-as.factor(VolCorn$reps)


##     yr reps dens y.pct y.kg
## 1 2009    1    0 3.284 6273
## 2 2009    2    0 9.467 5872
## 3 2009    3    0 4.625 6186


The variable 'yr' is the year the study was completed (either 2008 or 2009), reps denotes the replicate (1 through 4), 'dens' is the volunteer corn density in plants/$m^2$ (0 to 2.4), 'y.pct' is the percentage dry bean yield loss as compared with the zero volunteer corn density, and 'y.kg' is the dry bean yield in kg/ha.

Cousens 1985 proposed a re-parameterization of a rectangular hyperbola (perhaps better known as Michaelis-Menten) model as a tool to analyze competition experiments, and the drc is well suited for this type of analysis Ritz and Streibig 2005. This package must be loaded with the code:

library(drc)


The drm() function can be used to fit a variety of non-linear models, including the Michaelis-Menten model. The code to fit the Michaelis-Menten model to the volunteer corn data is for one of the two years, 2009, by using the argument data=VolCorn, subset=yr==2009,

drm(y.pct ~ dens, fct=MM.2(names=c("Vmax","K")), data=VolCorn,
subset=yr==2009, na.action=na.omit)->y2009.MM2
summary(y2009.MM2)

##
## Model fitted: Michaelis-Menten (2 parms)
##
## Parameter estimates:
##
##                  Estimate Std. Error t-value p-value
## Vmax:(Intercept)   66.741     73.875   0.903    0.38
## K:(Intercept)       2.649      4.821   0.549    0.59
##
## Residual standard error:
##
##  13.32 (18 degrees of freedom)


Obviously, the Vmax and K parameter of the Michaelis-Menten model were non-significant, the reason is that the range of density of weeds were not large enough, we only catch the linear part ( Figure 2.

plot(y2009.MM2,log="", ylab= "Percent Yield Loss",xlab="Density of Weeds",ylim=c(0,70),
xlim=c(0,10))
abline(h=66.7, lty=2)


Figure 2. Yield loss curve with a two parameter Michaelis-Menten's curve (the argument in drm() is fct=MM.2(). The Vmax is the upper limit, shown with the broken horizontal line, and k is the rate constant.

All functions in drc package are defined in the getMeanFunctions() by writing ?MM.3 or ?MM.2 you can see the help on the function below.

The suffix 3 or 2 defines how many asymptotes we use. The most common one is MM.2 where there only are one upper limit Vmax. If you compare the model above with the log-logistic models, used in the selectivity and dose-response chapter, they look almost identical except that the b in the log-logistic is not existing in MM.2 or MM.3 or to paraphrase b=1. The standard errors reveals that none of the two parameters are significantly different from zero. This is due to the fact that only the first part of the curve is supported by experimental data as seen in the graph in Figure 2 above, no data to support the upper limit of the curve.

The tradition in weed science, as mentioned above, is to reparametrizise the Michaelis-Menten model and use:

where A now is the upper limit and I is the initial slope of the curve as shown in the first graph in this chapter. However, in the drc package we have this reparametrization as shown below:

yieldLoss(y2009.MM2, "as")

##
## Estimated A parameters
## (Asymptotically-based confidence interval(s))
##
##   Estimate Std. Error Lower Upper
## 1     66.7       73.9 -88.5   222
##
##
## Estimated I parameters
## (Asymptotically-based confidence interval(s))
##
##   Estimate Std. Error Lower Upper
## 1     25.2       18.7 -14.2  64.6


The upper limit, which is called Vmax in the Michaelis-Menten and A in the yieldLoss function is the same 67% and the rate constant in the Michaelis-Menten is 2.64 (corresponding to ED50 in the Log-logistic), but for the rectangular hyperbola the initial slope is 25. Of course the parameters of the YieldLoss function were not different from zero either.

Looking at the yield loss curve one could suspect that a straight line relationship applies. We test the regression against the most general model an ANOVA:

lm.Y2009<-lm(y.pct ~ dens,data=VolCorn,     subset=yr==2009)
ANOVA.Y2009<-lm(y.pct ~ factor(dens),data=VolCorn,     subset=yr==2009)
anova(lm.Y2009,ANOVA.Y2009)

## Analysis of Variance Table
##
## Model 1: y.pct ~ dens
## Model 2: y.pct ~ factor(dens)
##   Res.Df  RSS Df Sum of Sq    F Pr(>F)
## 1     18 3242
## 2     15 2637  3       605 1.15   0.36


The test for lack of fit is non-significant (p=0.36) so we can confidently assume the straight line gives a good description of the relationship.

summary(lm.Y2009)

##
## Call:
## lm(formula = y.pct ~ dens, data = VolCorn, subset = yr == 2009)
##
## Residuals:
##    Min     1Q Median     3Q    Max
## -30.01  -6.50   1.46   7.37  27.64
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)     2.50       4.28    0.58   0.5672
## dens           13.00       3.48    3.74   0.0015 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.4 on 18 degrees of freedom
## Multiple R-squared:  0.437,  Adjusted R-squared:  0.406
## F-statistic:   14 on 1 and 18 DF,  p-value: 0.00149


The yield loss percent is increasing 13% with each volunteer corn plant being added into the system

library(dplyr)
Yield.m<-subset(VolCorn,yr==2009) %.%
group_by(dens) %.%
summarise(Yieldloss=mean(y.pct))

plot(Yieldloss ~ dens,data=Yield.m, ylab="Percent Yield Loss")
abline(lm.Y2009)


Figure 3. The straight line fit with a slope of 13 (3.4) meaning that for every unit density of corn we get a yield decrease of 13%. There is no indications that the density of the corn is so high that the intraspecific competition among corn plants has anything to say. The broken line is the nonlinear fit from the previous graph.

The assumption of a straight line relationship in Figure 3 is justified by the test for lack of fit and we can conclude we loose 13% yield per each volunteer corn plant.

However, yield loss in percent of yield with no weed competition is difficult to find, and in practice so dependent on the distribution of actual yields in response to weed density. Therefore, one can get an estimated weed-free yield by using the following re-parameterization of the rectangular hyperbolic

where Y0 is the intercept with the yield axis when weed density is zero. Of course it requires that the curve is well described over the range of relevant weed densities Ritz, Kniss and Streibig 2015.

## Replacement Series

The replacement series can assess interference, niche differentiation, resource utilization, and productivity in simple mixtures of two species. Of course more than two species can be used as long as the total density remains the same, but the interpretation of results becomes very difficult. The species are growing at the same total density, but the proportion between the two species vary. One of the good things about replacement series is that if the replacement graphs looks like the one in Figure 5, it could be the reference, because with linear relationships in Figure 5 shows no competition; the two species do not interfere with each others growth. As discuss earlier, when there is a straight line relationship between yield and density of a species ( Figure 1), the second species does not interfere. If there is a curved relationship there is intraspecific and/or inter specific competition. The intraspecific competition can only be assessed if a species is grown in pure stand.

The density dependence, maximum density determined by experimenter, impedes generalization for a replacement series. And as always, we only get a snapshot of what is going on in an otherwise dynamic competition scenario. In order to avoid criticisms, however, researchers should appreciate the assumptions and limitations of this methodology Jollife 2000

Figure 5. The frame of reference for replacement series with two species where there are no competition between species and the Yield Total (YT) does not change which ever the combination of the two species is. The two species do not need to have the same maximum yield in monoculture.

The philosophy of the replacement series is that the carrying capacity, in terms of say biomass, on a unit of land is constant whatever the proportion of the species.

The example uses a barley crop grown together with the weed Amsinckia menziesii. The experiment was run in greenhouse with the intention of having 20 plants in total in pots of 20 cm in diameter. After 20 days the plants were harvested and the actual number of plants were counted and the biomass per species measured.

The dataset Replacement series.csv is a mixture of csv and csv2 files, because the students who did the experiments came form continental Europe or Australia. It means the continental Europeans use the semicolon as variable separator mixed with the Australian's decimal separator of dot. This mess is taken care of by using the read.table() function.

As we have to use the percent of either species as independent variable to fit regression models we have to define new variables Pct.Amsinckia and Pct.Barley.

Replace.1<-read.table("http://rstats4ag.org/data/Replacement series.csv",header=TRUE,sep=";",dec=".")
Replace.1$Pct.Amsinckia<-with(Replace.1,D.Amsinckia*100/(D.Amsinckia+D.Barley)) Replace.1$Pct.Barley<-   with(Replace.1,D.Barley*100/(D.Amsinckia+D.Barley))

##   D.Amsinckia D.Barley Harvest.time B.Amsinckia B.Barley Pct.Amsinckia
## 1           0       14   14.05.2014           0     98.1             0
## 2           0       32   14.05.2014           0    112.7             0
## 3           0       21   14.05.2014           0     94.5             0
##   Pct.Barley
## 1        100
## 2        100
## 3        100


First we try straight lines relationships and illustrate the fit and with an analysis of residuals.

op <- par(mfrow = c(2, 2))
lm.Amsinckia<-lm(B.Amsinckia ~Pct.Amsinckia,data=Replace.1)
plot(B.Amsinckia ~Pct.Amsinckia,data=Replace.1, ylab="Yield of Amsinckia",
xlab="Percent of Amsinckia")
abline(lm.Amsinckia)

lm.Barley<-lm(B.Barley ~Pct.Barley,data=Replace.1)
plot(B.Barley ~Pct.Barley,data=Replace.1,ylab="Yield of Barley",
xlab="Percent of Barley",pch=16)
abline(lm.Barley,lty=2)
plot(lm.Amsinckia,which=1)
plot(lm.Barley,which=1)


Figure 6. Straight line relationships do not look very good to capture the variation in the species. The analysis of residuals shows systematic departure form the expected shotgun distribution of residuals,

Obviously, the relationships in Figure 6 for both species look like being a curved relationship. We are not sure of which relationship to use and resort to a second degree polynomial.

where a is the intercept with the y-axis and b and c are parameters for the x and the $x^2$.
The biological meanings of polynomial parameters in general are not often of interest because they can be hard to interpret. To fit the polynomial we use the lm()function because it is essentially a linear model we are fitting by adding a parameter for the $x^2$ by writing I(x^2). However, we cannot use the abline() function.

First we generate predicted values to get smooth fits. It is done with the predict()function predict(Pol.B.Amsinckia,data.frame(Pct.Amsinckia=seq(0,100,by=1)))
where

• Pol.B.Amsinckia the results of the second degree polynomial fit for Amsinckia
• data.frame(Pct.Amsinckia=seq(0,100,by=1)) generates a total of 101 predicted values for Amsinckia.

The same procedure applies to barley.

op <- par(mfrow = c(1, 2))
Pol.B.Amsinckia<-lm(B.Amsinckia ~Pct.Amsinckia+I(Pct.Amsinckia^2),data=Replace.1)
line.Pol.B.Amsinckia<-predict(Pol.B.Amsinckia,data.frame(Pct.Amsinckia=seq(0,100,by=1)))

Pol.B.Barley<-lm(B.Barley ~Pct.Barley+I(Pct.Barley^2),data=Replace.1)
line.Pol.B.Barley<-predict(Pol.B.Barley,data.frame(Pct.Barley=seq(0,100,by=1)))

plot(B.Amsinckia ~Pct.Amsinckia,data=Replace.1,ylim=c(0,100))
lines(line.Pol.B.Amsinckia~seq(0,100,by=1),col="red")
plot(B.Barley ~Pct.Barley,data=Replace.1,pch=16,ylim=c(0,100))
lines(line.Pol.B.Barley~seq(0,100,by=1),col="red",lty=2)


Figure 7. Fitting second degree polynomials to data.

Overall the second degree polynomials describe the variation reasonably well (Figure 7). However, there is a catch to it, the minimum for the Amsinckia second degree polynomial is at low percent of Amsinckia and the maximum for the second degree polynomial for barley is at rather high proportion of barley. If we can live with that we can use the fit to summarize the experiment by calculation the Yield Total (YT) as shown in the graph in b>Figure 1.

This is a good example of the problem with polynomials. A Second degree polynomial is symmetric with either a minimum or a maximum depending of the parameters. A third degree polynomial does not have this symmetric property, it becomes much more complex, and it is definitely not advisable to use with interpolations.

In order to summarize the experiment on the basis of the fits, we can combine the two curves for the two species and calculate the YT.

We want to use the line.Pol.B.Amsinckia and the line.Pol.B.Barley that define the smooth lines and calculate the YT. The maximum of Amsinckia is the 101 element. In the line.Pol.B.Amsinckia[101] the maximum yield is 65 and for line.Pol.B.Barley[101] is is 101. In order to calculate the YT we just have to sum the two but for line.Pol.B.Barley, however, we must reverse the order.

op <- par(mfrow = c(1, 1))
plot(c(line.Pol.B.Amsinckia[101],line.Pol.B.Barley[101])~c(100,0),type="l",ylim=c(0,120), xlab="", ylab="Yield of Barley")
mtext("Pct of Amsinckia --->",side=1,line=2)
mtext("<---Pct of  Barley",side=1,line=3)
mtext("Yield of Amsinckia",side=4)
Pred.RYT<-line.Pol.B.Amsinckia[1:101]+line.Pol.B.Barley[seq(101,1,by=-1)]
lines(Pred.RYT~seq(0,100,by=1),lty=3,col="red",lwd=2)
lines(line.Pol.B.Barley~seq(100,0,by=-1),col="red",lty=2)
lines(line.Pol.B.Amsinckia~seq(0,100,by=1),col="red",lty=1)

points(B.Amsinckia ~Pct.Amsinckia,data=Replace.1,ylim=c(0,100))
points(B.Barley ~Pct.Amsinckia,data=Replace.1,pch=16,ylim=c(0,100))


Figure 8. Summary of the replacement series experiment with barley and Amsinckia.

It seems in Figure 8 that if Amsinckia gains with its convex curve, and barley suffers with its concave curve, but gain and suffering still leave the YT below the theoretical value of YT.

Another issue is that that we do not test the regressions statistically, but use the fit to illustrate the relationships. In order to apply test statistics it requires more systematic designs with fixed number of plants per unit area, which unfortunately was not the case here. In the first example we had genuine replication with several replicates of the number of volunteer corn per unit area and therefore we could test which model could be used.

It also shows that when fitting curves we can derive predicted values, which can be used to calculate derived parameters such as YT. In replacement analysis content one often scale the results so that the theoretical maximum yield is equal to 1.0 and then calculate the Relative yield Total (RYT), which in mixed cropping research is called Land Equivalent Ratio (LER).