2
$\begingroup$

In my experiment I want to figure out how the size of different planting containers, i.e. their volume, affects the number of regenerated plant shoots from root fragments (terminology here is root suckering). So in other words, how many shoots are regenerated from the initial root material given the container they grow in. After some time, I counted the number of shoots per planting container and scaled the count up to 100 ml to make it easier to compare across container types.

I have the following example dataset:

df <- structure(list(Container = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("A", "B", 
"C", "D"), class = "factor"), Rep = c(1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), StandCount = c(2.24, 
2.88, 0.96, 1.92, 3.2, 2.88, 2.56, 2.24, 7.04, 3.2, 4.8, 3.52, 
8.96, 7.68, 8.96, 6.72, 3, 3.6, 4.2, 3.6, 3, 10.8, 4.8, 4.2, 
9, 10, 11, 32, 21, 24, 20, 19)), .Names = c("Container", "Rep", 
"StandCount"), class = "data.frame", row.names = c(NA, -32L))

The data looks as follows:

#require(ggplot2)
#ggplot(df, aes(x = Container, y = StandCount)) + geom_boxplot()

enter image description here

Now I fitted two models (m1 and m2). The first one is a generalized least squares model using the gls() function including the varIdent() function to account for the observed heteroskedasticity.

The second one is a generalized linear model using the Gamma distribution since Poisson requires integer data.

require(nlme)
m1 <- gls(StandCount ~ Container, weights = varIdent(form= ~1|Container), data = df)
#plot(m1)
m2 <- glm(StandCount ~ Container, family = Gamma, data = df)
#plot(m2)

Here are the estimated marginal means and standard errors. As it can be seen those SEs differ slightly resulting in different significant differences in the Tukey adjusted multiple mean comparison (cf. A-D for example):

require(emmeans)
emmeans(m1, pairwise~ Container)
$emmeans
 Container emmean        SE df  lower.CL  upper.CL
 A           2.36 0.2488832 28  1.850186  2.869814
 B           6.36 0.8055525 28  4.709901  8.010099
 C          18.25 2.8014653 28 12.511459 23.988541
 D           4.65 0.9053412 28  2.795493  6.504507

Confidence level used: 0.95 

$contrasts
 contrast estimate        SE df t.ratio p.value
 A - B       -4.00 0.8431238 28  -4.744  0.0003
 A - C      -15.89 2.8124990 28  -5.650  <.0001
 A - D       -2.29 0.9389279 28  -2.439  0.0927
 B - C      -11.89 2.9149824 28  -4.079  0.0018
 B - D        1.71 1.2118406 28   1.411  0.5032
 C - D       13.60 2.9441213 28   4.619  0.0004

P value adjustment: tukey method for comparing a family of 4 estimates 

> emmeans(m2, pairwise~ Container, type="response")
$emmeans
 Container response        SE df asymp.LCL asymp.UCL
 A             2.36 0.3513056 NA  1.826969  3.332189
 B             6.36 0.9467387 NA  4.923527  8.979967
 C            18.25 2.7166633 NA 14.128046 25.767986
 D             4.65 0.6921885 NA  3.599752  6.565532

Confidence level used: 0.95 
Intervals are back-transformed from the inverse scale 

$contrasts
 contrast    estimate         SE df z.ratio p.value
 A - B     0.26649611 0.06727806 NA   3.961  0.0004
 A - C     0.36893429 0.06360075 NA   5.801  <.0001
 A - D     0.20867505 0.07073414 NA   2.950  0.0168
 B - C     0.10243818 0.02478594 NA   4.133  0.0002
 B - D    -0.05782106 0.03965611 NA  -1.458  0.4631
 C - D    -0.16025924 0.03303521 NA  -4.851  <.0001

P value adjustment: tukey method for comparing a family of 4 estimates 

My question now is what is the correct approach to analyze those data? I think the gamma model is most appropriate also given the fact that there won't be any negative values, however there might be zeros though. In this case, which analysis would be most appropriate?

$\endgroup$
  • 1
    $\begingroup$ Arguably, none of the above is the "correct approach" (although one or more might work). Go back to the original counts. You can handle container type as an offset (in the form of amount of original root material). $\endgroup$ – whuber Jan 18 '18 at 21:18
  • $\begingroup$ Thanks @whuber I see! So assuming the root material was a constant (say 10 grams), but the volume of the containers is different, I can add the container volume as an offset? I will have to read more about it. $\endgroup$ – Stefan Jan 18 '18 at 22:07
  • 1
    $\begingroup$ The volume might not be the right offset. Ideally, you want to use values that are proportional to the expected amount of suckering. That would likely be the number of original root fragments in each container. If you don't know that, one proxy might be the surface area of the container. $\endgroup$ – whuber Jan 18 '18 at 22:31
  • $\begingroup$ Thanks for your response @whuber ! That really helps a lot! $\endgroup$ – Stefan Jan 18 '18 at 22:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.