3
$\begingroup$

Let's say I am looking at how unicorn herd size changes with food quantity. Unicorn herd sizes were surveyed at seven localities over the course of twelve months. Food quantity was assessed monthly, so for each locality at each month we have one value of food quantity (continuous predictor, let's say grass cover). For each locality on each month, unicorn herds were re-surveyed over four consecutive days to account for both unicorn movement and possible sampling error. On each sampling day there is a variable number of data points, each referring to the size of a different herd. To measure the correlation between herd size and food quantity I also need to account for:

  • the random effect of locality
  • the non-independence of measures done on different days of each survey month at each locality

(The random effect of survey month is negligible and I will not include it in the model).

Given the simulated dataset:

unicorns <- read.csv(text= RCurl::getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/unicorns.csv"))

I can model the data using glmmTMB:

unicorns$survey_day_f <- glmmTMB::numFactor(unicorns$survey_day_n)
levels(unicorns$survey_day_f)

unicorns_glmmTMB_ou <- glmmTMB(Herd_size_n ~ food.quantity        # fixed effect
                        + (1 | Locality)                          # random effect of locality
                        + ou(survey_day_f + 0 | Locality),        # correction for non-independence
                        family="poisson",
                        data=unicorns)
unicorns_glmmTMB_ou$sdr$pdHess==T # the optimiser converges
summary(unicorns_glmmTMB_ou)$coef

#                 Estimate Std. Error   z value    Pr(>|z|)
# (Intercept)   -0.1253359  0.7874638 -0.159164 0.873539653
# food.quantity  0.7145391  0.2730064  2.617298 0.008862902

An alternative approach came to my mind to account for the temporal autocorrelation of repeated measures. I subset the dataset by picking the data of a single survey day at random for each locality in each survey month and analysed the resulting subset of the original dataset as a study in itself. I repeated this operation N times and fitted on each of the N datasets the model:

glmmTMB1 <- glmmTMB(Herd_size_n ~ food.quantity
                        + (1 | Locality),
                        family="poisson",
                        data= subset)

Finally I averaged the N slope estimates, their N estimated standard errors, and the N associated p-values to obtain a single mean estimate.

Here is the code, for reproducibility, with N=200:

# create empty dataframes to store results:
unicorn.results <- data.frame(fit.n=rep(NA, imax),
                            slope=rep(NA, imax),
                            SE=rep(NA, imax),
                            pval=rep(NA, imax)
                            )
unicorn.results.mean <- data.frame(N=rep(NA, imax),
                            mean.slope=rep(NA, imax),
                            mean.SE=rep(NA, imax),
                            mean.pval=rep(NA, imax))
library(dplyr)
imax=200 # N of iterations
unicorns$DateLocality <- paste(unicorns$Date, unicorns$Locality)
for(i in 1:imax){
# subsampling unicorns by picking a day at random from each Year_Month*Locality:
new_df <- as.data.frame(unicorns %>% group_by(DateLocality) %>% sample_n(1))	
glmmTMB1 <- glmmTMB(Herd_size_n ~ food.quantity
						+ (1 | Locality),
						family="poisson",
						data= new_df)
if (glmmTMB1$sdr$pdHess==T){ # save estimates (only IF/WHEN they converge)
	unicorn.results$fit.n[i]=i
	unicorn.results$slope[i]=summary(glmmTMB1)$coef[[1]][2,1]
	unicorn.results$SE[i]=summary(glmmTMB1)$coef[[1]][2,2]
	unicorn.results$pval[i]=summary(glmmTMB1)$coef[[1]][2,4]
    
    unicorn.results.mean$N[i]=i
	unicorn.results.mean$mean.slope[i]=mean(unicorn.results$slope[1:i])
	unicorn.results.mean$mean.SE[i]=mean(unicorn.results$SE[1:i])
	unicorn.results.mean$mean.pval[i]=mean(unicorn.results$pval[1:i])
}
}

enter image description here

> unicorn.results.mean[200,]
      N mean.slope    mean.SE   mean.pval
200 200  0.6386316 0.08494007 9.80296e-06

For comparison, here are the estimates from unicorns_glmmTMB_ou again:

summary(unicorns_glmmTMB_ou)$coef

#                 Estimate Std. Error   z value    Pr(>|z|)
# (Intercept)   -0.1253359  0.7874638 -0.159164 0.873539653
# food.quantity  0.7145391  0.2730064  2.617298 0.008862902

What caused these discrepancies? My guess is that resampling the original dataset 200 times produced an unwarranted, artificial increase in the power of the analysis, which might explain the narrower StdErr and smaller p-value. Any more accurate and less hand-waving explanation is very welcome. Moreover I don't understand why the slope estimates obtained with the two approaches do not match.

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.