PREPARE DATASET FOR RUNNING MODELS                     

#z.transform quantitative vars
AKG.data$z.stoneweight=as.vector(scale(AKG.data$Impactor.weight.kg))
AKG.data$z.dbh=as.vector(scale(AKG.data$dbh))

#dummy code and center cateogorical vars
AKG.data$ast.code=as.numeric(AKG.data$ASTspecies.yn==levels(AKG.data$ASTspecies.yn)[2])
AKG.data$reclevel.code=as.numeric(AKG.data$recording.level==levels(AKG.data$recording.level)[2])
AKG.data$stonetyp.code=as.numeric(AKG.data$Impactor.Type==levels(AKG.data$Impactor.Type)[2])
AKG.data$hollow.code=as.numeric(AKG.data$throw.target=="hollow")
AKG.data$trunk.code=as.numeric(AKG.data$throw.target=="trunk")

#for random slopes center categorical vars
AKG.data$ctr.ast=AKG.data$ast.code-mean(AKG.data$ast.code)
AKG.data$ctr.rec=AKG.data$reclevel.code-mean(AKG.data$reclevel.code)
AKG.data$ctr.stonetyp=AKG.data$stonetyp.code-mean(AKG.data$stonetyp.code)
AKG.data$ctr.hollow=AKG.data$hollow.code-mean(AKG.data$hollow.code)
AKG.data$ctr.trunk=AKG.data$trunk.code-mean(AKG.data$trunk.code)

###############################################################################
RESPONSE 1:   ATTACK TIME
###############################################################################
resp=as.vector(AKG.data$AttackTime*1000)  #multiply by 1000 for easier interpretation of estimates

akg.attack=lmer(resp~ ASTspecies.yn  + throw.target + z.dbh + z.stoneweight  + Impactor.Type + recording.level +
           (1+z.dbh+ ctr.hollow+ ctr.trunk + ctr.ast+ctr.stonetyp+ctr.rec||Impactor)+ 
           (1+z.dbh +z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ ctr.rec||Latin.name)+  
           (1+z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ctr.rec||TreeID),                
            data=AKG.data, REML=F)

null.akg.attack=lmer(resp~ z.stoneweight  + Impactor.Type + recording.level +
           (1+z.dbh+ ctr.hollow+ ctr.trunk + ctr.ast+ctr.stonetyp+ctr.rec||Impactor)+ 
           (1+z.dbh +z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ ctr.rec||Latin.name)+  
           (1+z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ctr.rec||TreeID),              
            data=AKG.data, REML=F)

#Likelihood ratio test (LRT) to compare the null with the full model 
anova(null.akg.attack, akg.attack, test="Chisq")
                 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
null.akg.attack 25 585.15 655.86 -267.58   535.15                             
akg.attack      29 574.32 656.34 -258.16   516.32 18.837      4  0.0008461 ***

#full model summary
summary(akg.attack)  
Number of obs: 125, groups:  TreeID, 27; Latin.name, 13; Impactor, 10
Fixed effects:
                      Estimate Std. Error t value
(Intercept)             4.2152     0.7466   5.646
ASTspecies.ynyes        0.9711     0.6071   1.599                           
throw.targethollow     -1.4017     0.6664  -2.104
throw.targettrunk      -2.0682     0.5543  -3.731
z.dbh                   0.4215     0.2639   1.597
z.stoneweight          -0.0428     0.1699  -0.252
Impactor.Typelaterite  -0.1040     0.3226  -0.322
recording.levelthree    0.2394     0.4759   0.503

#LRT for all fixed effects in the model
drop1(akg.attack, test="Chisq")  
                Df    AIC     LRT  Pr(Chi)   
<none>             574.32                    
ASTspecies.yn    1 575.18  2.8642 0.090573 .                          
throw.target     2 581.18 10.8586 0.004386 **
z.dbh            1 574.60  2.2787 0.131159   
z.stoneweight    1 572.38  0.0628 0.802057   
Impactor.Type    1 572.42  0.1008 0.750893   
recording.level  1 572.52  0.2066 0.649474   

POSTHOC COMPARISON   #done only if one of the two categorical predictors of interest shows an LRT <0.05 in the full model

library(multcomp)
summary(glht(akg.attack, mcp(throw.target="Tukey"), covariate_average=TRUE), test=adjusted(type="none"))
Linear Hypotheses:
                       Estimate Std. Error z value Pr(>|z|)    
hollow - buttress == 0  -1.4017     0.6664  -2.104  0.03542 *  
trunk - buttress == 0   -2.0682     0.5543  -3.731  0.00019 ***
trunk - hollow == 0     -0.6665     0.6164  -1.081  0.27962    

###############################################################################
RESPONSE 2: SPECTRAL CENTROID (Repeat same as above with this response)
###############################################################################

akg.centroid=lmer(SpectralCentroid~ ASTspecies.yn  + throw.target + z.dbh + z.stoneweight  + Impactor.Type + recording.level +
           (1+z.dbh+ ctr.hollow+ ctr.trunk + ctr.ast+ctr.stonetyp+ctr.rec||Impactor)+ 
           (1+z.dbh +z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ ctr.rec||Latin.name)+  
           (1+z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ctr.rec||TreeID),                
            data=AKG.data, REML=F)

###############################################################################
RESPONSE 3: DAMPING COEFFICIENT  (Repeat same as above for this response)
###############################################################################
AKG.data$NNDamping=(AKG.data$Damping)*-1  #turn non-negative fo reasier interpretation of estimates and plots

akg.NNdamp=lmer(NNDamping~ ASTspecies.yn  + throw.target + z.dbh + z.stoneweight  + Impactor.Type + recording.level +
           (1+z.dbh+ ctr.hollow+ ctr.trunk + ctr.ast+ctr.stonetyp+ctr.rec||Impactor)+ 
           (1+z.dbh +z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ ctr.rec||Latin.name)+  
           (1+z.stoneweight+ctr.hollow+ctr.trunk+ctr.stonetyp+ctr.rec||TreeID),                
            data=AKG.data, REML=F)