Hi Juliet,

I had the same idea as you and created the following code for linear mixed
models.

Have fun~

cahn

### This fits Rail with nlm() (linear mixed models).

data(Rail, package="nlme")
lme(travel~1, random=~1|Rail, data=Rail)

*fixed: 66.5
random: 24.80547 4.020779 (standard deviations)*

### The followings set up REML and maximize it with optim()
### the last two elements of optim output needs to be sqrt(log())

# 6 Rails have 3 repeatitions each.
n<-6; r<-3;

# -log Restricted Maximum Likelihood
# beta[1]: mean

# exp(beta[2]): sigma_a

# exp(beta[3]): sigma
# exp is used to prevent negative estimator for variance

### Design Matrix

y<-Rail$travel
X<-rep(1,n*r)
Z<-kronecker(diag(n),rep(1,r))

### Restricted Maximum Likelihood
reml<-function(beta){
  V<-exp(beta[2])*Z%*%t(Z)+exp(beta[3])*diag(n*r)

ml<-log(det(V))/2+t(y-beta[1])%*%solve(V)%*%(y-beta[1])/2+log(det(t(X)%*%solve(V)%*%X))/2
  return(ml)
}

#### optim()

optim(c(66.5,5,2), reml, control=list(reltol=1e-20))$par





On Thu, Jan 14, 2010 at 2:48 PM, Juliet Hannah <juliet.hannah@gmail.com>wrote:

> Hi List,
>
> I'm trying to learn more about mixed models in detail. Does anyone know of
> any examples where code is provided showing simple examples of doing
> this from scratch. For example, setting up a likelihood, using optim, and
> things
> like that.
>
> Thanks,
>
> Juliet
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 
안재형
4033 Remington Oaks Circle
Cary, NC 27519
http://cahn88.blogspot.com/

	[[alternative HTML version deleted]]


