In rpkgs/JOPSbook: P-Spline

线性混合模型 {#Mixed-Model}

案例1

两组CO2实验，高CO2浓度和低CO2浓度。每组3棵树（tree）、每棵树测4片叶子的气孔面积（area）。

问题：CO2对气孔导度的影响

错误的方法，固定效应模型

library(JOPSbook)
library(gamair)
library(magrittr)
library(data.table)

data("stomata")
stomata %<>% data.table()

summary(stomata)
print(stomata)

$$ y_i = CO2 * α_j + tree * β_k + ǫ_i $$

其中i代表第i个观测数据，j代表CO2浓度，k代表树。

m1 <- lm(area ~ CO2 + tree, stomata)
m0 <- lm(area ~ CO2, stomata)
anova(m0, m1)

通过anova方差分析可以看到，tree对气孔的影响显著。

m2 <- lm(area ~ tree, stomata)
anova(m2, m1)

m2和m1的方差分析显示，CO2对气孔的影响也显著。CO2、tree对气孔是协同影响。

st = stomata[, .(area = mean(area)), .(tree, CO2)]
st

m3 <- lm(area ~ CO2, st)
summary(m3)
# anova(m3)

summary(m3)$sigma^2 - summary(m1)$sigma^2 / 4

案例2

library(nlme)

data(Machines)
Machines %<>% data.table()

summary(Machines)
head(Machines)

attach(Machines)
interaction.plot(Machine, Worker, score)

llm <- function(theta, X, Z, y) {
  ## untransform parameters...
  sigma.b <- exp(theta[1])
  sigma <- exp(theta[2])
  ## extract dimensions...
  n <- length(y)
  pf <- ncol(X)
  pr <- ncol(Z)
  ## obtain \hat \beta, \hat b...
  X1 <- cbind(X, Z)
  ipsi <- c(rep(0, pf), rep(1 / sigma.b^2, pr))
  b1 <- solve(
    crossprod(X1) / sigma^2 + diag(ipsi),
    t(X1) %*% y / sigma^2
  )
  ## compute log|Z’Z/sigma^2 + I/sigma.b^2|...
  ldet <- sum(log(diag(chol(
    crossprod(Z) / sigma^2 + diag(ipsi[-(1:pf)])
  ))))
  ## compute log profile likelihood...
  l <- (-sum((y - X1 %*% b1)^2) / sigma^2 - sum(b1^2 * ipsi) -
    n * log(sigma^2) - 
    pr * log(sigma.b^2) - 
    2 * ldet - 
    n * log(2 * pi)) / 2
  attr(l, "b") <- as.numeric(b1) ## return \hat beta and \hat b
  -l
}

# 将常数项去掉n
ormal.lik2<-function(theta, y){
  mu <- theta[1]
  sigma2 <- theta[2]
  n = length(y)
  z = (y-mu)/sigma  # logl = -0.5*n*log(sigma^2) - (1/(2*sigma^2)*sum((y-mu)^2)) 
  logl = 
    -1 / 2 * n * log(2 * pi) - 
     1 / 2 * n * log(sigma2) - 
     (1 / (2 * sigma2)) * sum((y - mu)**2)
  # `-1 / 2 * n * log(2 * pi)`: 常数项
  # logl = - 1 / 2 * n * log(sigma2) - (1 / (2 * sigma2)) * sum((y - mu)**2)
  -logl
}

# https://cloud.tencent.com/developer/article/1542008