library("nls2")
## Loading required package: proto
dados<-data.frame(
x = rep(c( 5, 10, 15, 17, 20, 22, 24, 25, 27, 30, 35), each=4 ),
y = c(1.2, 1.1, 1.5, 1.9,
4.1, 3.95, 3.75, 3.4,
5.15, 4.6, 4.75, 4.15,
3.7, 4.1, 3.9, 4.5,
7, 6.7, 6.7, 6.95,
7, 4.45, 6.15, 6.4,
7, 6.6, 6.7, 7,
4.5, 4.7, 5.75, 4.35,
5.4, 5.15, 5.7, 5.7,
0.001, 0.001, 0.5,
0.001, 0.001, 0.001, 0.001, 0.001)
)
str(dados)
## 'data.frame': 44 obs. of 2 variables:
## $ x: num 5 5 5 5 10 10 10 10 15 15 ...
## $ y: num 1.2 1.1 1.5 1.9 4.1 3.95 3.75 3.4 5.15 4.6 ...
with(dados, c(plot(y~x, xaxt="n", xlim=c(0,35)),
axis(1, at = round(x,0), las=2, cex=0.3),
abline(v=x, lty=3, col="grey"),
lines(lowess(x=x, y=y), lwd=1)))
## [1] 5 5 5 5 10 10 10 10 15 15 15 15 17 17 17 17 20 20 20 20 22 22 22
## [24] 22 24 24 24 24 25 25 25 25 27 27 27 27 30 30 30 30 35 35 35 35
- Criar a função desejada considerando nossos dados.
f_full <-with(dados,
y ~ Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1
)
- Criar um grid com um set de possiveis starters
st1 <- expand.grid(Yopt = seq(4, 8, len = 4),
Tmin = seq(0, 5, len = 4),
Topt = seq(15, 25, len = 4),
Tmax= seq(28, 38, len = 4),
b1 = seq(0, 4, len = 4))
mod_full <- nls2(f_full, start = st1, algorithm = "brute-force")
summary(mod_full)
##
## Formula: y ~ Yopt * ((x - Tmin)/(Topt - Tmin))^(b1 * (Topt - Tmin)/(Tmax -
## Topt)) * ((Tmax - x)/(Tmax - Topt))^b1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Yopt 6.6667 0.4127 16.154 < 2e-16 ***
## Tmin 0.0000 8.8976 0.000 1.000000
## Topt 18.3333 0.9408 19.487 < 2e-16 ***
## Tmax 38.0000 9.1180 4.168 0.000165 ***
## b1 2.6667 2.9963 0.890 0.378931
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.432 on 39 degrees of freedom
##
## Number of iterations to convergence: 1024
## Achieved convergence tolerance: NA
round(coef(mod_full),2)
## Yopt Tmin Topt Tmax b1
## 6.67 0.00 18.33 38.00 2.67
- Diagnostico grafico de nosso ajuste
#plot(mod_full) # Ver residuos
- Plotear os dados com o modelo ajustado
with(dados, c(plot(y~x))); points(fitted(mod_full)~I(dados$x), pch=19)
## NULL
with(as.list(coef(mod_full)),
curve(
Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1,
add=TRUE, col="red"))
Dados reduzidos
- Tirar o “#” daqueles niveis de x que queremos excluir, e o “&” ao primeiro
newdat <- subset(dados, #x!=5
#& x!=10
#& x!=15
#x!=17
# x!=20
#x!=22
# x!=24
# x!=25
x!=27
#& x!=30
&x!=35
)
with(newdat, c(plot(y~x, xaxt="n"),
axis(1, at = round(x,0), las=2, cex=0.3),
abline(v=x, lty=3, col="grey"),
lines(lowess(x=x, y=y), lwd=1)))
## [1] 5 5 5 5 10 10 10 10 15 15 15 15 17 17 17 17 20 20 20 20 22 22 22
## [24] 22 24 24 24 24 25 25 25 25 30 30 30 30
- Editar a função desejada considerando o novos dados.
f_redu <-with(newdat,
y ~ Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1
)
mod_redu <- nls2(f_redu, start = st1, algorithm = "brute-force")
summary(mod_redu)
##
## Formula: y ~ Yopt * ((x - Tmin)/(Topt - Tmin))^(b1 * (Topt - Tmin)/(Tmax -
## Topt)) * ((Tmax - x)/(Tmax - Topt))^b1
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Yopt 6.6667 0.3847 17.327 < 2e-16 ***
## Tmin 0.0000 13.0295 0.000 1.000
## Topt 21.6667 0.8187 26.464 < 2e-16 ***
## Tmax 31.3333 2.3438 13.369 2.07e-14 ***
## b1 1.3333 1.3294 1.003 0.324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.227 on 31 degrees of freedom
##
## Number of iterations to convergence: 1024
## Achieved convergence tolerance: NA
round(coef(mod_redu),2)
## Yopt Tmin Topt Tmax b1
## 6.67 0.00 21.67 31.33 1.33
- Plotear os dados com o modelo ajustado
with(newdat, c(plot(y~x))); points(fitted(mod_redu)~I(newdat$x), pch=19)
## NULL
with(as.list(coef(mod_redu)),
curve(
Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1,
add=TRUE, col="red"))
