Блог

У меня есть набор данных о ВЫСОТЕ и ДИАМЕТРЕ деревьев. Я хочу найти регрессионную зависимость между ними и построить ее. Например, я хочу попробовать a * DIAMETER b * DIAMETER^2 C и показать его кривую на диаграмме рассеяния. Следуя приведенной ниже инструкции, я достигаю нескольких строк, но мне нужна только линия тренда, связанная с разработанной моделью. что я должен делать?

 setwd('D:\PhD\Data\Field Measurments\Data Analysis\')

dat1 = read.table('Fagus.csv', header = TRUE, sep =',')


# fit a non-linear regression

Height = dat1$Height

Diameter = dat1$Diameter 

plot(Diameter, Height, main="Height Curve", xlab="Diameter", ylab="Height", pch=19)

nls1 <- nls(Height ~ a*(Diameter)^2 b*Diameter c, data = dat1, start = list(a =a, b=b,c=c), algorithm="port")

lines(fitted(nls1) ~ Diameter, lty = 1, col = "red") # solid red line
  

Является ли приведенная выше инструкция неправильной для моей цели?

Как указано выше, вы не должны вводить коэффициенты в свои формулы. Попробуйте:

 nls1 <- nls(Height ~ I(Diameter^2)   Diameter, data = dat1,  algorithm="port")
  

Что касается I(Diameter ^2) :
«Чтобы избежать этой путаницы, функция I() может использоваться для заключения в скобки тех частей формулы модели, где операторы используются в их арифметическом смысле. Например, в формуле y ~ a I(b c) термин b c следует интерпретировать как сумму b и c.» ~ formula{stats} документация

Я не запускал остальное (на мобильных устройствах), но ваш код на первый взгляд выглядит нормально.

Похоже, здесь существует недопонимание относительно линейных и нелинейных моделей. Линейная модель является линейной по коэффициентам. Нелинейная модель таковой не является. Является ли модель линейной по предикторным переменным (диаметр в вашем случае), не имеет значения. Итак, в вашем случае модель вида:

Высота = a * Диаметр b * Диаметр ^2 c

является линейной моделью. Вам не нужно использовать nls(...) . Вы можете указать формулу модели любым из двух способов, оба из которых приводят к идентичным результатам:

 Height~Diameter   I(Diameter^2)
  

или

 Height~poly(Diameter,2,raw=TRUE)
  

Вторая форма использует poly(...) функцию для создания многочлена порядка 2. raw=T указывает poly(...) генерировать необработанные многочлены, а не ортогональные многочлены (по умолчанию). Первая форма немного проще, если только вам не нужны полиномы порядка больше 2. Вот пример использования обеих форм.

 set.seed(1)    # for reproducible example
df <- data.frame(Diameter=sample(1:50,50))
df$Height <- with(df,2*Diameter   .5*Diameter^2   4   rnorm(50,sd=30))

fit <- lm(Height~Diameter   I(Diameter^2),df)
summary(fit)
# ...
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   -6.85088   12.26720  -0.558  0.57917    
# Diameter       3.31030    1.10964   2.983  0.00451 ** 
# I(Diameter^2)  0.47717    0.02109  22.622  < 2e-16 ***

fit.poly<- lm(Height~poly(Diameter,2,raw=TRUE),df)
summary(fit.poly)
# Coefficients:
#                                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                    -6.85088   12.26720  -0.558  0.57917    
# poly(Diameter, 2, raw = TRUE)1  3.31030    1.10964   2.983  0.00451 ** 
# poly(Diameter, 2, raw = TRUE)2  0.47717    0.02109  22.622  < 2e-16 ***
  

Для построения данных и кривой тренда:

 df$pred <- predict(fit)
with(df,plot(Height~Diameter))
with(df[order(df$Diameter),],lines(pred~Diameter,col="red",lty=2))
  

Ваша проблема заключается в вашем start= параметре. Вам необходимо указать фактические значения для a , b и c параметров. Вот воспроизводимый пример

 #sample data
dat<-data.frame(Diameter = runif(50, 2, 6))
dat<-transform(dat,Height=2*Diameter   .75 * Diameter^2  4   rnorm(50))
dat<-dat[order(dat$Diameter), ]

#now fit the model
mynls<-nls(Height ~ a*I(Diameter^2)   b*Diameter   c, dat,
    start=list(a=1, b=1, c=1), algorithm="port")
  

Обратите внимание, как мы устанавливаем значения по умолчанию, равные 1 для каждого из коэффициентов. Вы можете задать все, что считаете наиболее подходящим. И как мы можем построить необработанные значения с подобранными результатами

 plot(Height~Diameter,dat, main="Height Curve", 
    xlab="Diameter", ylab="Height", pch=19)
lines(fitted(mynls)~ dat$Diameter, col="red")
  

Это дает