#binomial likelihood function for magnolia presence-absence #deterministic component includes intercept and linear and quadratic coefficients for elevation #scaling each variable by max cover = mf\$cover elev = mf\$elev/1000 #elevation in easier units (km) streamdist = mf\$streamdist/1000 #streamdist in easier units (km) #likelihood function fn = function(par) {-sum(dbinom(as.numeric(cover>0),prob=par[1]+par[2]*elev+par[3]*(elev^2),size=1,log=T))} #optimization out = optim(fn,par=c(.1,0.0001,0)) #starting values for binomial likelihoods can be tricky due to 0-1 constraint plot(elev,cover>0) x = seq(0,2,length=100) lines(x,out\$par[1]+out\$par[2]*x+out\$par[3]*(x^2),col="purple",lwd=2) #likelihood profile of the quadratic term coefs = seq(-.35,-.1,length=40) #vector of possible slope values surrounding the MLE lik.out = numeric(length(coefs)) #holding vector (currently zeros) fn2 = function(par,m) { #new simplified likelihood function (only 2 fitted pars) -sum(dbinom(as.numeric(cover>0),prob=par[1]+par[2]*elev+m*(elev^2),size=1,log=T))} for(i in 1:length(coefs)) { #calculates new likelihoods for each new slope model lik.out[i] = optim(fn2,par=c(0,.7),m=coefs[i])\$value} plot(coefs,lik.out,type="l"); abline(v=out\$par[3],col="blue",lwd=2) #profile with MLE #examining the right side (one-tailed) 95% CI to assess overlap with zero upper.s = coefs[which.min(lik.out):length(lik.out)] #x values, right side upper.k = lik.out[which.min(lik.out):length(lik.out)] #y values, right side upperCI = approx(upper.k,upper.s,xout=min(lik.out)+qchisq(0.95,1)/2)\$y abline(v=upperCI,lty=2,lwd=2,col="red")