6
$\begingroup$

I am trying to do break stick linear regression to do two things,

1) calculate the rate of phosphorus needed to ahieve maximum yeild (critical value (CV)) for 4 different cultivars of clover

2) determine statistically if the CV differ between cultivars.

Here is a published example of what I am try to do....

example graph

and here are the methods of that example....

"Shoot dry matter at 250 mg P kg−1 was assumed to represent the maximum growth of the cultivars based on previous experience growing T. subterraneum in this soil (Haling et al. 2016a). The intersection of maximum shoot yield and the linear response to P-application rates between 10 and 80 mg P kg−1 soil was considered to represent the critical external requirement for P. This was analysed by means of a broken stick linear regression with the ‘right side’ regression forced to be a horizontal line (Eq. 2)

equation 2

Where y is the shoot dry mass, x is the P application, k is the maximum shoot dry mass, c is shoot dry mass at 0 mg P kg−1 and m is the gradient. The Breakpoint of the broken stick model (i.e. the critical P requirement) has a distribution which is approximately an F distribution on 1 and n-3 degrees of freedom. Utilising this, error estimates and confidence intervals for the Breakpoint were formed and generated in R2LINES. To test whether the breakpoints were pairwise significantly different, a difference of unpaired means t test was used."

Here is my attempt at this. I have made a rough and ready graph that calculates the critical value, i.e phosphorus applied to achieve maximum yield.

Critical value of one culitvar

So this graph is just of one cultivar. There is a regression line for the first 4 points, a vertical line running from maximum growth, and the intersection of these two lines is the 'critical value'. However, this is a little crude and does not achieve what I ultimately need.

Here is my data...

structure(list(cultivar = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Dinninup", 
"Riverina", "Seaton Park", "Yarloop"), class = "factor"), P = c(12.1, 
12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 18.24, 
18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 48.35, 
12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 
18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 
48.35, 12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 
18.24, 18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 
48.35, 48.35, 12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 
18.24, 18.24, 18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 
48.35, 48.35, 48.35), shoot = c(1.24, 1.12, 1.28, 1.28, 1.37, 
1.4, 1.39, 1.34, 1.34, 1.53, 1.25, 1.4, 1.44, 1.83, 1.65, 1.71, 
1.52, 1.75, 1.63, 1.7, 1.23, 1.22, 1.26, 0.89, 1.2, 1.55, 1.4, 
1.19, 1.75, 1.92, 1.63, 1.64, 1.34, 1.54, 1.66, 1.88, 1.9, 2.18, 
2.03, 1.68, 0.9, 1.49, 1.41, 1.57, 0.94, 1.83, 1.6, NA, 1.98, 
2.04, 1.64, 1.71, 1.97, 1.97, 1.87, 2.21, 2.1, 2.25, 2.1, 2.24, 
1.23, 1.32, 1.47, 1.54, 1.38, 1.09, 1.41, NA, 1.23, 1.14, 1.63, 
1.61, 1.42, 1.12, 1.74, 1.89, 1.4, 1.58, 1.71, 1.64)), class = "data.frame", row.names = c(5L, 
6L, 7L, 8L, 13L, 14L, 15L, 16L, 21L, 22L, 23L, 24L, 29L, 30L, 
31L, 32L, 37L, 38L, 39L, 40L, 45L, 46L, 47L, 48L, 53L, 54L, 55L, 
56L, 61L, 62L, 63L, 64L, 69L, 70L, 71L, 72L, 77L, 78L, 79L, 80L, 
85L, 86L, 87L, 88L, 93L, 94L, 95L, 96L, 101L, 102L, 103L, 104L, 
109L, 110L, 111L, 112L, 117L, 118L, 119L, 120L, 125L, 126L, 127L, 
128L, 133L, 134L, 135L, 136L, 141L, 142L, 143L, 144L, 149L, 150L, 
151L, 152L, 157L, 158L, 159L, 160L))

Can anyone please show me how I can do the desired analysis to emulate the published example above???

Based on your advice, I was able to graph it in ggplot.....

fit1 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
                data = subset(isosub, cultivar == "Yarloop"),
                start = list(c = 1, m = 0.05, bp = 25), na.action = na.omit)
    mm1 <- data.frame(P = seq(0, max(yar$P), length.out = 100))
    mm1$shoot <- predict(fit1, newdata = mm1)

    fit2 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
                data = subset(isosub, cultivar == "Dinninup"),
                start = list(c = 1, m = 0.05, bp = 25), na.action = na.omit)
    mm2 <- data.frame(P = seq(0, max(din$P), length.out = 100))
    mm2$shoot <- predict(fit2, newdata = mm2)

    fit3 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
                data = subset(isosub, cultivar == "Riverina"),
                start = list(c = 1, m = 0.05, bp = 25), na.action = na.omit)
    mm3 <- data.frame(P = seq(0, max(riv$P), length.out = 100))
    mm3$shoot <- predict(fit3, newdata = mm3)

    fit4 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
                data = subset(isosub, cultivar == "Seaton Park"),
                start = list(c = 1, m = 0.05, bp = 25), na.action = na.omit)
    mm4 <- data.frame(P = seq(0, max(seat$P), length.out = 100))
    mm4$shoot <- predict(fit4, newdata = mm4)

    tgq <- summarySE(isosub, measurevar="shoot", 
    groupvars=c("P","cultivar"),na.rm = TRUE)
    ggplot(tgq, aes(x = P, y = shoot)) +
      geom_point(aes(shape=cultivar,colour=cultivar),size=4)+
      scale_color_manual(values = c("#009E73", "#F0E442", "#0072B2", 
      "#D55E00"))+
      theme_bw()+
      geom_line(data = mm1, aes(x = P, y = shoot), colour = "#D55E00")+
      geom_line(data = mm2, aes(x = P, y = shoot), colour = "#009E73")+
      geom_line(data = mm3, aes(x = P, y = shoot), colour = "#F0E442")+
      geom_line(data = mm4, aes(x = P, y = shoot), colour = "#0072B2")

ggplot of broken stick model

This code you gave gives us the y value for the break point right?

coef(fit1)[["c"]] + coef(fit1)[["m"]] * coef(fit1)[["bp"]]

How do we do something similar to get the x value for the break point, as that is the "critical value' we neeed?

First we fit a simple nls model to get decent starting values:

fit1 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
        data = isosub,
        start = list(c = 1, m = 0.05, bp = 25), na.action = na.omit)

library(nlme)

Then we fit the same model with the gnls function:

fita <- gnls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
         data = isosub, 
         params = c+ m +bp ~ 1, start = as.list(coef(fit1)), na.action = 
na.omit)

Now, we stratify bp by cultivar:

fitb <- gnls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
         data = isosub, 
         params = list(bp ~ cultivar,c+ m +bp ~ 1), start = c(coef(fit1)[1], 
0, 0, 0, coef(fit1)[2]), 
         na.action = na.omit)

If you are interested in pairwise comparisons, do the last fit using subsets of two cultivars from your data:

  fitab <- gnls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), data = 
isosub[isosub$cultivar %in% c("Dinninup", "Yarloop"),], 
            params = list(bp ~ cultivar), start = c(coef(fit1)[1], 0, 
coef(fit1)[2], 0), 
            na.action = na.omit)
summary(fitab)$tTable

I am trying to create a break stick model with this data....

structure(list(pot = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 41L, 42L, 43L, 
44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 
57L, 58L, 59L, 60L, 81L, 82L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 
91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 121L, 122L, 
123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L, 132L, 133L, 
134L, 135L, 136L, 137L, 138L, 140L), rep = c(1L, 2L, 3L, 4L, 
1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 
1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 
1L, 2L, 3L, 4L, 1L, 2L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 4L), cultivar = structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Seaton 
Park", 
"Dinninup", "Yarloop", "Riverina"), class = "factor"), Waterlogging = 
structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Non- 
waterlogged", 
"Waterlogged"), class = "factor"), P = c(12.1, 12.1, 12.1, 12.1, 
15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 18.24, 18.24, 24.39, 
24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 48.35, 12.1, 12.1, 
12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 18.24, 
18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 48.35, 
12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 18.24, 
18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 48.35, 
12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 
18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35
), form = c(1.65, 0.61, 0.47, 0.57, 0.52, 0.61, 0.48, 0.8, 0.69, 
0.63, 0.39, 0.68, 0.66, 0.51, 0.4, 0.55, 0.45, 0.41, 0.47, 0.54, 
1.7, 1.78, 1.6, 2.34, 1.52, 1.88, 1.67, 1.7, 1.88, 1.59, 1.97, 
1.6, 1.97, 2.13, 1.52, 2.5, 1.88, 1.61, 1.61, 1.65, 0.05, 0.05, 
0.02, 0.05, 0.31, 0, 0.07, 0.12, 0, 0, 0, 0, 0, 0, 0.05, 0, 0, 
0, 0.03, 0.04, 0.08, 0.08, 0.06, 0.05, 0.12, 0.1, 0.13, 0.05, 
0.07, 0.06, 0.09, 0.05, 0.12, 0.05, 0.1, 0.06, 0.05, 0.06), G = c(0.4, 
0.23, 0.19, 0.12, 0.26, 0.25, 0.19, 0.23, 0.25, 0.4, 0.18, 0.26, 
0.39, 0.38, 0.21, 0.22, 0.28, 0.28, 0.25, 0.28, 1.02, 0.67, 0.8, 
0.78, 0.76, 0.66, 0.79, 0.81, 0.94, 0.61, 0.74, 0.64, 0.99, 0.85, 
0.86, 1, 0.86, 0.75, 0.91, 0.66, 0.91, 0.42, 0.43, 1.02, 1.48, 
0.53, 0.89, 0.7, 0.59, 0.61, 0.42, 1.04, 0.75, 0.59, 0.52, 0.84, 
0.43, 0.53, 0.66, 0.35, 0.19, 0.31, 0.21, 0.27, 0.25, 0.31, 0.21, 
0.28, 0.1, 0.29, 0.09, 0.27, 0.2, 0.19, 0.21, 0.24, 0.11, 0), 
 BA = c(1.61, 1.17, 0.94, 0.98, 1.25, 1.27, 1.15, 1.31, 1.23, 
1.42, 0.91, 1.25, 1.43, 1.61, 1.07, 1.32, 1.48, 1.38, 1.25, 
1.48, 0.09, 0.19, 0.2, 0.16, 0.1, 0.19, 0.13, 0.21, 0.14, 
0.16, 0.2, 0.14, 0.2, 0.21, 0.2, 0.21, 0.21, 0.21, 0.16, 
0.17, 0.23, 0.1, 0.21, 0.27, 0.35, 0.1, 0.31, 0.29, 0.32, 
0.14, 0.21, 0.36, 0.38, 0.16, 0.31, 0.32, 0.21, 0.12, 0.33, 
3.49, 2.53, 2.34, 2.5, 3.54, 2.76, 1.56, 3.13, 2.63, 1.48, 
1.58, 2.34, 2.68, 2.96, 1.31, 3.54, 2.18, 1.5, 1.17), total = c(3.66, 
2.02, 1.59, 1.67, 2.03, 2.13, 1.83, 2.34, 2.17, 2.44, 1.49, 
2.19, 2.48, 2.49, 1.69, 2.1, 2.22, 2.07, 1.97, 2.3, 2.81, 
2.64, 2.59, 3.28, 2.38, 2.72, 2.58, 2.73, 2.95, 2.36, 2.91, 
2.38, 3.16, 3.2, 2.58, 3.71, 2.95, 2.57, 2.68, 2.48, 1.19, 
0.57, 0.66, 1.34, 2.14, 0.63, 1.27, 1.11, 0.91, 0.75, 0.63, 
1.41, 1.13, 0.75, 0.89, 1.16, 0.64, 0.64, 1.02, 3.88, 2.79, 
2.73, 2.77, 3.86, 3.13, 1.97, 3.46, 2.95, 1.65, 1.94, 2.53, 
3, 3.28, 1.55, 3.85, 2.48, 1.66, 1.23), F2 = c(1.97, 2.21, 
1.25, 1.53, NA, 1.27, 0.78, 0.66, 1.21, 1.8, 1.36, 1.61, 
0.71, 0.14, 2.01, 1.29, 1.18, 0.97, 0.55, 1.1, 2.76, 2.34, 
2.43, 1.81, 1.7, 1.44, 1.88, 1.65, 2.34, 0.88, 1.95, 1.88, 
2.01, 1.33, 1.88, 2.02, 3.61, 1.44, 2.08, 2.01, 0.18, 0.16, 
0.15, 0.49, 0.1, 0.3, 0.15, 0.3, 0.45, 0.03, 0.07, 0.24, 
0.16, 0.04, 0.09, 0.08, 0.09, 0.26, 0.09, 0.3, 0.1, 0.3, 
0.16, NA, 0.17, 0.35, 0.25, 0.11, 0.1, 0.02, 0.09, 0.09, 
0.2, 0.39, 0.03, 0.09, 0.27, 0.05), G2 = c(0.69, 0.88, 0.31, 
0.54, NA, 0.44, 0.39, 1.25, 0.36, 0.36, 0.26, 0.8, 0.28, 
0.76, 0.76, 0.45, 0.35, 0.42, 0.23, 0.44, 0.55, 0.76, 0.69, 
0.97, 0.68, 0.87, 0.56, 0.99, 0.7, 0.47, 0.72, 0.94, 0.67, 
0.87, 0.63, 0.94, 0.72, 0.72, 0.69, 1.34, 0.58, 0.94, 0.7, 
1.16, 0.94, 0.87, 0.82, 1.14, 1.05, 0.63, 0.97, 0.6, 1.09, 
0.6, 0.59, 0.82, 0.85, 0.68, 0.94, 0.3, 0.31, 0.42, 0.25, 
NA, 0.39, 0.41, 0.5, 0.16, 0.29, 0.25, 0.29, 0.45, 0.35, 
0.39, 0.11, 0.18, 0.38, 0.21), BA2 = c(1.97, 1.76, 1.88, 
2.14, NA, 1.54, 1.72, 1.39, 1.69, 2.45, 1.94, 1.93, 1.14, 
0.56, 2.08, 2.07, 1.67, 1.94, 1.56, 1.32, 0.11, 0.23, 0.14, 
0.06, 0.17, 0.29, 0.14, 0.11, 0.16, 0.12, 0.14, 0.07, 0.13, 
0.29, 0.13, 0.07, 0.07, 0.14, 0.14, 0.2, 0.36, 0.38, 0.29, 
0.54, 0.33, 0.33, 0.35, 0.4, 0.38, 0.35, 0.35, 0.24, 0.39, 
0.3, 0.18, 0.33, 0.43, 0.26, 0.38, 4.23, 2.6, 4.66, 3.75, 
NA, 2.76, 4.1, 4.25, 1.71, 2.79, 2.47, 2.46, 2.68, 1.58, 
3.88, 1.39, 2.23, 4.13, 2.14), total2 = c(4.63, 4.85, 3.44, 
4.21, NA, 3.25, 2.89, 3.3, 3.26, 4.61, 3.56, 4.34, 2.13, 
1.46, 4.85, 3.81, 3.2, 3.33, 2.34, 2.86, 3.42, 3.33, 3.26, 
2.84, 2.55, 2.6, 2.58, 2.75, 3.2, 1.47, 2.81, 2.89, 2.81, 
2.49, 2.64, 3.03, 4.4, 2.3, 2.91, 3.55, 1.12, 1.48, 1.14, 
2.19, 1.37, 1.5, 1.32, 1.84, 1.88, 1.01, 1.39, 1.08, 1.64, 
0.94, 0.86, 1.23, 1.37, 1.2, 1.41, 4.83, 3.01, 5.38, 4.16, 
NA, 3.32, 4.86, 5, 1.98, 3.18, 2.74, 2.84, 3.22, 2.13, 4.66, 
1.53, 2.5, 4.78, 2.4), Shoot.bag = c(3.83, 3.89, 3.98, 3.7, 
3.94, 4.41, 4.81, 4.41, 4.13, 4.26, 4.59, 3.78, 3.95, 4.35, 
4.92, 4.15, 4.37, 4.54, 4.91, 4.44, 3.62, 3.7, 4.37, 4.63, 
4.91, 4.21, 4.94, 4.39, 4.27, 4.66, 4.89, 4.77, 4.77, 4.8, 
5.23, 4.74, 4.66, 4.42, 5.09, 4.82, 4.73, 4.62, 4.81, 4.85, 
4.68, 4.85, 4.83, 5.08, 4.87, 4.9, 5.36, 4.54, 5.35, 4.65, 
5.04, 5.05, 5.2, 5.21, 4.61, 4.25, 4.09, 3.76, 4.04, 3.77, 
3.84, 4.28, 4.66, 3.94, 4.21, 4, 4.66, 3.85, 4.32, 4.47, 
4.26, 4.95, 5.06, 4.75), shoot = c(0.37, 0.43, 0.52, 0.33, 
0.48, 0.95, 1.35, 0.95, 0.67, 0.8, 1.13, 0.32, 0.58, 0.98, 
1.46, 0.69, 1, 1.17, 1.45, 0.98, 0.25, 0.24, 0.91, 1.17, 
1.54, 1.01, 1.48, 0.93, 0.9, 1.29, 1.43, 1.31, 1.31, 1.43, 
1.77, 1.28, 1.29, 1.05, 1.63, 1.36, 1.36, 1.16, 1.35, 1.39, 
1.22, 1.39, 1.37, 1.71, 1.67, 1.44, 1.9, 1.08, 1.89, 1.19, 
1.58, 1.68, 2, 1.75, 1.24, 0.88, 0.72, 0.3, 0.58, 0.4, 0.47, 
0.82, 1.2, 0.57, 0.84, 0.54, 1.29, 0.48, 0.95, 1.01, 0.8, 
1.58, 1.6, 1.38), root.bag = c(2.98, 2.99, 2.91, 2.95, 3.16, 
3.01, 3.01, 3.01, 3, 2.98, 2.97, 2.98, 3.02, 3.03, 3.17, 
3.14, 2.96, 3.15, 2.93, 3.16, 2.84, 2.98, 3.06, 3.08, 3.03, 
3, 3.06, 3.05, 2.99, 3.01, 3.05, 3.05, 3.08, 3.14, 3.13, 
3.06, 3.01, 3.09, 3.08, 3.04, 3.12, 3.11, 3.24, 3.16, 3.18, 
3.16, 3.1, 3.22, 3.1, 3.08, 3.29, 3, 3.17, 3.04, 3.11, 3.21, 
3.14, 3.04, 3.23, 3.03, 2.97, 2.94, 3, 3, 3.04, 3.04, 3.02, 
3, NA, 3.02, 3.14, 2.98, 3.05, 3.01, 2.88, 2.95, 3.03, 3.04
), root = c(0.11, 0.12, 0.04, 0.08, 0.29, 0.14, 0.14, 0.14, 
0.13, 0.11, 0.1, 0.11, 0.15, 0.16, 0.3, 0.27, 0.09, 0.28, 
0.02, 0.29, 0.02, 0.11, 0.19, 0.21, 0.16, 0.13, 0.19, 0.18, 
0.12, 0.14, 0.18, 0.18, 0.21, 0.27, 0.26, 0.19, 0.14, 0.22, 
0.21, 0.17, 0.25, 0.24, 0.37, 0.29, 0.31, 0.29, 0.23, 0.35, 
0.23, 0.21, 0.42, 0.13, 0.3, 0.17, 0.24, 0.34, 0.27, 0.17, 
0.36, 0.16, 0.1, 0.07, 0.13, 0.13, 0.17, 0.17, 0.15, 0.13, 
0.19, 0.15, 0.27, 0.11, 0.18, 0.14, 0.18, 0.08, 0.16, 0.17
), S.R = c(0.229166667, 0.218181818, 0.071428571, 0.195121951, 
0.376623377, 0.128440367, 0.093959732, 0.128440367, 0.1625, 
0.120879121, 0.081300813, 0.255813953, 0.205479452, 0.140350877, 
0.170454545, 0.28125, 0.082568807, 0.193103448, 0.013605442, 
0.228346457, 0.074074074, 0.314285714, 0.172727273, 0.152173913, 
0.094117647, 0.114035088, 0.113772455, 0.162162162, 0.117647059, 
0.097902098, 0.111801242, 0.120805369, 0.138157895, 0.158823529, 
0.128078818, 0.129251701, 0.097902098, 0.173228346, 0.114130435, 
0.111111111, 0.155279503, 0.171428571, 0.215116279, 0.172619048, 
0.202614379, 0.172619048, 0.14375, 0.169902913, 0.121052632, 
0.127272727, 0.181034483, 0.107438017, 0.136986301, 0.125, 
0.131868132, 0.168316832, 0.118942731, 0.088541667, 0.225, 
0.153846154, 0.12195122, 0.189189189, 0.183098592, 0.245283019, 
0.265625, 0.171717172, 0.111111111, 0.185714286, 0.184466019, 
0.217391304, 0.173076923, 0.186440678, 0.159292035, 0.12173913, 
0.183673469, 0.048192771, 0.090909091, 0.109677419), SPAD_17NOV = c(43, 
39.9, 45, 46, 41, 41.3, 43.5, 43.2, 40, 39.6, 42.9, 43.9, 
42.6, 40.3, 38.4, 39.4, 41.6, 38.2, 36.5, 40.4, 42.6, 43.6, 
48, 43.2, 43, 45.3, 45.2, 48.5, 44.2, 46.8, 47.4, 48.7, 47.7, 
47.4, 43.1, 45.7, 43.9, 44.9, 47.9, 43.9, 52, 47.4, 51.2, 
47.4, 44.8, 47.7, 45.2, 44.2, 44.6, 48.1, 41.5, 44.8, 45.3, 
43.3, 46.6, 44.8, 42.1, 40.6, 46.8, 42.5, 46.7, 44.5, 45.3, 
43.9, 42.2, 43.5, 45.9, 41.1, 44.6, 46.7, 45.8, 42.8, 39, 
43.6, 43.4, 38.5, 39.5, 38.2), plant.height = c(60L, 80L, 
90L, 70L, 130L, 120L, 100L, 120L, 140L, 100L, 110L, 110L, 
130L, 160L, 140L, 130L, 160L, 150L, 170L, 190L, 30L, 140L, 
80L, 70L, 150L, 110L, 110L, 90L, 128L, 120L, 110L, 140L, 
120L, 150L, 130L, 120L, 180L, 160L, 150L, 160L, 80L, 110L, 
70L, 120L, 60L, 90L, 90L, 130L, 150L, 90L, 165L, 140L, 140L, 
150L, 130L, 170L, 210L, 200L, 160L, 50L, 60L, 40L, 40L, 110L, 
90L, 70L, 90L, 80L, 100L, 100L, 120L, 130L, 120L, 120L, 110L, 
140L, 160L, 150L), leaf.discolour.19NOV = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 
6L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
3L, 1L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 1L, 3L, 
3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 5L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
5L, 5L), .Label = c("", " 1 D", "1", "1  D", "2", "D"), class = "factor"), 
deformation.26NOV = structure(c(2L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 4L, 1L, 
1L, 1L, 1L, 2L, 1L, 1L, 4L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 4L, 1L, 1L, 1L, 
1L, 2L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("", 
"D", "D (bad) 1", "D 1"), class = "factor"), herb.dmg.30NOV = structure(c(2L, 
1L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 
2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
3L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 
3L, 1L), .Label = c("", "1", "D"), class = "factor"), herb.dmg.11.DEC = c(3L, 
2L, 4L, 4L, 4L, 4L, 3L, 4L, 4L, 4L, 3L, 4L, 4L, 4L, 3L, 4L, 
4L, 4L, 4L, 4L, 2L, 3L, 2L, 0L, 0L, 2L, 1L, 3L, 3L, 2L, 0L, 
2L, 0L, 2L, 1L, 2L, 3L, 4L, 2L, 2L, 0L, 0L, 1L, 1L, 1L, 0L, 
1L, 2L, 1L, 0L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 3L, 2L, 1L, 1L, 
3L, 2L, 4L, 4L, 2L, 3L, 3L, 4L, 4L, 3L, 4L, 3L, 3L, 4L, 3L, 
3L, 4L), X.plant.pot = structure(c(4L, 3L, 5L, 2L, 1L, 4L, 
4L, 3L, 4L, 3L, 4L, 3L, 2L, 4L, 5L, 2L, 4L, 4L, 4L, 2L, 2L, 
3L, 4L, 4L, 5L, 4L, 5L, 4L, 4L, 3L, 4L, 6L, 4L, 5L, 5L, 4L, 
3L, 4L, 4L, 4L, 4L, 8L, 4L, 4L, 5L, 5L, 5L, 5L, 4L, 4L, 8L, 
4L, 4L, 5L, 4L, 4L, 4L, 5L, 4L, 4L, 4L, 3L, 3L, 9L, 2L, 6L, 
4L, 2L, 3L, 9L, 4L, 2L, 3L, 2L, 2L, 3L, 4L, 5L), .Label = c("2", 
"3", "4", "5", "6", "7", "8", "9", "dead"), class = "factor"), 
nod = c(2, 3, 3, 1, 2, 3, 2, 0.5, 2, 2, 2, 2, 2, 2, 2, 1, 
3, 2, 1, 3, 1, 2, 3, 2, 2, 2, 3, 2, 2, 2, 3, 1, 2, 2, 3, 
3, 2, 2, 3, 2, 4, 2, 1, 3, 3, 3, 2, 3, 3, 3, 2, 2, 2, 2, 
2, 3, 3, 2, 2, 4, 2, 2, 3, 0, 0, 3, 3, 0.5, 2, 0, 2, 2, 2, 
0.5, 2, 0.5, 2, 2), root.dis = c(2L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 3L, 2L, 3L, 2L, 2L, 3L, 2L, 
2L, 2L, 2L, NA, 2L, 1L, 2L, 4L, 2L, 1L, 2L, 2L), surface.root = c(2L, 
2L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 
1L, 1L, 1L, 0L, 0L, 0L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 1L, 1L, 2L, 0L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 
1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 
1L, 1L, 0L, 0L, 2L, 2L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 
2L, 1L), X = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA), big.bag = c(3.37, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), 
med.bad = c(3.2, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA), small.bag = c(3.46, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA
)), row.names = c(1L, 2L, 3L, 4L, 9L, 10L, 11L, 12L, 17L, 
18L, 19L, 20L, 25L, 26L, 27L, 28L, 33L, 34L, 35L, 36L, 41L, 42L, 
43L, 44L, 49L, 50L, 51L, 52L, 57L, 58L, 59L, 60L, 65L, 66L, 67L, 
68L, 73L, 74L, 75L, 76L, 81L, 82L, 84L, 89L, 90L, 91L, 92L, 97L, 
98L, 99L, 100L, 105L, 106L, 107L, 108L, 113L, 114L, 115L, 116L, 
121L, 122L, 123L, 124L, 129L, 130L, 131L, 132L, 137L, 138L, 139L, 
140L, 145L, 146L, 147L, 148L, 153L, 154L, 156L), class = "data.frame")

and I am using this model to create the break point...

fit4 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
        data = subset(isosub, cultivar == "Seaton Park"),
        start = list(c = 1.5, m = 0.1, bp = 70), na.action = na.omit)
mm4 <- data.frame(P = seq(0, max(seat$P), length.out = 100))
mm4$shoot <- predict(fit4, newdata = mm4)

However I am having trouble with the estimates. I am using 70 for the bp as I have already calculated it but cant get the c and m to fit. Can you please offer any help?

Improve this question
$\endgroup$
14
  • 1
    $\begingroup$ The published plot does not appear to be an analysis of data: the points instead appear to be specific locations along a set of given graphs. Without a doubt the actual data did not all fall on those graphs. $\endgroup$
    – whuber ♦
    Commented Jan 13, 2020 at 13:35
  • $\begingroup$ Is that because all the points fall exactly on the line? $\endgroup$
    – Eliott Reed
    Commented Jan 13, 2020 at 23:52
  • $\begingroup$ The output lists bp, presumably for "breakpoint". I might call this "critical x value". The code you list is for the y value at the breakpoint. I might call this the "plateau value". $\endgroup$
    – Sal Mangiafico
    Commented Jan 14, 2020 at 3:35
  • $\begingroup$ Pairwise comparisons for what? $\endgroup$
    – Sal Mangiafico
    Commented Jan 14, 2020 at 3:36
  • 1
    $\begingroup$ As to @whuber 's about the artificial nature of these models, sometimes a quadratic-plateau model is used, which is similar in spirit, but makes a little more sense as a biological response. One practical consideration here is that the linear-plateau and quadratic-plateau models will give different critical x values. If we have enough data, sometimes we use a Cate-Nelson approach and set the plateau value to something that we are actually interested in, say 90% of maximum crop yield. $\endgroup$
    – Sal Mangiafico
    Commented Jan 15, 2020 at 11:59

2 Answers 2

Reset to default
6
$\begingroup$

Here is how to do this for one cultivar:

plot(shoot ~ P, data = subset(DF, cultivar == "Dinninup"))

resulting plot illustrating relationship between "shoot" and "P"

fit1 <- nls(shoot ~ ifelse(P < bp, m * P + c, m * bp + c), 
            data = subset(DF, cultivar == "Dinninup"),
            start = list(c = 1, m = 0.05, bp = 25), na.action = na.omit)
summary(fit1)
#Formula: shoot ~ ifelse(P < bp, m * P + c, m * bp + c)
#
#Parameters:
#    Estimate Std. Error t value Pr(>|t|)    
#c   0.831689   0.105243   7.903 4.31e-07 ***
#m   0.033129   0.005829   5.684 2.69e-05 ***
#bp 24.700463   2.193671  11.260 2.65e-09 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 0.106 on 17 degrees of freedom
#
#Number of iterations to convergence: 2 
#Achieved convergence tolerance: 2.205e-08

curve(predict(fit1, newdata = data.frame(P = x)), add = TRUE)

plot illustrating the fit

#calculate k:
coef(fit1)[["c"]] + coef(fit1)[["m"]] * coef(fit1)[["bp"]]
#[1] 1.65

You can then create a combined model using the approach from my answer to your question at Stack Overflow: https://stackoverflow.com/a/59677502/1412059

The model might be sensitive to starting values (in particular for the break point). You should take care and fit the model repeatedly with slightly different starting values.

Improve this answer
$\endgroup$
13
  • $\begingroup$ Btw., I assume you have a fully randomized experimental design. If that's not the case (for example if you have a block design), you'd probably need a nonlinear mixed-effects model. $\endgroup$
    – Roland
    Commented Jan 13, 2020 at 8:02
  • $\begingroup$ yes it has four reps, each one a block. So you are saying that it needs to be included as a random effect? $\endgroup$
    – Eliott Reed
    Commented Jan 13, 2020 at 9:07
  • $\begingroup$ You need to consider this, yes. $\endgroup$
    – Roland
    Commented Jan 13, 2020 at 12:27
  • $\begingroup$ Thank you Roland. I have updated the post with your advice, though I don't understand how to implement the ANOVA and pairwise comparisons from where I am up to. $\endgroup$
    – Eliott Reed
    Commented Jan 14, 2020 at 1:04
  • $\begingroup$ Study my Stack Overflow answer. The approach is exactly the same. $\endgroup$
    – Roland
    Commented Jan 14, 2020 at 6:43
6
$\begingroup$

The package mcp was made just for scenarios like this. See below how I structured your data as df later.

Fit a change point model

First, let's define a slope followed by a joined plateau. We add varying (random) change point locations (the left-hand side of the equation):

model = list(
  shoot ~ 1 + P,  # intercept and slope
  1 + (1|cultivar) ~ 0  # joined plateau
)

Now we fit the model with default priors:

library(mcp)
fit = mcp(model, data = df, iter = 5000)

Inspect the fit

Let's inspect the full fit for each cultivar:

plot(fit, facet_by = "cultivar", cp_dens = FALSE)

enter image description here

You can see raw parameter estimates using summary(fit) and the corresponding plot_pars(fit) (population-level). To focus on the varying change points (i.e., how each cultivar group deviates from the population-level change point (cp_1)), do ranef(fit) and plot_pars(fit, "varying").

enter image description here

Testing

Here are two ideas how to test hypotheses. If you want to test whether a given change point occurs later than another, do this to obtain Bayes Factors:

hypothesis(fit, "`cp_1_cultivar[Dinninup]` < `cp_1_cultivar[Yarloop]`")

I get a BF of around 16 for this one. If you want to test whether a varying change point improves predictive performance in general, fit a null model and use cross-validation:

model_null = list(shoot ~ 1 + P, ~ 0)
fit_null = fit = mcp(model, data = df, iter = 5000)

# Leave-one-out cross-validation
fit$loo = loo(fit)
fit_null$loo = loo(fit_null)
loo::loo_compare(fit$loo, fit_null$loo)

You can read about mcp on the mcp website and the underlying models in the associated preprint.

Perhaps it would be appropriate to use an informative prior that the first slope is positive.

Data

I removed cases with NA values, and put the rest in a data.frame:

df = data.frame(
  P = c(12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 18.24, 
        18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 48.35, 
        12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 18.24, 
        18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 48.35, 
        48.35, 12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 18.24, 
        18.24, 18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 48.35, 
        48.35, 48.35, 12.1, 12.1, 12.1, 12.1, 15.17, 15.17, 15.17, 15.17, 
        18.24, 18.24, 18.24, 18.24, 24.39, 24.39, 24.39, 24.39, 48.35, 
        48.35, 48.35, 48.35),
  shoot = c(1.24, 1.12, 1.28, 1.28, 1.37,
            1.4, 1.39, 1.34, 1.34, 1.53, 1.25, 1.4, 1.44, 1.83, 1.65, 1.71, 
            1.52, 1.75, 1.63, 1.7, 1.23, 1.22, 1.26, 0.89, 1.2, 1.55, 1.4, 
            1.19, 1.75, 1.92, 1.63, 1.64, 1.34, 1.54, 1.66, 1.88, 1.9, 2.18, 
            2.03, 1.68, 0.9, 1.49, 1.41, 1.57, 0.94, 1.83, 1.6, NA, 1.98, 
            2.04, 1.64, 1.71, 1.97, 1.97, 1.87, 2.21, 2.1, 2.25, 2.1, 2.24, 
            1.23, 1.32, 1.47, 1.54, 1.38, 1.09, 1.41, NA, 1.23, 1.14, 1.63, 
            1.61, 1.42, 1.12, 1.74, 1.89, 1.4, 1.58, 1.71, 1.64),
  cultivar = rep(c("Dinninup", "Yarloop", "Riverina", "Seaton"), each = 20)
)

df = df[complete.cases(df),]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.