NonDirected Activity
Hannah & Marcus
Thursday, February 19, 2015
NON DIRECTED ACTIVTY TEST
- Is there a difference between urban and non-urban populations
- Is there a difference between populations
- Is activity level repeatable
- Does behavioural variation vary between populations (i.e are skinks more consistent in their behaviour in different populations)
Populations are nested within regions and skinks are the replicates, with repeated nested within these replicates
Testing Assumptions
library(car)
library(plyr)
library(multcomp)
library(nlme)
library(reshape)
library(MASS)
library(Hmisc)
library(Rmisc)
library(ggplot2)
library(gridExtra)
library(lme4)
setwd("//ad.monash.edu/home/User081/mmichel/Desktop/Urban paper")
activity<-read.csv("NonDirected.csv", header=T, strip.white=T)
head(activity)## Repeat Region Population VIE Trans LogTrans SqrtTrans
## 1 A N NP 1 RR.. 47 1.6720979 6.855655
## 2 A N NP 1 OO.. 20 1.3010300 4.472136
## 3 A N NP 1 YY.. 24 1.3802112 4.898979
## 4 A N NP 1 ..RR 9 0.9542425 3.000000
## 5 A N NP 1 ..OO 111 2.0453230 10.535654
## 6 A N NP 1 ..YY 60 1.7781513 7.745967str(activity)## 'data.frame': 118 obs. of 7 variables:
## $ Repeat : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
## $ Region : Factor w/ 2 levels "N","U": 1 1 1 1 1 1 1 1 1 1 ...
## $ Population: Factor w/ 4 levels "NP 1","NP 2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ VIE : Factor w/ 59 levels "..OO","..OR",..: 42 27 56 4 1 8 39 34 53 46 ...
## $ Trans : int 47 20 24 9 111 60 56 92 134 47 ...
## $ LogTrans : num 1.672 1.301 1.38 0.954 2.045 ...
## $ SqrtTrans : num 6.86 4.47 4.9 3 10.54 ...activity.ag<-ddply(activity, ~Region+Population, numcolwise(mean))
activity.ag## Region Population Trans LogTrans SqrtTrans
## 1 N NP 1 67.41176 1.674703 7.640612
## 2 N NP 2 70.10714 1.741573 7.921416
## 3 U UP 1 62.85294 1.669405 7.439638
## 4 U UP 2 100.95455 1.911891 9.598595#Means for each population, obviously something funny is happening with UP 2 across both trials
ggplot(activity.ag, aes(y=Trans, x=Region))+geom_boxplot()
#Shows that the urban populations are highly varied, whereas the non urban are both consistentURBAN vs. NON URBAN
activity.null<-lme(Trans~1, random=~1|Region/Population/VIE/Repeat, data=activity, method='ML')
activity.lme<-lme(Trans~Region, random=~1|Region/Population/VIE/Repeat, data=activity, method='ML')
anova(activity.null, activity.lme)## Model df AIC BIC logLik Test L.Ratio
## activity.null 1 6 1242.883 1259.507 -615.4417
## activity.lme 2 7 1244.216 1263.611 -615.1080 1 vs 2 0.6674704
## p-value
## activity.null
## activity.lme 0.4139AIC(activity.lme)## [1] 1244.216AIC(activity.null)## [1] 1242.883#region has no impactattach(activity)
activity.ag2<-ddply(activity, ~Population+VIE, numcolwise(mean))
activity.ag2## Population VIE Trans LogTrans SqrtTrans
## 1 NP 1 ..OO 123.5 2.0894309 11.098779
## 2 NP 1 ..OR 64.0 1.6758015 7.463328
## 3 NP 1 ..OY 39.0 1.5611079 6.140055
## 4 NP 1 ..RR 36.0 1.3767915 5.468627
## 5 NP 1 ..RY 96.0 1.9493626 9.617546
## 6 NP 1 ..YO 17.5 0.7720340 2.958040
## 7 NP 1 ..YR 49.5 1.6690289 6.934301
## 8 NP 1 ..YY 67.5 1.8266063 8.203110
## 9 NP 1 O..O 48.5 1.4493626 6.190416
## 10 NP 1 OO.. 61.0 1.6548151 7.285820
## 11 NP 1 OY.. 94.0 1.9730295 9.694811
## 12 NP 1 R.Y. 34.0 1.2491553 5.086661
## 13 NP 1 RO.. 94.5 1.9360198 9.507939
## 14 NP 1 RR.. 69.5 1.8179428 8.223659
## 15 NP 1 RY.. 74.5 1.8403490 8.477580
## 16 NP 1 YO.. 147.5 2.1669653 12.132207
## 17 NP 1 YY.. 29.5 1.4621396 5.407530
## 18 NP 2 .O.O 33.0 1.4929377 5.661833
## 19 NP 2 .O.R 47.0 1.6237411 6.672615
## 20 NP 2 .OO. 30.0 1.4237863 5.316625
## 21 NP 2 .R.O 95.0 1.9440617 9.563334
## 22 NP 2 .R.R 125.0 2.0933228 11.157324
## 23 NP 2 .RR. 110.5 1.9628623 10.057439
## 24 NP 2 .Y.Y 45.5 1.5836587 6.474588
## 25 NP 2 O..R 28.5 1.4466034 5.313392
## 26 NP 2 O.O. 67.0 1.7241986 7.745662
## 27 NP 2 OOOO 68.5 1.8257359 8.229451
## 28 NP 2 R..R 60.5 1.7722820 7.736103
## 29 NP 2 RRRR 108.0 2.0297438 10.370360
## 30 NP 2 Y..Y 87.5 1.5871753 7.942049
## 31 NP 2 Y.Y. 75.5 1.8719116 8.659043
## 32 UP 1 .OR. 152.0 2.1807034 12.320744
## 33 UP 1 .OY. 77.0 1.8861609 8.773299
## 34 UP 1 .R.Y 51.5 1.6963485 7.113336
## 35 UP 1 .RO. 66.5 1.8183442 8.133816
## 36 UP 1 .Y.R 31.5 1.4542425 5.475422
## 37 UP 1 .YY. 68.0 1.8317562 8.242641
## 38 UP 1 O..Y 80.0 1.8747135 8.801704
## 39 UP 1 OORR 39.5 1.5491488 6.120125
## 40 UP 1 OOYY 24.5 1.3883506 4.947426
## 41 UP 1 OROR 54.0 1.6751240 7.117887
## 42 UP 1 ORRO 17.5 1.2108020 4.107802
## 43 UP 1 R.O. 64.0 1.8048506 7.993885
## 44 UP 1 ROOR 54.5 1.7332855 7.369226
## 45 UP 1 RROO 17.5 1.2204545 4.129967
## 46 UP 1 Y..O 115.0 2.0583205 10.709160
## 47 UP 1 YYOO 38.0 0.9404068 4.358899
## 48 UP 1 YYYY 117.5 2.0568714 10.758514
## 49 UP 2 .YO. 97.0 1.9867717 9.848858
## 50 UP 2 OOOY 80.0 1.8441654 8.655918
## 51 UP 2 OYYO 120.0 2.0631631 10.854830
## 52 UP 2 RORO 40.5 1.6058272 6.358006
## 53 UP 2 RRRY 80.0 1.8598313 8.729704
## 54 UP 2 RYRY 160.5 2.2039333 12.657631
## 55 UP 2 Y..R 81.5 1.8422431 8.690263
## 56 UP 2 Y.O. 192.0 2.2828236 13.852598
## 57 UP 2 YOOY 15.0 1.1721961 3.864328
## 58 UP 2 YRYR 115.5 2.0618932 10.742834
## 59 UP 2 YYYR 128.5 2.1079509 11.329576g0<-ggplot(activity.ag2, aes(y=Trans, x=Population))+geom_boxplot()
g0
g1<-ggplot(activity, aes(y=Trans, x=Repeat, color=Population))+scale_color_discrete(guide=FALSE)+geom_boxplot()+theme_classic()
g1
tapply(Trans, list(Population, Repeat), mean)## A B
## NP 1 51.00000 83.82353
## NP 2 82.78571 57.42857
## UP 1 59.88235 65.82353
## UP 2 98.81818 103.09091#Interesting that National park looked more variable between trials
model<-glmer(Trans~Population+(1|Region/Population/VIE/Repeat), family=poisson)
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population + (1 | Region/Population/VIE/Repeat)
##
## AIC BIC logLik deviance df.resid
## 1259.6 1281.8 -621.8 1243.6 110
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.04306 -0.13941 0.04127 0.08181 0.18023
##
## Random effects:
## Groups Name Variance Std.Dev.
## Repeat:(VIE:(Population:Region)) (Intercept) 5.029e-01 0.7091473
## VIE:(Population:Region) (Intercept) 1.370e-01 0.3701208
## Population:Region (Intercept) 0.000e+00 0.0000000
## Region (Intercept) 1.124e-10 0.0000106
## Number of obs: 118, groups:
## Repeat:(VIE:(Population:Region)), 118; VIE:(Population:Region), 59; Population:Region, 4; Region, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.91528 0.15415 25.399 <2e-16 ***
## PopulationNP 2 0.10556 0.22863 0.462 0.6443
## PopulationUP 1 -0.02848 0.21768 -0.131 0.8959
## PopulationUP 2 0.49418 0.24450 2.021 0.0433 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1
## PopulatnNP2 -0.674
## PopulatnUP1 -0.707 0.477
## PopulatnUP2 -0.630 0.425 0.446#Population:Region is the fixed effect (seen under random effects)
#Most variance is seen between repeats within skinks
#Second most variance between skinks...i.e. there is more within-individual variance than between individual variance. However when we look at means, its seems that this variance exists mainly in the national park populations, but not the urban populations
model.null<-glmer(Trans~1+(1|Region/Population/VIE/Repeat), family=poisson)
summary(model.null)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ 1 + (1 | Region/Population/VIE/Repeat)
##
## AIC BIC logLik deviance df.resid
## 1258.8 1272.6 -624.4 1248.8 113
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05110 -0.11735 0.03779 0.09120 0.18184
##
## Random effects:
## Groups Name Variance Std.Dev.
## Repeat:(VIE:(Population:Region)) (Intercept) 0.5025 0.70886
## VIE:(Population:Region) (Intercept) 0.1725 0.41539
## Population:Region (Intercept) 0.0023 0.04796
## Region (Intercept) 0.0000 0.00000
## Number of obs: 118, groups:
## Repeat:(VIE:(Population:Region)), 118; VIE:(Population:Region), 59; Population:Region, 4; Region, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.02619 0.09277 43.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Value for population as a random effect is very small, so should be considred a fixed effect, as well as region.
model3<-glmer(Trans~Population*Repeat + (1|VIE), data=activity, family=poisson)
summary(model3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population * Repeat + (1 | VIE)
## Data: activity
##
## AIC BIC logLik deviance df.resid
## 2027.7 2052.6 -1004.8 2009.7 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.2470 -1.5157 -0.1769 1.5146 11.2239
##
## Random effects:
## Groups Name Variance Std.Dev.
## VIE (Intercept) 0.3115 0.5581
## Number of obs: 118, groups: VIE, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.80208 0.14000 27.157 < 2e-16 ***
## PopulationNP 2 0.51012 0.20694 2.465 0.01370 *
## PopulationUP 1 0.12478 0.19768 0.631 0.52788
## PopulationUP 2 0.63112 0.22148 2.850 0.00438 **
## RepeatB 0.49689 0.04300 11.556 < 2e-16 ***
## PopulationNP 2:RepeatB -0.86260 0.06283 -13.728 < 2e-16 ***
## PopulationUP 1:RepeatB -0.40229 0.06097 -6.598 4.17e-11 ***
## PopulationUP 2:RepeatB -0.45456 0.06038 -7.528 5.16e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1 PplUP2 RepetB PNP2:R PUP1:R
## PopulatnNP2 -0.676
## PopulatnUP1 -0.708 0.479
## PopulatnUP2 -0.632 0.428 0.448
## RepeatB -0.191 0.129 0.135 0.121
## PpltnNP2:RB 0.131 -0.155 -0.093 -0.083 -0.684
## PpltnUP1:RB 0.135 -0.091 -0.177 -0.085 -0.705 0.483
## PpltnUP2:RB 0.136 -0.092 -0.096 -0.155 -0.712 0.487 0.502AIC(model3)## [1] 2027.666plot(model3)
#Probably need to split analysis by trialRepeatA <- subset(activity, Repeat=="A")
RepeatA## Repeat Region Population VIE Trans LogTrans SqrtTrans
## 1 A N NP 1 RR.. 47 1.6720979 6.855655
## 2 A N NP 1 OO.. 20 1.3010300 4.472136
## 3 A N NP 1 YY.. 24 1.3802112 4.898979
## 4 A N NP 1 ..RR 9 0.9542425 3.000000
## 5 A N NP 1 ..OO 111 2.0453230 10.535654
## 6 A N NP 1 ..YY 60 1.7781513 7.745967
## 7 A N NP 1 RO.. 56 1.7481880 7.483315
## 8 A N NP 1 OY.. 92 1.9637878 9.591663
## 9 A N NP 1 YO.. 134 2.1271048 11.575837
## 10 A N NP 1 RY.. 47 1.6720979 6.855655
## 11 A N NP 1 ..RY 60 1.7781513 7.745967
## 12 A N NP 1 ..OR 107 2.0293838 10.344080
## 13 A N NP 1 ..YO 0 0.0000000 0.000000
## 14 A N NP 1 ..OY 53 1.7242759 7.280110
## 15 A N NP 1 ..YR 33 1.5185139 5.744563
## 16 A N NP 1 O..O 9 0.9542425 3.000000
## 17 A N NP 1 R.Y. 5 0.6989700 2.236068
## 18 A N NP 2 Y.Y. 88 1.9444827 9.380832
## 19 A N NP 2 O.O. 108 2.0334238 10.392305
## 20 A N NP 2 .R.R 141 2.1492191 11.874342
## 21 A N NP 2 .O.O 44 1.6434527 6.633250
## 22 A N NP 2 Y..Y 9 0.9542425 3.000000
## 23 A N NP 2 RRRR 122 2.0863598 11.045361
## 24 A N NP 2 OOOO 83 1.9190781 9.110434
## 25 A N NP 2 .O.R 68 1.8325089 8.246211
## 26 A N NP 2 R..R 73 1.8633229 8.544004
## 27 A N NP 2 .Y.Y 70 1.8450980 8.366600
## 28 A N NP 2 .OO. 16 1.2041200 4.000000
## 29 A N NP 2 .RR. 172 2.2355284 13.114877
## 30 A N NP 2 O..R 34 1.5314789 5.830952
## 31 A N NP 2 .R.O 131 2.1172713 11.445523
## 32 A U UP 1 O..Y 52 1.7160033 7.211103
## 33 A U UP 1 Y..O 103 2.0128372 10.148892
## 34 A U UP 1 R.O. 69 1.8388491 8.306624
## 35 A U UP 1 .R.Y 38 1.5797836 6.164414
## 36 A U UP 1 .OY. 80 1.9030900 8.944272
## 37 A U UP 1 YYYY 89 1.9493900 9.433981
## 38 A U UP 1 OORR 22 1.3424227 4.690416
## 39 A U UP 1 .RO. 76 1.8808136 8.717798
## 40 A U UP 1 ROOR 48 1.6812412 6.928203
## 41 A U UP 1 .Y.R 45 1.6532125 6.708204
## 42 A U UP 1 RROO 23 1.3617278 4.795832
## 43 A U UP 1 OROR 28 1.4471580 5.291503
## 44 A U UP 1 .YY. 72 1.8573325 8.485281
## 45 A U UP 1 .OR. 163 2.2121876 12.767145
## 46 A U UP 1 OOYY 23 1.3617278 4.795832
## 47 A U UP 1 YYOO 76 1.8808136 8.717798
## 48 A U UP 1 ORRO 11 1.0413927 3.316625
## 49 A U UP 2 .YO. 97 1.9867717 9.848858
## 50 A U UP 2 OOOY 41 1.6127839 6.403124
## 51 A U UP 2 RYRY 147 2.1673173 12.124356
## 52 A U UP 2 YYYR 120 2.0791812 10.954451
## 53 A U UP 2 OYYO 152 2.1818436 12.328828
## 54 A U UP 2 YOOY 17 1.2304489 4.123106
## 55 A U UP 2 RORO 37 1.5682017 6.082763
## 56 A U UP 2 YRYR 122 2.0863598 11.045361
## 57 A U UP 2 Y..R 39 1.5910646 6.244998
## 58 A U UP 2 Y.O. 201 2.3031961 14.177447
## 59 A U UP 2 RRRY 114 2.0569049 10.677078RepeatB <- subset(activity, Repeat=="B")
#Here I just created subsets of the data so that I could focus on each trial seperately
glmer.model<-glmer(Trans~Population + (1|VIE), data=RepeatA, family=poisson)
plot(glmer.model)
summary(glmer.model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population + (1 | VIE)
## Data: RepeatA
##
## AIC BIC logLik deviance df.resid
## 626.8 637.2 -308.4 616.8 54
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.83446 -0.15072 0.04455 0.08392 0.17243
##
## Random effects:
## Groups Name Variance Std.Dev.
## VIE (Intercept) 0.6863 0.8284
## Number of obs: 59, groups: VIE, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5232 0.2082 16.920 < 2e-16 ***
## PopulationNP 2 0.6633 0.3062 2.166 0.03028 *
## PopulationUP 1 0.3715 0.2917 1.274 0.20284
## PopulationUP 2 0.8507 0.3273 2.599 0.00936 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1
## PopulatnNP2 -0.679
## PopulatnUP1 -0.713 0.484
## PopulatnUP2 -0.635 0.432 0.453#Large signficance with repeat A
summary(glht(glmer.model, linfct=mcp(Population='Tukey')))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmer(formula = Trans ~ Population + (1 | VIE), data = RepeatA,
## family = poisson)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## NP 2 - NP 1 == 0 0.6633 0.3062 2.166 0.1324
## UP 1 - NP 1 == 0 0.3715 0.2917 1.274 0.5790
## UP 2 - NP 1 == 0 0.8507 0.3273 2.599 0.0459 *
## UP 1 - NP 2 == 0 -0.2918 0.3038 -0.961 0.7713
## UP 2 - NP 2 == 0 0.1874 0.3382 0.554 0.9453
## UP 2 - UP 1 == 0 0.4792 0.3251 1.474 0.4523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)library(phia)
testInteractions(glmer.model)## Chisq Test:
## P-value adjustment method: holm
## Value Df Chisq Pr(>Chisq)
## NP 1-NP 2 0.51516 1 4.6935 0.15138
## NP 1-UP 1 0.68973 1 1.6218 0.60852
## NP 1-UP 2 0.42712 1 6.7537 0.05613 .
## NP 2-UP 1 1.33887 1 0.9226 0.67357
## NP 2-UP 2 0.82911 1 0.3071 0.67357
## UP 1-UP 2 0.61926 1 2.1732 0.56172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1glmer.model2<-glmer(Trans~Population + (1|VIE), data=RepeatB, family=poisson)
summary(glmer.model2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population + (1 | VIE)
## Data: RepeatB
##
## AIC BIC logLik deviance df.resid
## 629.7 640.1 -309.9 619.7 54
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.14429 -0.19594 0.03942 0.08433 0.18925
##
## Random effects:
## Groups Name Variance Std.Dev.
## VIE (Intercept) 0.5142 0.7171
## Number of obs: 59, groups: VIE, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2834 0.1766 24.250 <2e-16 ***
## PopulationNP 2 -0.4308 0.2639 -1.632 0.103
## PopulationUP 1 -0.3933 0.2516 -1.563 0.118
## PopulationUP 2 0.1682 0.2816 0.597 0.550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1
## PopulatnNP2 -0.669
## PopulatnUP1 -0.701 0.470
## PopulatnUP2 -0.627 0.420 0.440lme.model<-lme(Trans~Population, random=~1|VIE, data=RepeatA, correlation=corAR1(form=~VIE), method='REML')
summary(lme.model)## Linear mixed-effects model fit by REML
## Data: RepeatA
## AIC BIC logLik
## 600.3822 614.4335 -293.1911
##
## Random effects:
## Formula: ~1 | VIE
## (Intercept) Residual
## StdDev: 42.46657 15.92496
##
## Correlation Structure: AR(1)
## Formula: ~VIE | VIE
## Parameter estimate(s):
## Phi
## 0
## Fixed effects: Trans ~ Population
## Value Std.Error DF t-value p-value
## (Intercept) 51.00000 11.00004 55 4.636348 0.0000
## PopulationNP 2 31.78571 16.36858 55 1.941874 0.0573
## PopulationUP 1 8.88235 15.55640 55 0.570977 0.5703
## PopulationUP 2 47.81818 17.54999 55 2.724685 0.0086
## Correlation:
## (Intr) PplNP2 PplUP1
## PopulationNP 2 -0.672
## PopulationUP 1 -0.707 0.475
## PopulationUP 2 -0.627 0.421 0.443
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.63341894 -0.28553532 -0.01407598 0.18733617 0.79831487
##
## Number of Observations: 59
## Number of Groups: 59lme.model2<-lme(Trans~Population, random=~1|VIE, data=RepeatB, correlation=corAR1(form=~VIE), method='REML')
summary(lme.model2)## Linear mixed-effects model fit by REML
## Data: RepeatB
## AIC BIC logLik
## 598.3508 612.4021 -292.1754
##
## Random effects:
## Formula: ~1 | VIE
## (Intercept) Residual
## StdDev: 41.68951 15.63356
##
## Correlation Structure: AR(1)
## Formula: ~VIE | VIE
## Parameter estimate(s):
## Phi
## 0
## Fixed effects: Trans ~ Population
## Value Std.Error DF t-value p-value
## (Intercept) 83.82353 10.79876 55 7.762332 0.0000
## PopulationNP 2 -26.39496 16.06906 55 -1.642595 0.1062
## PopulationUP 1 -18.00000 15.27175 55 -1.178647 0.2436
## PopulationUP 2 19.26738 17.22885 55 1.118320 0.2683
## Correlation:
## (Intr) PplNP2 PplUP1
## PopulationNP 2 -0.672
## PopulationUP 1 -0.707 0.475
## PopulationUP 2 -0.627 0.421 0.443
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.71046486 -0.25967762 -0.03803878 0.15411609 0.85620387
##
## Number of Observations: 59
## Number of Groups: 59#When looking at each trial seperately, UP2 differs from the other parks but only in Trial 1. In trial 2 there is no significant difference...not againest NP 1 anyway
glmer.model3<-glmer(Trans~Population+(1|VIE/Repeat), data=activity, family=poisson)
summary(glmer.model3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population + (1 | VIE/Repeat)
## Data: activity
##
## AIC BIC logLik deviance df.resid
## 1255.6 1272.2 -621.8 1243.6 112
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.04305 -0.13941 0.04127 0.08181 0.18022
##
## Random effects:
## Groups Name Variance Std.Dev.
## Repeat:VIE (Intercept) 0.5029 0.7092
## VIE (Intercept) 0.1370 0.3701
## Number of obs: 118, groups: Repeat:VIE, 118; VIE, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.91530 0.15415 25.399 <2e-16 ***
## PopulationNP 2 0.10553 0.22863 0.462 0.6444
## PopulationUP 1 -0.02851 0.21768 -0.131 0.8958
## PopulationUP 2 0.49415 0.24450 2.021 0.0433 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1
## PopulatnNP2 -0.674
## PopulatnUP1 -0.707 0.477
## PopulatnUP2 -0.630 0.425 0.446anova(glmer.model3)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Population 3 5.4164 1.8055 1.8055#Model withouth Repeat as a factor
glmer.model4<-glmer(Trans~Population*Repeat+(1|VIE/Repeat), data=activity, family=poisson)
summary(glmer.model4)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population * Repeat + (1 | VIE/Repeat)
## Data: activity
##
## AIC BIC logLik deviance df.resid
## 1252.0 1279.7 -616.0 1232.0 108
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15646 -0.13529 0.04344 0.08250 0.26810
##
## Random effects:
## Groups Name Variance Std.Dev.
## Repeat:VIE (Intercept) 0.4054 0.6367
## VIE (Intercept) 0.1853 0.4305
## Number of obs: 118, groups: Repeat:VIE, 118; VIE, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5386 0.1936 18.281 < 2e-16 ***
## PopulationNP 2 0.6558 0.2846 2.304 0.02121 *
## PopulationUP 1 0.3572 0.2713 1.317 0.18794
## PopulationUP 2 0.8350 0.3043 2.744 0.00607 **
## RepeatB 0.7435 0.2265 3.282 0.00103 **
## PopulationNP 2:RepeatB -1.0862 0.3352 -3.241 0.00119 **
## PopulationUP 1:RepeatB -0.7531 0.3200 -2.354 0.01859 *
## PopulationUP 2:RepeatB -0.6705 0.3576 -1.875 0.06080 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1 PplUP2 RepetB PNP2:R PUP1:R
## PopulatnNP2 -0.679
## PopulatnUP1 -0.713 0.484
## PopulatnUP2 -0.636 0.432 0.453
## RepeatB -0.605 0.411 0.432 0.385
## PpltnNP2:RB 0.409 -0.596 -0.291 -0.260 -0.676
## PpltnUP1:RB 0.430 -0.292 -0.597 -0.273 -0.709 0.479
## PpltnUP2:RB 0.384 -0.260 -0.273 -0.596 -0.633 0.428 0.449#This is the final model, with Repeat as a factor...as you can see Repeat:VIE shows the largest variance, whereas difference between skinks is less. This suggests that are repeatability isnt very good.
anova(glmer.model4)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Population 3 5.3226 1.7742 1.7742
## Repeat 1 1.3617 1.3617 1.3617
## Population:Repeat 3 11.4223 3.8074 3.8074anova(glmer.model3, glmer.model4)## Data: activity
## Models:
## glmer.model3: Trans ~ Population + (1 | VIE/Repeat)
## glmer.model4: Trans ~ Population * Repeat + (1 | VIE/Repeat)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## glmer.model3 6 1255.6 1272.2 -621.81 1243.6
## glmer.model4 10 1252.0 1279.7 -616.00 1232.0 11.624 4 0.02037 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Compare the two models, glmer.model4 clearly the stronger one. Essentially alot variation is explained by repeat
library(phia)
testInteractions(glmer.model4)## Chisq Test:
## P-value adjustment method: holm
## Value Df Chisq Pr(>Chisq)
## NP 1-NP 2 : A-B 0.33748 1 10.5014 0.007157 **
## NP 1-UP 1 : A-B 0.47091 1 5.5393 0.092972 .
## NP 1-UP 2 : A-B 0.51144 1 3.5154 0.243201
## NP 2-UP 1 : A-B 1.39538 1 0.9906 0.787487
## NP 2-UP 2 : A-B 1.51547 1 1.2555 0.787487
## UP 1-UP 2 : A-B 1.08607 1 0.0535 0.817143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(glht(glmer.model4, linfct=mcp(Population='Tukey')))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmer(formula = Trans ~ Population * Repeat + (1 | VIE/Repeat),
## data = activity, family = poisson)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## NP 2 - NP 1 == 0 0.6558 0.2846 2.304 0.0969 .
## UP 1 - NP 1 == 0 0.3572 0.2713 1.317 0.5514
## UP 2 - NP 1 == 0 0.8350 0.3043 2.744 0.0308 *
## UP 1 - NP 2 == 0 -0.2986 0.2825 -1.057 0.7150
## UP 2 - NP 2 == 0 0.1792 0.3143 0.570 0.9408
## UP 2 - UP 1 == 0 0.4779 0.3022 1.581 0.3886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)summary(glht(glmer.model4, linfct=mcp(Repeat='Tukey')))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmer(formula = Trans ~ Population * Repeat + (1 | VIE/Repeat),
## data = activity, family = poisson)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## B - A == 0 0.7435 0.2265 3.282 0.00103 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)plot(glmer.model4)
glmer.model5<-glmer(Trans~Population*Repeat+(1|VIE), data=activity, family=poisson)
summary(glmer.model5)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Trans ~ Population * Repeat + (1 | VIE)
## Data: activity
##
## AIC BIC logLik deviance df.resid
## 2027.7 2052.6 -1004.8 2009.7 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.2470 -1.5157 -0.1769 1.5146 11.2239
##
## Random effects:
## Groups Name Variance Std.Dev.
## VIE (Intercept) 0.3115 0.5581
## Number of obs: 118, groups: VIE, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.80208 0.14000 27.157 < 2e-16 ***
## PopulationNP 2 0.51012 0.20694 2.465 0.01370 *
## PopulationUP 1 0.12478 0.19768 0.631 0.52788
## PopulationUP 2 0.63112 0.22148 2.850 0.00438 **
## RepeatB 0.49689 0.04300 11.556 < 2e-16 ***
## PopulationNP 2:RepeatB -0.86260 0.06283 -13.728 < 2e-16 ***
## PopulationUP 1:RepeatB -0.40229 0.06097 -6.598 4.17e-11 ***
## PopulationUP 2:RepeatB -0.45456 0.06038 -7.528 5.16e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PplNP2 PplUP1 PplUP2 RepetB PNP2:R PUP1:R
## PopulatnNP2 -0.676
## PopulatnUP1 -0.708 0.479
## PopulatnUP2 -0.632 0.428 0.448
## RepeatB -0.191 0.129 0.135 0.121
## PpltnNP2:RB 0.131 -0.155 -0.093 -0.083 -0.684
## PpltnUP1:RB 0.135 -0.091 -0.177 -0.085 -0.705 0.483
## PpltnUP2:RB 0.136 -0.092 -0.096 -0.155 -0.712 0.487 0.502testInteractions(glmer.model5)## Chisq Test:
## P-value adjustment method: holm
## Value Df Chisq Pr(>Chisq)
## NP 1-NP 2 : A-B 0.42206 1 188.460 < 2.2e-16 ***
## NP 1-UP 1 : A-B 0.66879 1 43.532 1.251e-10 ***
## NP 1-UP 2 : A-B 0.63473 1 56.667 2.581e-13 ***
## NP 2-UP 1 : A-B 1.58456 1 53.391 1.094e-12 ***
## NP 2-UP 2 : A-B 1.50387 1 42.724 1.260e-10 ***
## UP 1-UP 2 : A-B 0.94908 1 0.745 0.3881
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(glht(glmer.model5, linfct=mcp(Population='Tukey')))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmer(formula = Trans ~ Population * Repeat + (1 | VIE), data = activity,
## family = poisson)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## NP 2 - NP 1 == 0 0.5101 0.2069 2.465 0.0653 .
## UP 1 - NP 1 == 0 0.1248 0.1977 0.631 0.9218
## UP 2 - NP 1 == 0 0.6311 0.2215 2.850 0.0226 *
## UP 1 - NP 2 == 0 -0.3853 0.2066 -1.865 0.2427
## UP 2 - NP 2 == 0 0.1210 0.2295 0.527 0.9524
## UP 2 - UP 1 == 0 0.5063 0.2212 2.289 0.1001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)anova(glmer.model4, glmer.model5)## Data: activity
## Models:
## glmer.model5: Trans ~ Population * Repeat + (1 | VIE)
## glmer.model4: Trans ~ Population * Repeat + (1 | VIE/Repeat)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## glmer.model5 9 2027.7 2052.6 -1004.8 2009.7
## glmer.model4 10 1252.0 1279.7 -616.0 1232.0 777.67 1 < 2.2e-16
##
## glmer.model5
## glmer.model4 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Again, glmer.model4 clearly the stronger one.
tapply(Trans, list(Population, Repeat), mean)## A B
## NP 1 51.00000 83.82353
## NP 2 82.78571 57.42857
## UP 1 59.88235 65.82353
## UP 2 98.81818 103.09091tapply(Trans, list(Population, Repeat), var)## A B
## NP 1 1611.750 1734.279
## NP 2 2318.489 1685.648
## UP 1 1434.985 1934.404
## UP 2 3424.764 2842.091#this provides us with means and variances, as you can see the NP's differ alot between trialsactivity.lme3<-lme(Trans~Population*Repeat, random=~1|VIE/Repeat, correlation=corAR1(form=~1|VIE/Repeat), method='REML', data=activity)
newdata<-expand.grid(Population=levels(activity$Population), Repeat=levels(activity$Repeat))
coefs<-fixef(activity.lme3)
Xmat<-model.matrix(~Population*Repeat, data=newdata)
fit<-as.vector(Xmat %*% coefs)
newdata$fit<-fit
head(newdata)## Population Repeat fit
## 1 NP 1 A 51.00000
## 2 NP 2 A 82.78571
## 3 UP 1 A 59.88235
## 4 UP 2 A 98.81818
## 5 NP 1 B 83.82353
## 6 NP 2 B 57.42857se<-sqrt(diag(Xmat %*% vcov(activity.lme3) %*% t(Xmat)))
newdata$lower<-fit-2*se
newdata$upper<-fit+2*se
head(newdata)## Population Repeat fit lower upper
## 1 NP 1 A 51.00000 29.20028 72.79972
## 2 NP 2 A 82.78571 58.76359 106.80784
## 3 UP 1 A 59.88235 38.08263 81.68208
## 4 UP 2 A 98.81818 71.71757 125.91879
## 5 NP 1 B 83.82353 62.02381 105.62325
## 6 NP 2 B 57.42857 33.40645 81.45070newdata$Population<-factor(newdata$Population, levels=c("NP 1", "NP 2", "UP 1", "UP 2"), labels=c("Lane Cove (NP 1)", "Kurin-gai(NP 2)", "Newtown (UP 1)", "Padstow (UP 2)"))
head(newdata)## Population Repeat fit lower upper
## 1 Lane Cove (NP 1) A 51.00000 29.20028 72.79972
## 2 Kurin-gai(NP 2) A 82.78571 58.76359 106.80784
## 3 Newtown (UP 1) A 59.88235 38.08263 81.68208
## 4 Padstow (UP 2) A 98.81818 71.71757 125.91879
## 5 Lane Cove (NP 1) B 83.82353 62.02381 105.62325
## 6 Kurin-gai(NP 2) B 57.42857 33.40645 81.45070nondirectactivity<-ggplot(newdata, aes(y=fit, x=Population, fill=Repeat))+
geom_bar(position=position_dodge(.9), stat="identity", color="black")+
geom_errorbar(aes(ymin = lower, ymax = upper), width=.2, position=position_dodge(0.9))+
coord_cartesian(ylim=c(0,150))+
scale_fill_manual(values=c("#404040", "#C0C0C0"))+
theme_classic()+
scale_y_continuous("Number of Transitions")+
scale_x_discrete("")+
theme(text = element_text(size = 10),axis.title.y = element_text(vjust=2, size=rel(2)), axis.title.x = element_text(vjust=-1, size=rel(2)), axis.text= element_text(size = 14), legend.position=("none"))
nondirectactivity