# fits logistic and probit models to Meron data m <- read.table(file="meron.200.glim.txt", header=TRUE ) m[1:5,] attach(m) N = length(tree.92) N # define appropriate indicator variables goats = 1 * (graze.type == 2) # so will be a numeric variable cattle.low = (graze.type == 1)*(graze.intensity ==1) cattle.med.high = (graze.type == 1)*(graze.intensity >= 2) # compute a linear regression linreg.1 = lm(tree.92 ~ alt + stream.dist + road.dist + water.dist + slope + ns.aspect + ew.aspect + ns.aspect*slope + ew.aspect*slope + cattle.low + cattle.med.high + goats ) summary(linreg.1) # fit logistic regression. logit.1 = glm(tree.92 ~ alt + stream.dist + road.dist + water.dist + slope + ns.aspect + ew.aspect + ns.aspect*slope + ew.aspect*slope + cattle.low + cattle.med.high + goats, family = binomial(link=logit)) summary(logit.1) # fit probit regression probit.1 = glm(tree.92 ~ alt + stream.dist + road.dist + water.dist + slope + ns.aspect + ew.aspect + ns.aspect*slope + ew.aspect*slope + cattle.low + cattle.med.high + goats, family = binomial(link=probit)) summary(probit.1) # fit probit regression probit.1 = glm(tree.92 ~ alt + stream.dist + road.dist + water.dist + slope + ns.aspect + ew.aspect + ns.aspect*slope + ew.aspect*slope, family = binomial(link=probit)) summary(probit.1) # compare probit and logistic link functions p = seq(0.005, 0.995, by=0.005) logit = log(p / (1 - p)) probit = qnorm(p) plot(logit, p, type="n", xlab="logit(p) (---) and probit(p) (solid)") lines(probit, p, lty=1) lines(logit, p, lty=2)