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
read.table()
orreadr::read_table()
on it... could you share it withdput(callerydata)
so it's copy/pasteable? And maybe drop the columns that are allNA
$\endgroup$