4
$\begingroup$

After looking through several questions on stackexchange, and even applying many of the suggested methods I have come to a point where I would like to have an expert's advice if the methods applied on my dataset are correct. This is for my master-dissertation and up untill 1,5 months ago I never worked with R.

I studied the micro-habitat characteristics for the ovipostion of the Pyrgus Malvae butterfly. I measured several environmental variables as the percentage of the hostplant site area, more specifically with an example: the percentage of dwarf-shrub cover per host-plant site, or the percentage of wild-boar digging per hostplant site. Each site equals a circle with a radius of 25 cm with the hostplant in the center. For each occupied hostplant I selected a paired unoccupied hostplant.

The test I chose to use is a GLMM (glmer from the package lme4) since I want to account for several random effects such as the hostplant species (HP_spp), the date on which I measured these variables (VS_Date) and the pair number (Pair_nr) of each occupied and unoccupied hostplant pair. In the GLMM I want to use my environmental variables as the response variable and use Occupancy (0 = unoccupied = no egg found; 1 = occupied) as the independent variable to check for differences within each environmental variable between the two levels of occupancy.

I have performed 3 different kind of glmer: I will use one of my variables as an example here: VS_G = vegetation structure grass = the percentage of grass cover within a hostplant site.

1) standard glmer with family = binomial
Note: VS_G is expressed as a decimal (i.e. 0.7 for 70% grass cover)

GLMMS106_VS_G_Occ <- glmer(VS_G~Occupancy + 
                             (1|VS_Date) + (1|Pair_nr) + (1|HP_spp), family = binomial, data = PM_data106)

with this test I obviously get the error message:

Warning message: In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

And a resulting p-value of:

#R              Estimate Std. Error z   value Pr(>|z|)    
#R (Intercept)    -0.665        0.290  -2.293 0.021869 *  
#R Occupancy      -2.574        0.777  -3.312 0.000925 ***

2) After looking up this error I found out I had to change my test since I am working with proportional data, so I applied 2 tests with the cbind command and the weights command:

2a) glmer with cbind
Note: VS_G is expressed as a integer (i.e. 70 for 70% grass cover), I created VS_G_inv = 100-VS_G; so in case of VS_G == 70, then VS_G_inv == 30.

GLMMS106_VS_G_Occ <- glmer(cbind(VS_G, VS_G_inv)~Occupancy +   
                              (1|VS_Date) + (1|Pair_nr) + (1|HP_spp), family = binomial, data = PM_data106)

Hooray! no error anymore however the p-value seems of:

#R             Estimate Std. Error z value Pr(>|z|)    
#R (Intercept) -0.58302    0.22478  -2.594  0.00949 **   
#R Occupancy   -0.96390    0.04594 -20.982  < 2e-16 ***

2b) glmer with weights
Note: VS_G is expressed as an integer (i.e. 70 for 70% grass cover) and VS_Weights was created (= 100 for each site)

GLMMS106_VS_G_Occ <- glmer(VS_G/VS_Weights~Occupancy +   
                              (1|VS_Date) + (1|Pair_nr) + (1|HP_spp), family = binomial, data = PM_data106, weights = VS_Weights)

Same result as the cbind test:

#R             Estimate Std. Error z value Pr(>|z|)    
#R (Intercept) -0.58302    0.22478  -2.594  0.00949 **  
#R Occupancy   -0.96390    0.04594 -20.982  < 2e-16 ***

So after discussing this with my supervisors I found out that via using the cbind or weights command my sample size is increased by a factor 100 in this case. This in it's turn results in these very low p-values. So now I am wondering if the first test is ok to use after all?

$\endgroup$
  • 4
    $\begingroup$ One thing people often do with outcomes which are proportions is to use beta regression. I am not an expert here and I do not know if R has a package letting you fit it with random effects but I am sure there must be one. $\endgroup$ – mdewey Dec 15 '16 at 14:26
  • $\begingroup$ I considered using the beta regression through the glmmADMB-package (only glmm package with this option to my knowledge), but far to many errors regarding my random effect to get it to work properly. An option I am considering is converting the proportional data to square centimeters, since each hostplant site has a constant area of 1963.5 cm2 (circle with r=25 cm), and then using a Poisson distribution since this is often done with continuous variables in glmm's. stats.stackexchange.com/questions/38530/… $\endgroup$ – Frederic de Schaetzen Dec 15 '16 at 15:18
  • $\begingroup$ Both 1) and 2a) and 2b) are wrong. The issue is similar to this question. Take the example you provide. In 1) you estimate a model where you observe a 70% head outcome with one flip of a coin. The likelihood glm uses will be wrong. In 2a) and 2b) you estimate a model as if there where 100 coins flips 70 of which had a head. The latter "fake" 100 observation yields the much lower p-values. You only have on observations so this is clearly wrong. $\endgroup$ – Benjamin Christoffersen Aug 5 '18 at 15:24
1
$\begingroup$

I think that You might use simply a discrete binomial GLMM rather than continuous, which is just slightly different from the model You described. By the way, the warning message You've mentioned is no error: it simply notifies You, that the binomial response variable was continuous (having non-integer, i.e. non-zero and non-one values). If the non-integer values are between zero and one, there should be no problem.

If I understand correctly, You are interested in whether the structural vegetation coverage influences the presence of P. malvae eggs, correct? In that case, Occupancy should be the response variable, because You expect change in the presence of eggs in response to other environmental variables. Having the cause and causation in the right order helps to make sense of the results of such ecological models. In my opinion, the model You might want to use would look something like this:

glmer(Occupancy~VS_G+HP_spp+(1|VS_Date)+(1|Pair_nr), family="binomial", data=PM_data106)

In this model, You can add all 3 variables representing vegetation-coverages, although You should keep an eye out for the possible interactions between the three coverage variables. Also, i think it would make more ecological sense, to specify the hostplant species as fixed effect, because it has relevance to see, whether one plant species or another is preferred by the animal for laying eggs.

Cheers,

ZR

$\endgroup$
  • $\begingroup$ Thank you for this indept answer! I already tested the effect of the different variables on the presence of P. malvae eggs the way you described it. I should have been more specific in my description but I sampled in three distinct habitats. And now I want to test whether there are differences between the different levels of Occupancy (Occupied, Unoccupied), species (Potentilla erecta and Potentilla anglica), habitat plots (Heathland, mixed heath/grassland and grassland) and the interactions between the three of them. In the end I want to do eight batches of tests: $\endgroup$ – Frederic de Schaetzen Dec 16 '16 at 11:37
  • $\begingroup$ Env_var~Occupancy; Env_var~HP_spp; Env_var~Plot; Env_var~Occupancy x Plot; Env_var~Occupancy x Species; Env_var~Plot x Species and Env_var~Occupancy x Plot x Species. And the only way to do this, by my knowledge, is by putting the environmental variable as the response variable. $\endgroup$ – Frederic de Schaetzen Dec 16 '16 at 11:38
  • $\begingroup$ Furthermore the preference of plant species by P. malvae in this study is related to the habitat type. P. erecta was only available in the Heathland habitat since it had the capability of growing in between the dwarf-shrubs (Calluna vulgaris). So in the end I might even add HP_spp as a random effect since it correlated with the different levels of the habitat plots. $\endgroup$ – Frederic de Schaetzen Dec 16 '16 at 11:38
  • $\begingroup$ Is the variable "Plot" a discrete (factor) variable? If so, then it might be a more simple but still quite reliable approach to use Conover-Iman tests (cran.r-project.org/web/packages/conover.test/index.html), where the variable in test is "Env_var", and the groups are the factor variables You described (each in a different conover test). This will test whether there are group-differences, and what are the direction of differences. However, You should post hoc adjust p-values resulted from the tests (e.g. with p.adjust(x, method="bonferroni"). Hope this helps. $\endgroup$ – Z. Radai Dec 16 '16 at 12:14
  • 3
    $\begingroup$ Thank you for the suggestions Z. Radai! However, I think I have found the most appropriate way of testing my data. I could finally get the glmmADMB package running without any errors, and therefore I can use a beta distribution glmm. In my study the inclusion of random factors is rather important, especially since the pair-wise setup of my sampling (For each occupied hostplant, I selected an unoccupied hostplant) so a mixed model is prefered. $\endgroup$ – Frederic de Schaetzen Dec 16 '16 at 18:03

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.