4
$\begingroup$

I have a small data set that I am trying to analyze and have looked at it with both a linear model and a generalized linear model. The data are from seed traps placed beneath shrubs. The "successes" are seeds likely dispersed by birds, and the "fails" are the seeds collected in the seed trap. X represents the total number of fruits on the shrub. Because the total number of seeds on an individual shrub weren't always over the seed trap, X can only be equal to or greater than successes + failures.

The results are different and I haven't come across this before. Here is code to replicate what I am talking about.

library(sjPlot)

success = c(0, 0, 8, 4, 9, 0, 42, 29, 1, 71, 21, 9, 65, 7, 0, 51, 31, 27, 
            29, 39, 32, 9, 8, 37, 15, 53, 102, 14, 10, 3, 9, 51, 46, 13, 
            28, 23, 11, 23, 28, 23)
fail = c(29, 12, 149, 0, 0, 20, 10, 22, 7, 2, 1, 1, 34, 4, 57, 1, 2, 
         6, 3, 14, 5, 4, 6, 2, 12, 0, 15, 10, 16, 10, 29, 14, 3, 22, 10, 
         24, 6, 14, 8, 12)
X = c(29, 12, 157, 50, 119, 26, 65, 104, 45, 73, 22, 77, 123, 50, 
      125, 149, 149, 180, 32, 53, 37, 66, 18, 39, 150, 124, 117, 24, 
      75, 13, 61, 65, 49, 35, 61, 77, 26, 37, 50, 44)
prop = success/(success+fail)
dat = data.frame(prop, success, fail, X)

model.lm = lm(prop~X, data=dat)
plot(dat$X, dat$prop)
plot_model(model.lm,type="pred")

model.glm = glm(cbind(success,fail)~X, family="binomial", data=dat)
plot_model(model.glm, type="pred")

The plot for the linear model shows an increase in the proportion of "successes" as X increases. However, the plot for the generalized model shows the opposite: a decrease in success as X increases. The GLM approach is appropriate, but the when you plot the points (proportion of success as a function of X), it "looks" like success increases as a function of X, even though the errors aren't evenly distributed. Can anyone help me with this?

$\endgroup$
3
  • $\begingroup$ You may want to say more about what these data are & what you want from them. It is a little suspicious that the total is often equal to X, & never exceeds it. $\endgroup$ Jun 3, 2021 at 16:11
  • $\begingroup$ Sure. I will add some additional information to the question. $\endgroup$
    – user44796
    Jun 3, 2021 at 17:55
  • $\begingroup$ Thanks, interesting situation. I don't think the data would suffer from an obvious problem due to this. I guess I don't know for sure, but I think you're OK to move forward w/ a standard logistic regression. $\endgroup$ Jun 3, 2021 at 19:24

1 Answer 1

5
$\begingroup$

The correlation between X and prop isn't really that strong. It also ignores the fact that the total amount of information contributing to the observed proportions differs across your data. Moreover, that total is confounded with X.

total = success + fail

cor(X, prop)   # [1] 0.1866438
cor(X, total)  # [1] 0.4947678
windows(width=7, height=4)
  layout(matrix(1:2, nrow=1))
  plot(X, prop)
  plot(X, total)

scatterplots of X vs prop & total

A linear model that does include weights to indicate the different amounts of information in each prop value does show a negative slope with X:

model.wlm = lm(prop~X, data=dat, weights=total)
summary(model.wlm)
# Call:
# lm(formula = prop ~ X, data = dat, weights = total)
# 
# Weighted Residuals:
#     Min      1Q  Median      3Q     Max 
# -6.2145 -1.0587  0.3519  1.2314  3.0663 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.716082   0.112251   6.379 1.72e-07 ***
# X           -0.001077   0.001135  -0.949    0.349    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 2.097 on 38 degrees of freedom
# Multiple R-squared:  0.02316,   Adjusted R-squared:  -0.002547 
# F-statistic: 0.9009 on 1 and 38 DF,  p-value: 0.3485

In general, I'm not a big fan of the linear probability model, and I wouldn't use it here.


Below, I've made a plot incorporating @whuber's suggestion to make the size of the points larger when the proportion is estimated from a larger total:

dot.cex = (((total - min(total))/max(total))*10) + 1
plot(X, prop, pch=16, col=rgb(0,0,0,alpha=.5), cex=dot.cex)

scatterplot with differently sized dots

$\endgroup$
4
  • $\begingroup$ +1 Modifying the first plot to make the dot areas proportional to total might help clarify your point about the weights. $\endgroup$
    – whuber
    Jun 3, 2021 at 15:33
  • 1
    $\begingroup$ This is what I was looking for. I was planning on using a generalized linear model, but just couldn't figure out why the linear model was incorrect. Great explanation. $\endgroup$
    – user44796
    Jun 3, 2021 at 16:05
  • $\begingroup$ You're welcome, @user44796. $\endgroup$ Jun 3, 2021 at 16:09
  • $\begingroup$ The modified plot is interesting. I usually avoid scaling the dots directly by another quantity, preferring to use the square root so the dot areas reflect the quantity (as in a quick and dirty ... cex=3*sqrt(total/150)...): but the visual exaggeration of not doing so helps drive the point home. $\endgroup$
    – whuber
    Jun 4, 2021 at 22:25

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.