R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model分析藻类数据实例（上）：https://developer.aliyun.com/article/1496401
右侧渐近线中的方差估计值是非零的。
加入随机效应后，参数根本就没有什么变化。
最大的比例差异是3.1%（在比例参数中）。
nlmefit2 <- update(list(asyR+xmd+scal+asp ~1),
start )
我们可以通过AIC或似然比检验来比较模型
AICtab(nlmefit1,nlmefit2,weights=TRUE)
anova(nlmefit1,nlmefit2)
可以做一个F测试而不是 LRT（即考虑到有限大小的修正）。
pchisq(iff,df=2,lower.tail=FALSE)
##分母非常大的F检验。
pf(diff/2,df1=2,df2=1000000,lower.tail=FALSE)
我们不知道真正相关的df，但上面的总结表明df是40。
nlmer
我想现在可以为nlmer得到正确的模型规范，但我找不到一个方便的语法来进行固定效应建模（即在这种情况下允许一些参数因组而异）--当我构建了正确的语法，nlmer无法得到答案。
基本的RE模型（没有群体效应）运行良好。
nlmer(
X ~ SSfpl(Day, asy, as, x, s) ~
asy|Indi,)
根据我的理解，人们只需要构建自己的函数来封装固定效应结构；为了与nlmer一起使用，该函数还需要计算相对于固定效应参数的梯度。这有点麻烦，但可以通过修改派生函数生成的函数，使之稍微自动化。
构建虚拟变量:
mm <- model.matrix(~Group,data=d)
grp2 <- mm\[,2\]
构建一个函数来评估预测值及其梯度；分组结构是硬编码的。
deriv(~A+((B0+B1\*grp2+B2\*grp3-A)/(1+exp((x-xmid)/scale)
通过插入与传递给函数的参数名称相匹配的行来查看所产生的函数，并将这些参数名称分配给梯度矩阵。
L1 <- grep("^ +\\\.value +<-")
L2 <- grep("^ +attr\\\(\\\.value",)
eval(parse(text))
尝试一下拟合:
nlmer( X ~ fpl(Day, asym, as, asymp, asR3, xmi, sca) ~ as|Indi, start = list(nlpars)),data=d)
失败了（但我认为这是由于nlmer本身造成的，而不是设置有什么根本性的问题）。为了确定，我应该按照同样的思路生成一个更大的人工数据集，看看我是否能让它工作起来。
现在我们可以用稳定版（lme4.0）得到一个答案。
结果不理想
fixef(nlmerfit2)
range(predict(nlmerfit2))
我不能确定，在nlmer中是否有更简单的方法来做固定效果。
AD模型生成器
我们还可以使用AD模型生成器来解决这个问题。它可以处理更复杂的模型，比如拟合更多参数的群体效应。
部分原因是我对ADMB的熟悉程度较低，这有点费劲，最后我通过循序渐进的步骤才成功。
最小的例子
首先尝试没有随机效应、分组变量等。(即等同于上面的nls拟合）。)
##设置数据：调整名称，等等 d0 <- c(list(nobs=nrow(d)),as.list(d0)) ##起始值：调整名称，增加数值 names(svec3) <- gsub("\\\.","",names(svec3)) ## 移除点 svec3$asympR <- 0.6 ## 单一值 ## 运行 do_admb("algae0", data, params, run.opts)
结果不错
固定效应模型
现在尝试用固定效应分组，使用上面构建的虚拟变量（也可以使用if语句，或者用R[Group[i]]的for循环中的R值向量，或者（最佳选择）为R传递一个模型矩阵...）。我们必须使用elem_div而不是/来对两个向量进行元素除法。
model1 <- "
参数部分
向量 pred(1,nobs) // 预测值
向量Rval(1,nobs) //预测值
过程部分
pred = as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal)
"
试着用模型矩阵来代替它。
model1B <- "
参数部分
向量 pred(1,nobs) // 预测值
向量Rval(1,nobs) //预测值
过程部分
pred = asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc)) 。
"
当然，在参数相同的情况下，也可以工作。
随机效应
现在添加随机效应。回归函数并没有完全实现随机效应模型（尽管这应该在即将到来的版本中被修复），所以我们用公式减去（n/2 log({RSS}/n)），其中RSS是残差平方和。
model2 <- "
参数部分
向量 pred(1,nobs) // 预测值
向量Rval(1,nobs) //预测值
过程部分
pred = asym+elem
f = 0.5\*no\*log(norm2(X-pr)/n)+norm2(R)。
"
由于ADMB不处理稀疏矩阵，也不惩罚循环，如果将随机效应实现为(i=1; i<=nobs; i++) Rval[i] += Rsigma*Ru[Group[i]]，效率会略高，但我是懒人/我喜欢矩阵表示的紧凑性和可扩展性.
现在我们终于可以测试R以外的参数的固定效应差异了。
model3 <- "
参数部分
向量 prd(1,nobs) // 预测值
向量Rl(1,nobs) // 预测值
向量 scalal(1,nobs)
向量xmal(1,nobs)
sdror opr(1,nobs) //输出预测值
程序部分
Rval = XR\*Rve+Rsma\*（Z*Ru）。
xmval = Xd*xdvec;....
f = 0.5\*nobs\*log(norm2(X-pred)/nobs)+norm2(Ru)
"
结果:
summary(admbfit3)
有一个非常大的AIC差异。如上文所示，对nlme拟合的似然比F测试是作为一种练习......
对于该图，最好是按组指定参数重新进行拟合，而不是按基线+对比度进行拟合。
fit3B <- do_admb(, data, params, re, run.opts=run.control) plot2(list(cc),intercept=TRUE)
现在我们对标准化的问题很困扰，所以（经过一番折腾）我们可以在不同的面板上重新画出群体变化的参数。
诊断图
##放弃条件模式/样本-R估计值
diagplot1 %+% dp2
也许这暗示了两个实验组中更大的差异？
拟合与残差
diagplot2 %+% dp2
叠加预测（虚线）:
g1 + geom_line
如果能生成平滑的预测曲线（即对中间的日值），那就更好了，但也更繁琐。
结论
- 从参数估计中得出的主要结论是，第三组下降得更早一些（xmidvec更小），同时下降得更远（Rvec更低）。
似然分析
计算一个( sigma^2_R ) 似然函数的代码并不难，但运行起来有点麻烦：它很慢，而且计算在置信度下限附近的几个点上出现了非正-无限矩阵；我运行了另一组值，试图充分覆盖这个区域。
lapply(Rsigmavec,fitfun)
## 尝试填补漏洞
lapply(Rsigmavec2,fitfun)
带有插值样条的剖面图和似然比检验分界线。
在sigma^2_R 上的95%剖面置信区间是{0.0386,0.2169}。
我没有计算过，但转换后的剖面图（在对应于偏离度与最小偏离度的平方根偏差的 y ）上，所以二次剖面将是一个对称的V）显示，二次近似对这种情况相当糟糕 ...
ggplot(sigma,sqrt(2*(NLL-min(NLL))+
geom_point()
扩展
- 更多地讨论分母df问题。参数引导法/MCMC？
- 我们可以尝试在xmid和scale参数中加入随机效应。
- 在组间或作为X的函数的方差（无论是残差还是个体间的方差）中可能有额外的模式。
