Tf3 = Tf2[,c("minT","Transect","logger","ELEV","LOG.STRDST","minSYN")] #our variables of interest from lme5 Tf3 = na.exclude(Tf3) #remove rows containing NAs (equivalent to na.method=na.exclude in lme and lmer) attach(Tf3) trans = as.numeric(Transect) #treat transect levels numerically logger2 = as.numeric(logger) #treat logger levels numerically T.effect = rep(0,dim(Tf3)[1]) #create empty vector for random transect effect l.effect = rep(0,dim(Tf3)[1]) #create empty vector for random logger effect ML1 = function(par) { #likelihood function for optim for(a in 1:length(unique(logger2))) { #generate a 'random effect' for each logger l.effect[logger2==a] = rnorm(length(logger2[logger2==a]),mean=0,sd=log(par[5])) } for(b in 1:length(unique(trans))) { #generate a 'random effect' for each transect T.effect[trans==b] = rnorm(length(trans[trans==b]),mean=0,sd=log(par[6])) } y = par[1] + par[2]*ELEV + par[3]*LOG.STRDST + par[4]*minSYN + T.effect + l.effect -sum(dnorm(minT,mean=y,sd=log(par[7]),log=T)) #variation within loggers (residual) } start = c(1,0,0,0,3,3,3) #initial parameter values for optim out1 = optim(ML1,par=start,method="BFGS") #I got BFGS to converge, but not SANN or Nelder-Mead #variance components log(out1\$par[5:7])^2 / sum(log(out1\$par[5:7])^2) ML2 = function(par) { #likelihood function for optim for(b in 1:length(unique(trans))) { #generate a 'random effect' for each transect T.effect[trans==b] = rnorm(length(trans[trans==b]),mean=0,sd=log(par[5])) } y = par[1] + par[2]*ELEV + par[3]*LOG.STRDST + par[4]*minSYN + T.effect -sum(dnorm(minT,mean=y,sd=log(par[6]),log=T)) #variation within loggers (residual) }