线性混合模型系列五：REML实战

# zhangqin laoshi PPT codeherd <- factor(c(1,2,2,1,1,2,1,2,2))
sire <- factor(c(1,1,1,2,2,3,4,4,4))
y <- c(240,190,170,180,200,140,170,100,130)
REML_data <- data.frame(herd,sire,y)
X <- model.matrix(y~herd-1)
Z <- model.matrix(y~sire-1)
A <- matrix(c(1,0.25,0,0,0.25,1,0,0,0,0,1,0,0,0,0,1),nr=4)

mixmodel.loglik<-function(theta,y){
  sigmag2<-theta[1]
  sigmae2<-theta[2]

  n = length(y)
  G = A*sigmag2
  R = diag(n)*sigmae2
  V = Z %*% G %*% t(Z) + R
  Vi = solve(V)

  P = Vi - Vi %*% X %*% solve(t(X) %*% Vi %*% X) %*% t(X) %*% Vi
  logl = -0.5*(log(det(V)) + log(det(t(X) %*% Vi %*% X)) + t(y) %*% P %*% y)  return(-logl)
}

optim(c(1,1),mixmodel.loglik,y=y)

mixmodel.loglik2 <-function(theta,y){
  sigmag2<-theta[1]
  sigmae2<-theta[2]

  n = length(y)
  G = A*sigmag2
  R = diag(n)*sigmae2
  V = Z %*% G %*% t(Z) + R
  Vi = solve(V)
  P = Vi - Vi %*% X %*% solve(t(X) %*% Vi %*% X) %*% t(X) %*% Vi
  C11 = t(X) %*%  solve(R) %*% X
  C12 = t(X) %*%  solve(R) %*% Z
  C21 = t(Z) %*%  solve(R) %*% X
  C22 = t(Z) %*%  solve(R) %*% Z + solve(G)
  C = rbind(cbind(C11,C12),cbind(C21,C22))
  logl2 = -0.5*(log(det(C)) + log(det(R)) + log(det(G)) + t(y) %*% P %*% y)  return(-logl2)
}
optim(c(1,1),mixmodel.loglik2,y=y)

herd <- factor(c(1,2,2,1,1,2,1,2,2))
sire <- factor(c(1,1,1,2,2,3,4,4,4))
y <- c(240,190,170,180,200,140,170,100,130)
REML_data <- data.frame(herd,sire,y)
X <- model.matrix(y~herd-1)
Z <- model.matrix(y~sire-1)
XX <- crossprod(X)
XZ <- crossprod(X,Z)
ZX <- t(XZ)
ZZ <- crossprod(Z)
LHS <- rbind(cbind(XX,XZ),cbind(ZX,ZZ))
Xy <- crossprod(X,y)
Zy <- crossprod(Z,y)
RHS <- rbind(Xy,Zy)
A <- matrix(c(1,0.25,0,0,0.25,1,0,0,0,0,1,0,0,0,0,1),nr=4)
Ainv <- solve(A)
yy <- crossprod(y)
N <- length(y)
rankX <-qr(X)$rank
q <- length(unique(sire))
nh <- length(unique(herd))
k0 <- 12threshold <- 0.00000001write.table(t(c("sigmaE","sigmaS","k")),file="results.txt",row.names = FALSE, col.names = FALSE)repeat{
    LHS1 <- LHS
    k <- k0
    LHS1[(nh+1):(nh+q), (nh+1):(nh+q)] <- LHS1[(nh+1):(nh+q), (nh+1):(nh+q)]+Ainv*k
    sol <- solve(LHS1,RHS)
  C <- solve(LHS1)
  Czz <- C[(nh+1):(nh+q), (nh+1):(nh+q)]
    sigmaE <- (yy - crossprod(sol,RHS))/(N-rankX)
    sAs <- t(sol[(nh+1):(nh+q)]) %*% Ainv %*% sol[(nh+1):(nh+q)]
    trCA <- sum(diag(Czz %*% Ainv))
    sigmaS <- (sAs + trCA*sigmaE)/q
    k0 <- as.numeric(sigmaE/sigmaS)
    out <- c(sigmaE,sigmaS,k0)
    write.table(t(out), file="results.txt", append = TRUE, row.names = FALSE, col.names = FALSE)
    if(abs(k-k0) < threshold) break}
results <- read.table("results.txt",header = TRUE)
results

# 自然对数的3次方exp(3)

# 对上面结果求自然对数log(exp(3))

log(10)

log10(10)