Come inserire un'equazione complicata in una formula R?

Question

Abbiamo il diametro degli alberi come predittore e altezza dell'albero come variabile dipendente. Esistono diverse equazioni per questo tipo di dati e cerchiamo di modellarne alcuni e confrontare i risultati.

Tuttavia, non riusciamo a capire come inserire correttamente un'equazione nel corrispondente formato Rformula.

Il set di dati trees in R può essere utilizzato come esempio.

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 

In primo luogo, l'esempio di un'equazione che sembra funzionare bene:.

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 

Coefficienti a, b e c sono stimati, che è quello che ci interessa

Ora l'equazione problematica:

Cercando di adattarsi in questo modo:

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

dà un errore:

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

Credo che questo è dovuto al fatto / viene interpretato come un modello nidificato e non un operatore aritmetico ?

questo non dà un errore:

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

Ma il risultato non è quello che vogliamo:

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

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

un solo coefficiente è dato per l'intera durata all'interno del esterna I(), che sembra essere logico.

Come possiamo adattare la seconda equazione ai nostri dati?

Answer 1

Supponendo che si sta utilizzando nls la formula R può utilizzare una funzione ordinaria R, H(a, b, c, D), quindi la formula può essere solo h ~ H(a, b, c, dbh) e questo funziona:

# 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 

graficamente l'uscita:

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

Answer 2

Hai un paio di problemi. (1) Ti manca parentesi per il denominatore di form2 (e R non ha modo di sapere che vuoi aggiungere una costante a al denominatore, o dove mettere uno dei parametri, in realtà), e molto più problematico: (2) il tuo 2o modello non è lineare, quindi lm non funzionerà.

di fissaggio (1) è facile:

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

di fissaggio (2), anche se ci sono molti modi per stimare i parametri per un modello non lineare, il nls (non lineari minimi quadrati) è un buon punto di partenza:

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

è necessario fornire a partire ipotesi per i parametri in nls.Ho scelto 1, ma dovresti usare una migliore ipotesi per capire quali potrebbero essere i parametri.

Answer 3

modificare: fisso, non è più in modo non corretto utilizzando compensati ...

Una risposta che integra @ Shujaa di:

È possibile trasformare il vostro problema da

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

a

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

Questo normalmente rovinerebbe le ipotesi del modello (ad esempio, se H fosse distribuito normalmente con varianza costante, quindi non sarebbe 1/(H-1.3). Tuttavia, proviamo lo stesso:

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 

Questi risultati sarebbero normalmente abbastanza buono per ottenere buoni valori di partenza per la nls forma. Tuttavia, puoi fare meglio di questo tramite glm, che utilizza una funzione di collegamento per consentire alcune forme di non-linearità. In particolare,

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

si può vedere che i risultati sono circa la stessa della forma lineare, ma non del tutto.

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

Usiamo predict, ma necessità di correggere la previsione per tenere conto del fatto che abbiamo sottratto un costante dal 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() 

Poiché abbiamo usato il costante il LHS (questo quasi, ma non del tutto, si inserisce nella struttura dell'utilizzo di uno scostamento - potremmo usare solo uno scostamento se la nostra formula fosse 1/H - 1.3 = a/D^2 + ..., cioè se il costante adeguamento erano sul link (inverso) scala piuttosto che la scala originale), questo non si adatta perfettamente nel quadro geom_smoothggplot s'

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)