Didrik Vanhoenacker, 2009
Mission 7. Meta population. Generalized Linear Models (4 November)
Learning goals
- Be able to use generalized linear models, glm(), to test binary response data.
- Be able to calculate some isolation measures that are very important in ecology.
- Get insights into meta population ecology.
- Be better at making logistic regression graphs.
- Be better at interpreting test tables and model comparisons.
Data collection
- Use the baggdata or vattenvaexter data set on
- Choose one species.
- Choose one site size measurement.
- Calculate one isolation measure.
Mission
- Calculate some isolation measure.
- Test if the presence of a species depends on site size and or isolation.
- Make one or two logistic regression graphs to show the relationship(s).
Hand in
- One Mini-report on one species.
Mail it to me, , and to your feed-back pal.
Feedback pal
Return the favour of your last feed back pal:
Niklas Shayer Esther Wang Joana Yao Nizar Karin J Charles Peter Yong Nok Li Rushana Malin Tahir Niklas
Meta population
These are data sets from PhD students here at Botan. Hannah Östergård has contributed her data set (”baggdata.xls”) on seed predators of Lathyrus vernus (=Vårärt). There is one weevil beetle and one bruchid beetle. As measures of habitat size she has recorded both patch area and the number of flowering plants in a patch.
Johan Dahlgrens has contributed his data set (”vattenvaexter.xls”) on lake plants, Potamogeton natans (gäddnate) and Isoëtes lacustris (braxengräs). As measures of habitat size he has recorded lake size and littoral size.
- Choose one data set and one species!
- Choose one measure of habitat size.
- Calculate one measure of isolation, either connectivity or closest neighbour (the measure you think is most relevant or the best predictor)
- Analyse the effect of population size and isolation.
Distance matrix
With the function dist you can calculate the distances between all combinations of sites. Divide the coordinates with 1000 to get them in km:s instead of meters. Km is the standard in connectivity analyses, while the coordinates in the Swedish grid is in meters.
distance<-dist(cbind(east/1000,north/1000))
Calculate connectivity
This rather complicated formula calculates the connectivity. Use copy-paste! I have copy-pasted it myself from Jari Oksanen. Alpha is a parameter for how rapidly the importance of a neighbour patch decreases with distance. The first line that starts connectmatrix is the actual formula for how the importance of a neighbour patch decreases with distance. The second line that starts with connectmatrix multiplies it with the site sizes. Then the importance of other populations are summed for each site.
distance<-dist(cbind(east/1000,north/1000))
alpha<-1
connectmatrix <-as.matrix(exp(-alpha*distance))
connectmatrix <-sweep(connectmatrix,2,area,"*")
connectivity<-rowSums(connectmatrix)
connectivity
Calculate nearest neighbour
To calculate the nearest neighbour population you can use the same connectivity matrix as the one used for connectivity. The for loop looks for the smallest value in each row, that is the shortest distance to any other population. The for loop does this for all rows (= populations). Because the distance from pop1 to pop1 (to itself!) of course is 0, you need the row diag(distmatrix)<-max(distmatrix) which substitutes all the zeroes with the max distance to another population. Otherwise all rows (= populations) would get 0 meters to nearest neighbour.
distance<-dist(cbind(east/1000,north/1000))
distmatrix<-as.matrix(distance)
diag(distmatrix)<-max(distmatrix)
nearest<-NULL
for (i in 1:length(east)){
nearest[i]<-min(distmatrix[,i])
}
nearest
Analyse the effect of pop size and isolation
Now you can build models with patch size and isolation (the measure you chose and calculated, respectively). Areas and population sizes are usually logged, but this is not done for connectivity or nearest neighbour (and area and pop size is not logged when calculating connectivity).
mod.int<-glm(brushid~log.area+connectivy+log.area:connectivy, family=binomial)
mod.both<-glm(brushid~log.area+connectivy, family=binomial)
anova(mod.int, mod.both, test="Chi")
If you don't want to do model comparison, make an anova table with the package car.
library(car)
Anova(glm(brushid~log.area*connectivy, family=binomial))
Remember that the * sign means both the main effects AND their interaction (y~x1*x2 = y~x1+x2+x1:x2).
Reduce the models by remove unnecessary explanatory variables.
Model assumptions
As this is generalized linear model with binomial errors you should of course not check for normally distributed residuals or non-constant variance.
Graphs
To illustrate the relationships graphically you can do logistic regression graphs [ Pimp your graph ]. As both these explanatory variables are continuous it is easiest to do two separate graphs, one for each explanatory. These graphs do not correspond exactly to your test (y~x1 and y~x2 rather than y~x1 given x2 and y~x2 given x1). But it is very hard to make sensible 3D-graph with a categorical response.
Map graphs (not necessary, just for fun)
You can also make pretty neat illustrations of the actual spatial pattern in R. Just use the coordinate variables east and north. Use asp to get the same scale in the x- and y-axis, or your map might become stretched.
plot(east,north,asp=1)
Or with some pimping:
plot(east,north,asp=1,cex=2*no.flowering/max(no.flowering)+2,pch=21, lwd=2,bg=ifelse(brushid=="YES","grey30","white"))
Probably a legend would be a good idea:
legend("bottomright", title="Brushid", pch=21, pt.cex=2 ,pt.lwd=2, legend=c("Present","Absent"),pt.bg=c("grey30","white"))
A really neat variant is to make a graph without axes and text and then paste it onto a map that you have downloaded from the digital map library, google maps, or eniro.
plot(east,north,asp=1,cex=2* no.flowering/max(no.flowering)+2,pch=21, lwd=2,bg=ifelse(brushid=="YES","grey30","white"), axes=F, xlab="",ylab="")
My primrose populations on Öland…
1