代码之家 › 专栏 › 技术社区 › Benjamin Christoffersen

`nlme`具有交叉随机效果

random-effects nlme r

6

Benjamin Christoffersen · 技术社区 · 6 年前

我在试着装一个十字架 非线性 随机效应模型 线性随机效应模型 如前所述 in this question 在这里面 mailing list post 使用 nlme 包裹。尽管如此,无论我尝试什么,我都会出错。下面是一个例子

library(nlme)

#####
# simulate data
set.seed(18112003)
na <- 30
nb <- 30

sigma_a <- 1
sigma_b <- .5
sigma_res <- .33

n <- na*nb

a <- gl(na,1,n)
b <- gl(nb,na,n)
u <- gl(1,1,n)

x <- runif(n, -3, 3)

y_no_noise <- x + sin(2 * x)
y <- 
  x + sin(2 * x) + 
  rnorm(na, sd = sigma_a)[as.integer(a)] + 
  rnorm(nb, sd = sigma_b)[as.integer(b)] + 
  rnorm(n, sd = sigma_res)

#####
# works in the linear model where we know the true parameter
fit <- lme(
  # somehow we found the right values
  y ~ x + sin(2 * x), 
  random = list(u = pdBlocked(list(pdIdent(~ a - 1), pdIdent(~ b - 1)))))
vv <- VarCorr(fit)
vv2 <- vv[c("a1", "b1"), ]
storage.mode(vv2) <- "numeric"
print(vv2,digits=4)
#R    Variance StdDev
#R a1    1.016 1.0082
#R b1    0.221 0.4701

#####
# now try to do the same with `nlme`
fit <- nlme(
  y ~ c0 + sin(c1),
  fixed = list(c0 ~ x, c1 ~ x - 1),
  random = list(u = pdBlocked(list(pdIdent(~ a - 1), pdIdent(~ b - 1)))), 
  start = c(0, 0.5, 1))
#R Error in nlme.formula(y ~ a * x + sin(b * x), fixed = list(a ~ 1, b ~  : 
#R    'random' must be a formula or list of formulae

这个 lme 这个例子类似于“S和S-PLUS混合效应模型”的163-166页,只有2个随机效应而不是3个。

1 回复 | 直到 6 年前

1

4

Benjamin Christoffersen 6 年前

我应该用双面公式 help("nlme")

fit <- nlme(
  y ~ c0 + c1 + sin(c2),
  fixed = list(c0 ~ 1, c1 ~ x - 1, c2 ~ x - 1),
  random = list(u = pdBlocked(list(pdIdent(c0 ~ a - 1), pdIdent(c1 ~ b - 1)))), 
  start = c(0, 0.5, 1))

# fixed effects estimates
fixef(fit)
#R c0.(Intercept)           c1.x           c2.x 
#R     -0.1788218      0.9956076      2.0022338

# covariance estimates
vv <- VarCorr(fit)
vv2 <- vv[c("c0.a1", "c1.b1"), ]
storage.mode(vv2) <- "numeric"
print(vv2,digits=4)
#R       Variance StdDev
#R c0.a1   0.9884 0.9942
#R c1.b1   0.2197 0.4688