1
$\begingroup$

I have been given a dataset to analyze looking at some herbicide treatments for invasive trees in three states with 3 sites in state 1, 2 sites in state 2, and 1 site in state 3. We hope to answer which of the 7 herbicide treatments differed significantly from the untreated control after one year post-treatment. There are 5 reps of each treatment at each site. The response variable is binary(0=alive, 1=dead) and treatment is a predictor with diameter of tree, and site within state as potential covariates. We want to analyze by state. Any suggestions are welcome.

Looking at the data, we expected a strong effect of treatment, with three treatments being significantly different from the control since all trees died, 2 might be significant since they worked intermediately, while 2 treatments worked very poorly and would likely not be significant. I have run the data many different ways, used similar stack questions to modify, coded the glm different ways, reduced the model, etc. and continue to get p values of almost 1 for everything. I originally thought there wasn't enough power, however I'm relatively new to glms and assume I'm doing something wrong, or could have a better approach.

Example paired down model with code and output for state 1 with three sites is below.

> callerydata<-read.csv("GLMKSoneyear.csv")
> callerydata
   Location Site Rep TREE TRT TREATMENT Crown crown.2  X dead Sprout   Date Lat long DBH MULTISTEM X.1
1     SN KS    1   1  102  NA         8     0      NA NA    1      n 7/7/22  NA   NA   7         n  NA
2     SN KS    1   2  109  NA         8     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
3     SN KS    1   5  121  NA         8     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
4     SN KS    1   4  128  NA         8     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
5     SN KS    1   3  142  NA         8     0      NA NA    1      n 7/7/22  NA   NA   5         n  NA
6     SN KS    1   1  101  NA         2     0      NA NA    1      n 7/7/22  NA   NA   8         n  NA
7     SN KS    1   2  111  NA         2     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
8     SN KS    1   5  123  NA         2     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
9     SN KS    1   4  129  NA         2     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
10    SN KS    1   3  131  NA         2     0      NA NA    1      n 7/7/22  NA   NA   5         n  NA
11    SN KS    1   2  107  NA         3     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
12    SN KS    1   1  115  NA         3     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
13    SN KS    1   5  124  NA         3     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
14    SN KS    1   4  127  NA         3     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
15    SN KS    1   3  134  NA         3     0      NA NA    1      n 7/7/22  NA   NA   8         n  NA
16    SN KS    1   1  104  NA         4    20      NA NA    0      n 7/7/22  NA   NA   3         n  NA
17    SN KS    1   2  110  NA         4     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
18    SN KS    1   5  120  NA         4    60      NA NA    0      n 7/7/22  NA   NA   3         n  NA
19    SN KS    1   4  126  NA         4    10      NA NA    0      n 7/7/22  NA   NA   3         n  NA
20    SN KS    1   3  132  NA         4    40      NA NA    0      n 7/7/22  NA   NA   4         n  NA
21    SN KS    1   5  105  NA         5     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
22    SN KS    1   1  106  NA         5     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
23    SN KS    1   2  113  NA         5     0      NA NA    1      n 7/7/22  NA   NA   4         n  NA
24    SN KS    1   4  130  NA         5    30      NA NA    0      n 7/7/22  NA   NA   3         n  NA
25    SN KS    1   3  136  NA         5     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
26    SN KS    1   1  103  NA         6    20      NA NA    0      n 7/7/22  NA   NA   7         n  NA
27    SN KS    1   2  108  NA         6    10      NA NA    0      n 7/7/22  NA   NA   4         n  NA
28    SN KS    1   5  118  NA         6     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
29    SN KS    1   4  125  NA         6     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
30    SN KS    1   3  137  NA         6    10      NA NA    0      n 7/7/22  NA   NA   4         n  NA
31    SN KS    1   2  112  NA         7    90      NA NA    0      n 7/7/22  NA   NA   3         n  NA
32    SN KS    1   1  116  NA         7     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
33    SN KS    1   5  119  NA         7   100      NA NA    0      n 7/7/22  NA   NA   7         n  NA
34    SN KS    1   3  138  NA         7     0      NA NA    1      n 7/7/22  NA   NA   3         n  NA
35    SN KS    1   4  139  NA         7    90      NA NA    0      n 7/7/22  NA   NA   3         n  NA
36    SN KS    1   2  114  NA         1   100      NA NA    0      n 7/7/22  NA   NA   3         n  NA
37    SN KS    1   1  115  NA         1   100      NA NA    0      n 7/7/22  NA   NA   4         n  NA
38    SN KS    1   1  117  NA         1   100      NA NA    0      n 7/7/22  NA   NA   3         n  NA
39    SN KS    1   5  122  NA         1   100      NA NA    0      n 7/7/22  NA   NA   6         n  NA
40    SN KS    1   3  133  NA         1   100      NA NA    0      n 7/7/22  NA   NA   3         n  NA
41    SN KS    1   3  135  NA         1   100      NA NA    0      n 7/7/22  NA   NA   5         n  NA
42    SN KS    1   4  140  NA         1   100      NA NA    0      n 7/7/22  NA   NA   3         n  NA
43    SN KS    1   4  141  NA         1   100      NA NA    0      n 7/7/22  NA   NA   5         n  NA
44    SG KS    2   1  201  NA         8     0      NA NA    1      n 7/5/22  NA   NA   6         n  NA
45    SG KS    2   2  213  NA         8     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
46    SG KS    2   3  216  NA         8     0      NA NA    1      n 7/5/22  NA   NA   4         n  NA
47    SG KS    2   5  226  NA         8     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
48    SG KS    2   4  235  NA         8     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
49    SG KS    2   1  204  NA         2     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
50    SG KS    2   2  209  NA         2     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
51    SG KS    2   3  220  NA         2     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
52    SG KS    2   5  229  NA         2     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
53    SG KS    2   4  233  NA         2     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
54    SG KS    2   1  202  NA         3     0      NA NA    1      n 7/5/22  NA   NA   5         n  NA
55    SG KS    2   2  212  NA         3     0      NA NA    1      n 7/5/22  NA   NA   3         n  NA
  
 [ reached 'max' / getOption("max.print") -- omitted 72 rows ]
> str(callerydata)
'data.frame':   127 obs. of  18 variables:
 $ Location    : chr  "SN KS" "SN KS" "SN KS" "SN KS" ...
 $ Site        : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Rep         : int  1 2 5 4 3 1 2 5 4 3 ...
 $ TREE        : int  102 109 121 128 142 101 111 123 129 131 ...
 $ TRT         : logi  NA NA NA NA NA NA ...
 $ TREATMENT   : int  8 8 8 8 8 2 2 2 2 2 ...
 $ Crown       : int  0 0 0 0 0 0 0 0 0 0 ...
 $ crown.2     : logi  NA NA NA NA NA NA ...
 $ X           : logi  NA NA NA NA NA NA ...
 $ dead        : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Sprout      : chr  "n" "n" "n" "n" ...
 $ Date        : chr  "7/7/22" "7/7/22" "7/7/22" "7/7/22" ...
 $ Lat         : logi  NA NA NA NA NA NA ...
 $ long        : logi  NA NA NA NA NA NA ...
 $ DBH         : int  7 4 4 3 5 8 4 4 3 5 ...
 $ MULTISTEM   : chr  "n" "n" "n" "n" ...
 $ X.1         : logi  NA NA NA NA NA NA ...
 $ X..Mortality: logi  NA NA NA NA NA NA ...
> site<-as.factor(callerydata$Site)
> site
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2
 [50] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [99] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Levels: 1 2 3
> treatment<-as.factor(callerydata$TREATMENT)
> treatment
  [1] 8 8 8 8 8 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 1 1 1 1 1 1 1 1 8 8 8 8 8 2
 [50] 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 1 1 1 1 1 1 8 8 8 8 8 2 2 2 2 2 3 3 3 3
 [99] 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 1 1 1 1 1 1 1 1
Levels: 1 2 3 4 5 6 7 8
> dbh<-as.numeric(callerydata$DBH)
> dbh
  [1] 7 4 4 3 5 8 4 4 3 5 4 4 4 3 8 3 3 3 3 4 3 3 4 3 3 7 4 3 3 4 3 3 7 3 3 3 4 3 6 3 5 3 5 6 3 4 3 3 3
 [50] 3 3 3 3 5 3 4 3 3 3 4 3 5 5 3 3 4 3 4 5 3 3 3 3 3 3 3 4 3 3 5 3 3 3 3 3 3 4 4 3 4 3 3 3 3 4 4 3 3
 [99] 3 3 3 3 3 3 3 3 3 4 3 3 3 4 4 3 4 3 3 3 3 3 4 3 4 3 4 4 4
> dead<-as.numeric(callerydata$dead)
> dead     
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 1 1 1 0 1 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1
 [50] 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [99] 1 0 0 0 1 0 0 1 1 1 1 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
> m1<-glm(dead~treatment+dbh, family=binomial(link="logit"), data=callerydata)
> coef(m1)
(Intercept)  treatment2  treatment3  treatment4  treatment5  treatment6  treatment7  treatment8 
-17.5640251  42.3426692  42.4068377  18.5110079  21.8141357  21.8365887  19.9715666  42.1111665 
        dbh 
 -0.8633232 
> confint(m1)
Waiting for profiling to be done...
                  2.5 %        97.5 %
(Intercept)          NA  244.72126524
treatment2  -202.891548 1788.62332384
treatment3  -235.238956 1630.39616904
treatment4  -185.246873            NA
treatment5  -162.467758 1102.08737891
treatment6  -154.115330 1147.77089247
treatment7  -154.123272 1156.08550946
treatment8  -215.458754 1816.80108631
dbh           -2.052214   -0.04449088
There were 50 or more warnings (use warnings() to see the first 50)
> summary(m1)

Call:
glm(formula = dead ~ treatment + dbh, family = binomial(link = "logit"), 
    data = callerydata)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.92514  -0.00006   0.00002   0.00018   1.90777  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -17.5640  3708.3908  -0.005   0.9962  
treatment2    42.3427  5482.5951   0.008   0.9938  
treatment3    42.4068  5516.9642   0.008   0.9939  
treatment4    18.5110  3708.3905   0.005   0.9960  
treatment5    21.8141  3708.3905   0.006   0.9953  
treatment6    21.8366  3708.3905   0.006   0.9953  
treatment7    19.9716  3708.3905   0.005   0.9957  
treatment8    42.1112  5652.9478   0.007   0.9941  
dbh           -0.8633     0.4977  -1.735   0.0828 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 171.106  on 126  degrees of freedom
Residual deviance:  60.047  on 118  degrees of freedom
AIC: 78.047

Number of Fisher Scoring iterations: 19
$\endgroup$
1
  • $\begingroup$ I tried importing your data, but just got errors about the wrong number of columns when I ran read.table() or readr::read_table() on it... could you share it with dput(callerydata) so it's copy/pasteable? And maybe drop the columns that are all NA $\endgroup$ Sep 19, 2023 at 14:18

1 Answer 1

0
$\begingroup$

Your conclusions rely on the asymptotic theory for maximum likelihood estimation, which may not be valid. In fact, it is not valid in this case as your standard errors tend to infinity as the maximum likelihood estimate goes away from 0. This is known as the Hauck and Donner phenomenon (see https://doi.org/10.1080/01621459.1977.10479969).

To assess the significance of the treatment variable, use a likelihood ratio test, which does not suffer from the issue as the z test. The function anova in R can do this for you.

$\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.