## Example of a step in the model selection procedure #### mod.NB.GC <- glm.nb(Non ~ offset(I(1*log(L))) + log(gc), data = data) # model WITH gc included mod.NB.0 <- glm.nb(Non ~ offset(I(1*log(L))), data = data) # model WITHout gc included require(lmtest) lrtest(mod.NB.GC, mod.NB.0) # likelihood ratio test comparing the two models # If lrtest is significant, keep dS in the model, otherwise discard. ## Estimating the P-value for significance of a predictor variable #### N = 1000 # number of permutations LLrt = NULL ## In this example, we estimate the P-value for %GC M <- glm.nb(Non ~ offset(I(1*log(L))) + log(gc), data = data) # model WITH gc M0 <- glm.nb(Non ~ offset(I(1*log(L))), data = data) # model WITHOUT gc LLrt.real = -2*(logLik(M0)[1]-logLik(M)[1]) # calculate the real LL-ratio for(i in 1:N){ data$permute.gc = sample(data$gc, size = length(data$gc), replace = F) # permute the gc data permute.M <- glm.nb(Non ~ offset(I(1*log(L))) + log(permute.gc), data = data) # re-run the model LLrt[i] = -2*(logLik(M0)[1]-logLik(permute.M)[1]) # collect the LL-ratio statistic for the permuted data } hist(LLrt) # look at the generated null distribution P.value = sum(abs(LLrt)>abs(LLrt.real))/N # calculate the P-value