r - 複雑な方程式を R 式に入れる方法は?

Question

木の直径を予測子として、木の高さを従属変数として使用します。この種のデータにはさまざまな方程式が存在し、それらのいくつかをモデル化して結果を比較しようとしています。

R formulaただし、1 つの方程式を対応する形式に正しく配置する方法を理解することはできません。

のtreesデータ セットをR例として使用できます。

data(trees)
df <- trees
df$h <- df$Height * 0.3048   #transform to metric system
df$dbh <- (trees$Girth * 0.3048) / pi   #transform tree girth to diameter

最初に、うまく機能すると思われる方程式の例を示します。

form1 <- h ~ I(dbh ^ -1) + I( dbh ^ 2)  
m1 <- lm(form1, data = df)
m1

Call:
lm(formula = form1, data = df)

Coefficients:
(Intercept)    I(dbh^-1)     I(dbh^2)  
27.1147      -5.0553       0.1124  

係数a、bおよびcが推定されます。これは、私たちが関心を持っていることです。

問題のある方程式は次のとおりです。

このように合わせようとしています：

form2 <- h ~ I(dbh ^ 2) / dbh + I(dbh ^ 2) + 1.3

エラーが発生します：

m1 <- lm(form2, data = df)
Error in terms.formula(formula, data = data) 
invalid model formula in ExtractVars

/これは、算術演算子ではなく、ネストされたモデルとして解釈されるためだと思いますか?

これはエラーになりません：

form2 <- h ~ I(I(dbh ^ 2) / dbh + I(dbh ^ 2) + 1.3)
m1 <- lm(form2, data = df)

しかし、結果は私たちが望むものではありません:

m1
Call:
lm(formula = form2, data = df)

Coefficients:
(Intercept)  I(I(dbh^2)/dbh + I(dbh^2) + 1.3)  
19.3883                            0.8727  

I()論理のように、outer 内の項全体に対して 1 つの係数のみが与えられます。

2 番目の方程式をデータにどのように当てはめることができるでしょうか?

Answer 1

いくつか問題があります。（1）の分母の括弧がform2ありません（そしてRには、分母に定数を追加したいことa、または実際にパラメーターのどこに配置するかを知る方法がありません）、そしてはるかに問題があります：（2 ）2番目のモデルは線形ではないため、機能しlmません。

修正（1）は簡単です：

form2 <- h ~ 1.3 + I(dbh^2) / (a + b * dbh + c * I(dbh^2))

（2）を修正すると、非線形モデルのパラメーターを推定する方法はたくさんありますが、nls（非線形最小二乗）から始めるのが適切です。

m2 <- nls(form2, data = df, start = list(a = 1, b = 1, c = 1))

のパラメータの開始推測を提供する必要がありますnls。私はちょうど1を選びました、しかしあなたはパラメータが何であるかもしれないか球場であるより良い推測を使うべきです。

Answer 2

R式を使用していると仮定するとnls、通常のR関数を使用できるH(a, b, c, D)ため、式は正しくh ~ H(a, b, c, dbh)なり、これは機能します。

# use lm to get startingf values
lm1 <- lm(1/(h - 1.3) ~ I(1/dbh) + I(1/dbh^2), df)
start <- rev(setNames(coef(lm1), c("c", "b", "a")))

# run nls
H <- function(a, b, c, D) 1.3 + D^2 / (a + b * D + c * D^2)
nls1 <- nls(h ~ H(a, b, c, dbh), df, start = start)

nls1 # display result

出力のグラフ化：

plot(h ~ dbh, df)
lines(fitted(nls1) ~ dbh, df)

Answer 3

edit :修正され、オフセットを誤って使用しなくなりました ...

@shujaaを補完する答え：

あなたはあなたの問題をから変えることができます

H = 1.3 + D^2/(a+b*D+c*D^2)

に

1/(H-1.3) = a/D^2+b/D+c

これは通常、モデルの仮定を台無しにします (つまり、H一定の分散で正規分布している場合は1/(H-1.3)そうではありません。しかし、とにかく試してみましょう:

data(trees)
df <- transform(trees,
            h=Height * 0.3048,   #transform to metric system
            dbh=Girth * 0.3048 / pi   #transform tree girth to diameter
            )
lm(1/(h-1.3) ~ poly(I(1/dbh),2,raw=TRUE),data=df)

## Coefficients:
##                    (Intercept)  poly(I(1/dbh), 2, raw = TRUE)1  
##                       0.043502                       -0.006136  
## poly(I(1/dbh), 2, raw = TRUE)2  
##                       0.010792  

これらの結果は、通常、適合の適切な開始値を取得するのに十分nlsです。glmただし、リンク関数を使用していくつかの形式の非線形性を可能にする を使用すると、それよりも優れた結果を得ることができます。具体的には、

(fit2 <- glm(h-1.3 ~ poly(I(1/dbh),2,raw=TRUE),
             family=gaussian(link="inverse"),data=df))

## Coefficients:
##                    (Intercept)  poly(I(1/dbh), 2, raw = TRUE)1  
##                       0.041795                       -0.002119  
## poly(I(1/dbh), 2, raw = TRUE)2  
##                       0.008175  
## 
## Degrees of Freedom: 30 Total (i.e. Null);  28 Residual
## Null Deviance:       113.2 
## Residual Deviance: 80.05     AIC: 125.4 
## 

結果は線形近似とほぼ同じですが、完全ではないことがわかります。

pframe <- data.frame(dbh=seq(0.8,2,length=51))

を使用predictしますが、LHS から定数を引いたという事実を説明するために予測を修正する必要があります。

pframe$h <- predict(fit2,newdata=pframe,type="response")+1.3
p2 <- predict(fit2,newdata=pframe,se.fit=TRUE) ## predict on link scale
pframe$h_lwr <- with(p2,1/(fit+1.96*se.fit))+1.3
pframe$h_upr <- with(p2,1/(fit-1.96*se.fit))+1.3
png("dbh_tmp1.png",height=4,width=6,units="in",res=150)
par(las=1,bty="l")
plot(h~dbh,data=df)
with(pframe,lines(dbh,h,col=2))
with(pframe,polygon(c(dbh,rev(dbh)),c(h_lwr,rev(h_upr)),
      border=NA,col=adjustcolor("black",alpha=0.3)))
dev.off()

LHS で定数を使用したため (これはオフセットを使用するフレームワークにほぼ適合しますが、完全には適合しません。式がの場合、つまり、定数調整がリンク上にある場合にのみオフセット1/H - 1.3 = a/D^2 + ...を使用できます (逆）元のスケールではなくスケール）、これはggplotのgeom_smoothフレームワークに完全には適合しません

library("ggplot2")
ggplot(df,aes(dbh,h))+geom_point()+theme_bw()+
   geom_line(data=pframe,colour="red")+
   geom_ribbon(data=pframe,colour=NA,alpha=0.3,
             aes(ymin=h_lwr,ymax=h_upr))

ggsave("dbh_tmp2.png",height=4,width=6)