library(nlme)
# 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
reml<-function(beta){
y<-Rail$travel
X<-rep(1,n*r)
Z<-kronecker(diag(n),rep(1,r))
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)
}
# minimize reml function with optim()
optim(c(66.5, 5, 2),reml)
$par
[1] 66.505645 6.421954 2.782963
> sqrt(exp(c(6.421954,2.782963)))
[1] 24.803307 4.020802
# ML: 22.62435 4.020779
# REML: 24.80547 4.020779
lme(travel~1, random=~1|Rail, data=Rail, method="ML")
lme(travel~1, random=~1|Rail, data=Rail)
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment