0

nlmeを使用してモデルを適合させましたが、残差が十分に見えないと思うので、次の方法で変換を試みています

fit.nlme6 <- update(fit.nlme5, weights = varPower())

私は得るError: Singularity in backsolve at level 0, block 1。とにかくあまり意味がないと思う他のクラスを試し、さまざまなforms = ~とを試しfixedました。すべて同じエラー メッセージが表示されますが、おまけの警告メッセージが表示されることもあります。

遺留分はこちら。私は完璧に動作するはずだと思うvarPowerのに、なぜうまくいかないのですか?

残差:

詳細は、 (最大バイオマス)、(成長終了の瞬間)、および(最大成長の瞬間) のfit.nlme53 つのパラメーターを持つベータ成長関数に適合するモデルです。モデルはこんな感じw.maxt.et.m

fit.nlme5 <- update(fit.nlme4, fixed = list(t.e ~ trt + ground + trt:ground, w.max + t.m ~ 1),
                start = c(cfsTD[1:4], rep(0,2) ,cfsTD[5:6]))

2 つの場所 (グラウンド) に 3 つの処理 (trt) があります。残差が修正されたら、いくつかの対比を実行して、さまざまな場所での処理を比較します。

ここにデータがありますhttps://dl.dropbox.com/u/21080842/cobsgddv8.txt

そして、最終モデルに到達するために必要なコードは次のとおりです。

  #You'll need these functions
## Implementing the beta growth function from (Yin et al 2003)

bgfInit <- function(mCall, LHS, data){

  xy <- sortedXyData(mCall[["time"]], LHS, data)
  if(nrow(xy) < 4){
    stop("Too few distinct input values to fit a bgf")
  }
  w.max <- max(xy[,"y"])
  t.e <- NLSstClosestX(xy, w.max)
  t.m <- NLSstClosestX(xy, w.max/2)
  value <- c(w.max, t.e, t.m)
  names(value) <- mCall[c("w.max","t.e","t.m")]
  value

}


bgf <- function(time, w.max, t.e, t.m){

  .expr1 <- t.e / (t.e - t.m)
  .expr2 <- (time/t.e)^.expr1
  .expr3 <- (1 + (t.e - time)/(t.e - t.m))
  .value <- w.max * .expr3 * .expr2

  ## Derivative with respect to t.e
  .exp1 <- ((time/t.e)^(t.e/(t.e - t.m))) * ((t.e-time)/(t.e-t.m) + 1)
  .exp2 <- (log(time/t.e)*((1/(t.e-t.m) - (t.e/(t.e-t.m)^2) - (1/(t.e - t.m)))))*w.max
  .exp3 <- (time/t.e)^(t.e/(t.e-t.m))
  .exp4 <- w.max * ((1/(t.e-t.m)) - ((t.e - time)/(t.e-t.m)^2))
  .exp5 <- .exp1 * .exp2 + .exp3 * .exp4 

  ## Derivative with respect to t.m
  .ex1 <- t.e * (time/t.e)^((t.e/(t.e - t.m))) * log(time/t.e) * ((t.e - time)/(t.e -      
t.m) + 1) * w.max
  .ex2 <- (t.e - time) * w.max * (time/t.e)^(t.e/(t.e-t.m))
  .ex3 <- (t.e - t.m)^2
  .ex4 <- .ex1 / .ex3 + .ex2 / .ex3

  .actualArgs <- as.list(match.call()[c("w.max", "t.e", "t.m")])

##  Gradient
  if (all(unlist(lapply(.actualArgs, is.name)))) {
    .grad <- array(0, c(length(.value), 3L), list(NULL, c("w.max", 
                                                      "t.e", "t.m")))
    .grad[, "w.max"] <- .expr3 * .expr2
    .grad[, "t.e"] <- .exp5
    .grad[, "t.m"] <- .ex4 
    dimnames(.grad) <- list(NULL, .actualArgs)
    attr(.value, "gradient") <- .grad
  }
    .value
}

SSbgf <- selfStart(bgf, initial = bgfInit, c("w.max", "t.e", "t.m"))

#Now for the data and fitting
grow<-read.table("cobsgddv8.txt", header=T)

library(nlme)

grow10<-subset(grow, grow$year == "2010")
grow10$EU<- with(grow10, factor(ground):factor(plot))
grow10G<-groupedData(mass ~ gdd | EU, data=grow10)

fit.beta.10 <- nlsList(mass ~ SSbgf(gdd, w.max, t.e, t.m), data = grow10G)

fit.nlme.10<-nlme(fit.beta.10, random=pdDiag(w.max ~1))

cfs <- fixef(fit.nlme.10)
fit.nlme2 <- update(fit.nlme.10, fixed = list(t.e ~ trt, w.max + t.m ~ 1),
                    start = c(cfs[2], rep(0,2), cfs[c(1,3)]))

cfsT <- fixef(fit.nlme2)


fit.nlme3 <- update(fit.nlme.10, fixed = list(t.e ~ ground, w.max + t.m ~ 1),
                    start = c(cfs[2], rep(0,1), cfs[c(1,3)]))

cfsG <- fixef(fit.nlme3)

fit.nlme4 <- update(fit.nlme.10, fixed = list(t.e ~ trt + ground, w.max + t.m ~ 1),
                    start = c(cfsT[1:2], cfsG[1:2], cfs[c(1,3)]))

cfsTD <- fixef(fit.nlme4)

fit.nlme5 <- update(fit.nlme4, fixed = list(t.e ~ trt + ground + trt:ground, w.max +     
t.m ~ 1),
                    start = c(cfsTD[1:4], rep(0,2) ,cfsTD[5:6]))


fit.nlme6 <- update(fit.nlme5, weights = varPower())
4

1 に答える 1

0

クリスが言ったように、これは十分な情報ではありません。非常に多くのポイントがゼロに近いので、あなたは与えるべきです

varConstPower(power=0.5,fixed=list(const=10))

試して、オフセットで遊んでください。

于 2012-08-17T17:50:21.223 に答える