Jak umieścić skomplikowane równanie w formule R?

Question

Mamy średnicę drzew jako predyktor i wysokość drzewa jako zmienną zależną. Istnieje wiele różnych równań dla tego rodzaju danych i staramy się modelować niektóre z nich i porównywać wyniki.

Jednak nie możemy dowiedzieć się, jak poprawnie umieścić jedno równanie w odpowiednim formacie Rformula.

Jako przykład można użyć zestawu danych trees w 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 

pierwsze, przykładem równania, które wydaje się działać dobrze.

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 

Współczynniki a, b i c szacuje, czyli to, co nas interesuje

Teraz problematyczne równanie:

Próbując dopasować go tak:

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

daje błąd:

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

to chyba dlatego / jest interpretowany jako zagnieżdżonego modelu, a nie operator arytmetyczny ?

To nie daje błąd:

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

Ale wynik nie jest jeden chcemy:

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

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

Tylko jeden współczynnik jest podawany dla całego wyrażenia wewnątrz zewnętrznej I(), co wydaje się być logiczne.

Jak możemy dopasować drugie równanie do naszych danych?

Answer 1

Zakładając, że za pomocą nls wzorze R może użyć zwykłej funkcji R, H(a, b, c, D), więc formuła może być tylko h ~ H(a, b, c, dbh) i to działa:

# 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 

Wykresy wyjście:

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

Answer 2

Masz kilka problemów. (1) Brakuje nawiasów w mianowniku form2 (a R nie może wiedzieć, że chcesz dodać stałe w metodzie stałej a lub gdzie umieścić którykolwiek z parametrów) i znacznie bardziej problematyczne: (2) twój drugi model nie jest liniowy, więc lm nie będzie działać.

mocujący (1) jest prosta:

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

Mocowanie (2), choć istnieje wiele sposobów szacowania parametrów nieliniowego modelu nls (nieliniowe najmniejszych kwadratów) jest dobrym miejscem, aby rozpocząć:

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

Należy podać początkowe wartości domyślne dla parametrów w nls.Właśnie wybrałam 1, ale powinieneś użyć lepszych przypuszczeń, które mogą być tymi parametrami.

Answer 3

edit: stałe, nie jest już nieprawidłowo offsetowego ...

odpowiedź, która uzupełnia @ shujaa użytkownika:

Można przekształcić swój problem z

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

do

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

Normalnie zepsułoby to założenia modelu (to znaczy, gdyby H były normalnie dystrybuowane ze stałą wariancją, to nie byłoby 1/(H-1.3). Jednak spróbujmy to tak:

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 

Wyniki te powinny być one wystarczająco dobre, aby uzyskać dobre wartości wyjściowych dla nls dopasowanie. Można jednak zrobić to lepiej za pomocą glm, która wykorzystuje funkcję łącza, aby umożliwić pewne formy nieliniowości. Konkretnie

(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 
## 

Widać, że wyniki są około taka sama jak dopasowanie liniowe, ale nie całkiem.

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

Używamy predict, ale trzeba skorygować prognozy w celu uwzględnienia faktu, że odejmuje stałe z 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() 

Ponieważ użyliśmy stała na LHS (to prawie, ale nie całkiem, pasuje do schematu użycia przesunięcie - możemy użyć tylko offsetu, jeśli nasza formuła byłaby 1/H - 1.3 = a/D^2 + ..., czyli jeśli Stała dostosowawcze odnośnik (odwrotny) skalowania zamiast oryginalnej skali), nie jest idealnie pasuje do ggplot jest geom_smooth ramach

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)