Hello-

I conducted a field experiment where I added seedlings near and far from 
conspecific trees at multiple plots at multiple sites on multiple 
islands and recorded survival, among other things. Here, I'm trying to 
determine the overall effect of distance (near vs far) on the proportion 
of seeds that survived in each plot (site is a random effect). I'd also 
like to determine whether the effect of distance varies by canopy 
openness. Finally, I would like to test whether the effect of distance 
differs between islands. I'm confused about how to set up the right 
contrasts given the interactions, and how to interpret the output. 
Treatment contrasts use a reference level, which is fine for distance 
where 'near' is a meaningful reference level, but for 'island', the 
reference level is not meaningful (i.e. it doesn't make sense to compare 
each island to one arbitrarily chosen island). A link to my data, and 
output from an example are below. This is for one species, but I have 
data for six species total, and ideally, I'd like to compare these 
effects between species (see question 2 below). I have three specific 
questions at the bottom.

Here's a link to the data: 
https://www.dropbox.com/s/fregbu154zw7whz/aglaifen.csv

 > aglaiafen<-read.csv("aglaifen.csv")
 > aglaiafen$dist<-factor(aglaiafen$dist, levels=c("near", "far"))
 >
 > #response is cbind(number of seedlings that survived/number of 
seedlings that died), family = binomial
 > #Three predictors- island, distance, and canopy openness
 > #island has 4 levels (A, B, C, D). Island A is reference level by 
default, but is not really meaningful as reference level.
 > #dist has 2 levels, near and far, with near set as reference level.
 > #centavgopen is centered at the mean canopy openness.
 > #site is a random effect (3-5 sites per island, 4 islands). There are 
usually 4 near plots and 4 far plots per site.
 >
 > #Here are three approaches to analysis that give qualitatively 
similar answers, but I'm not sure which (if any) is the best approach.
 >
 > #Approach 1: Use the standard treatment or poly contrasts (depending 
on whether unordered or ordered) for all values and without intercept 
suppressed
 > contrasts(aglaiafen$island)<-contr.treatment #just to make sure
 > 
survm1<-glmer(cbind(numalive,numseedplant-numalive)~island*dist+centavgopen*dist+(1|site), 
family=binomial, data=aglaiafen)
 > summary(survm1)
Generalized linear mixed model fit by maximum likelihood ['glmerMod']
Family: binomial ( logit )
Formula: cbind(numalive, numseedplant - numalive) ~ island * dist + 
centavgopen * dist + (1 | site)
Data: aglaiafen

AIC BIC logLik deviance
353.3382 382.2129 -165.6691 331.3382

Random effects:
Groups Name Variance Std.Dev.
site (Intercept) 0.2544 0.5044
Number of obs: 102, groups: site, 13

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.923415 0.316583 6.076 1.24e-09 ***
island2 -0.927197 0.480333 -1.930 0.0536 .
island3 -0.147921 0.494233 -0.299 0.7647
island4 1.922899 0.875871 2.195 0.0281 *
distfar 0.518164 0.337777 1.534 0.1250
centavgopen 0.096598 0.046130 2.094 0.0363 *
island2:distfar -0.467564 0.459629 -1.017 0.3090
island3:distfar -0.954684 0.472885 -2.019 0.0435 *
island4:distfar -0.534748 1.092362 -0.490 0.6245
distfar:centavgopen -0.006505 0.067601 -0.096 0.9233
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 > confint(survm1, method="Wald")
2.5 % 97.5 %
(Intercept) 1.30292450 2.54390628
island2 -1.86863116 0.01423767
island3 -1.11660052 0.82075815
island4 0.20622293 3.63957512
distfar -0.14386687 1.18019517
centavgopen 0.00618545 0.18701081
island2:distfar -1.36841927 0.43329197
island3:distfar -1.88152203 -0.02784606
island4:distfar -2.67573843 1.60624191
distfar:centavgopen -0.13899924 0.12598992

 > #The way I understand it, the intercept represents islandA, 
distancenear and mean openness (0), so the coefficient for distfar is 
based on islandA only, with the interactions of islandx:distfar 
describes how being far on island x, y or z differs from being far on 
islandA. Here, there is no effect of distance on IslandA, but there is 
an effect of distance on IslandC relative to IslandA. It doesn't make a 
lot of sense to me to relate one island back to a reference island- it 
would make more sense to think of each island independently, and get 
coefficients for dist (regardless of island), dist*islandA, 
dist*islandB... and so on with the comparisons to a global mean effect 
of distance.
 >
 > #Approach 2: suppress intercept by adding -1; contrasts for island 
are still treatment contrasts
 > contrasts(aglaiafen$island)<-contr.treatment
 > 
survm2<-glmer(cbind(numalive,numseedplant-numalive)~-1+island*dist+centavgopen*dist+(1|site), 
family=binomial, data=aglaiafen)
 > summary(survm2)
Generalized linear mixed model fit by maximum likelihood ['glmerMod']
Family: binomial ( logit )
Formula: cbind(numalive, numseedplant - numalive) ~ -1 + island * dist + 
centavgopen * dist + (1 | site)
Data: aglaiafen

AIC BIC logLik deviance
353.3382 382.2129 -165.6691 331.3382

Random effects:
Groups Name Variance Std.Dev.
site (Intercept) 0.2545 0.5044
Number of obs: 102, groups: site, 13

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
islandA 1.923429 0.316594 6.075 1.24e-09 ***
islandB 0.996259 0.359586 2.771 0.0056 **
islandC 1.775531 0.388569 4.569 4.89e-06 ***
islandD 3.846454 0.827065 4.651 3.31e-06 ***
distfar 0.518158 0.337779 1.534 0.1250
centavgopen 0.096612 0.046132 2.094 0.0362 *
island2:distfar -0.467552 0.459631 -1.017 0.3090
island3:distfar -0.954689 0.472887 -2.019 0.0435 *
island4:distfar -0.534941 1.092387 -0.490 0.6243
distfar:centavgopen -0.006516 0.067602 -0.096 0.9232
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 > confint(survm2, method="Wald")
2.5 % 97.5 %
islandA 1.302916127 2.54394152
islandB 0.291484255 1.70103419
islandC 1.013948776 2.53711283
islandD 2.225437148 5.46747149
distfar -0.143877567 1.18019287
centavgopen 0.006194755 0.18702884
island2:distfar -1.368411924 0.43330703
island3:distfar -1.881530349 -0.02784674
island4:distfar -2.675980039 1.60609838
distfar:centavgopen -0.139013810 0.12598121
 > #get similar qualitative results as survm1 above, but the 
coefficients for the main island effects differ because they are 
relative to zero, not relative to islandA. coefficients for distance, 
openness and interactions are the same as in approach 1.
 >
 > #Approach 3: change contrasts to simple contrasts
 > c<-contr.treatment(4)
 > my.coding<-matrix(rep(1/4, 12), ncol=3)
 > my.simple<-c-my.coding
 > contrasts(aglaiafen$island)<-my.simple
 >
 > 
survm3<-glmer(cbind(numalive,numseedplant-numalive)~island*dist+centavgopen*dist+(1|site), 
family=binomial, data=aglaiafen)
 > summary(survm3)
Generalized linear mixed model fit by maximum likelihood ['glmerMod']
Family: binomial ( logit )
Formula: cbind(numalive, numseedplant - numalive) ~ island * dist + 
centavgopen * dist + (1 | site)
Data: aglaiafen

AIC BIC logLik deviance
353.3382 382.2129 -165.6691 331.3382

Random effects:
Groups Name Variance Std.Dev.
site (Intercept) 0.2544 0.5044
Number of obs: 102, groups: site, 13

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.135430 0.261757 8.158 3.4e-16 ***
island2 -0.927143 0.480345 -1.930 0.0536 .
island3 -0.147931 0.494245 -0.299 0.7647
island4 1.923037 0.875906 2.195 0.0281 *
distfar 0.028877 0.301839 0.096 0.9238
centavgopen 0.096604 0.046131 2.094 0.0362 *
island2:distfar -0.467541 0.459629 -1.017 0.3091
island3:distfar -0.954656 0.472885 -2.019 0.0435 *
island4:distfar -0.534855 1.092395 -0.490 0.6244
distfar:centavgopen -0.006514 0.067601 -0.096 0.9232
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 > confint(survm3, method="Wald")
2.5 % 97.5 %
(Intercept) 1.622395665 2.64846473
island2 -1.868600916 0.01431566
island3 -1.116633278 0.82077127
island4 0.206292798 3.63978081
distfar -0.562716853 0.62047025
centavgopen 0.006189382 0.18701920
island2:distfar -1.368397627 0.43331609
island3:distfar -1.881494409 -0.02781798
island4:distfar -2.675909531 1.60619928
distfar:centavgopen -0.139009874 0.12598173
 > #get different coefficients for main effects of island and distance, 
but same values as both models above for openness and all interactions.

*Questions: *
 > # 1) Is one of these methods (or some other method) best for testing 
the effect of distance on survival relative to canopy openness and 
island? Can I conclude that, for this species, there is not a distance 
effect for island A, B or D, there isn't an interaction between distance 
and openness, but that there was lower seedling survival in far plots on 
island C? Is that lower relative only to near plots on islandC? Why 
don't the coefficients for the interactions change when contrasts are 
changed?

 > # 2) I have similar data for five other species, and I'd like to 
compare the magnitude of the distance effect between species. Would 
adding a species to the model (by adding species*island*distance and 
species*openness*distance) be advisable, or is it better to just analyze 
each species separately to avoid the challenges with interpreting 
three-way interactions, especially in glmer's where options for 
hypothesis testing are more limited than lm's or glm's? Again, I don't 
really want to know the effect of one species relative to a reference 
species, but simply compare the magnitude of the distance effect between 
species- how do I set up my model to do that? I can share the full 
dataset offline if that would help.

 > # 3) Is the use of confint(model, method="Wald") acceptable for 
testing this? I tried using bootMer a few different ways (e.g. 
bootMer(survm1, nsim=1000), using a couple different functions), but all 
bootstrap runs fail each time I've tried. Likelihood ratio tests can 
tell me whether a factor should be included in a model, but not 
differences between levels of factors (e.g. islands or species). In 
addition, results from LRT's are similar to those obtained using 
confint. Setting up the glht (multcomp package) for this, given all of 
the interactions and the continuous openness variable seems a bit 
overwhelming.

 > #Many thanks in advance for your help!

Haldre

-- 
*************
Huxley Fellow
Department of Ecology and Evolutionary Biology
MS170
Rice University
PO Box 1892
Houston, TX 77005

	[[alternative HTML version deleted]]

