आर सूत्र में जटिल समीकरण कैसे डालें?

Question

हमारे पास आश्रित और वृक्ष की ऊंचाई पर निर्भर चर के रूप में पेड़ का व्यास है। इस तरह के डेटा के लिए कई अलग-अलग समीकरण मौजूद हैं और हम उनमें से कुछ को मॉडल करने और परिणामों की तुलना करने की कोशिश करते हैं।

हालांकि, हम यह नहीं समझ सकते कि एक समीकरण को Rformula प्रारूप में सही तरीके से कैसे रखा जाए।

R में trees डेटा सेट को उदाहरण के रूप में उपयोग किया जा सकता है।

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 का अनुमान है, जो है क्या हम में रुचि रखते हैं

अब समस्याग्रस्त समीकरण:

इस तरह यह फिट करने के लिए कोशिश कर रहा है:

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

मुझे लगता है कि इस वजह से / और नहीं एक नेस्टेड मॉडल के रूप में व्याख्या की है एक अंकगणितीय ऑपरेटर:

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

त्रुटि देता ?

यह एक त्रुटि नहीं देता:

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() भीतर पूरी अवधि के लिए दिया जाता है, जो तर्क लगता है।

हम अपने डेटा के दूसरे समीकरण को कैसे फिट कर सकते हैं?

Answer 1

में पूरी तरह से फिट नहीं करता है मान लिया जाये कि आप nls आर सूत्र का उपयोग कर रहे हैं कर सकते हैं एक साधारण आर समारोह, 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 2

आपको कुछ समस्याएं हैं। (1) form2 के denominator के लिए आप कोष्ठक याद आ रहे हैं (और आर को यह जानने का कोई तरीका नहीं है कि आप denominator में निरंतर a जोड़ना चाहते हैं, या वास्तव में किसी भी पैरामीटर को कहां रखना है), और बहुत अधिक समस्याग्रस्त: (2) आपका दूसरा मॉडल रैखिक नहीं है, इसलिए lm काम नहीं करेगा।

फिक्सिंग (1) आसान है:

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

फिक्सिंग (2), हालांकि वहाँ एक nonlinear मॉडल के लिए मानकों का अनुमान लगाने के लिए कई तरीके हैं, nls (nonlinear कम से कम वर्गों) एक अच्छी जगह शुरू करने के लिए है:

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

आपको nls में पैरामीटर के लिए प्रारंभिक अनुमान प्रदान करने की आवश्यकता है।मैंने अभी 1 को चुना है, लेकिन आपको बेहतर अनुमानों का उपयोग करना चाहिए जो पैरामीटर हो सकते हैं।

Answer 3

संपादित: तय, अब गलत तरीके से ऑफसेट का उपयोग कर ...

एक जवाब है कि पूरक @ 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 उपयोग करें, लेकिन भविष्यवाणी तथ्य को ध्यान में है कि हम एलएचएस से एक निरंतर घटाया सही करने की आवश्यकता:

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

क्योंकि हम पर लगातार इस्तेमाल किया है एलएचएस (यह लगभग, लेकिन ऑफसेट का उपयोग करने के ढांचे में फिट नहीं है, लेकिन अगर हम फॉर्मूला 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)