Importação e preparo dos dados
# Leitura do área de transferência.
# tb <- read.table("clipboard", header = TRUE, sep = "\t", dec = ",")
# dput(round(tb[, -2], digits = 3))
tb <- data.frame(
kPa = c(1.01, 2, 4, 6, 10, 33, 100, 500, 1500, 9800, 2230, 1060,
15570, 26640, 13930, 18330, 38530, 39620, 36660, 49540,
51640, 60320, 34190, 65420, 4700, 3490, 1400, 2660, 24300,
17140, 24580, 29830, 59310, 27190, 98180, 25810, 97690,
78300, 69730, 48380, 2270, 18680, 13090, 8710, 11230, 24040,
11120, 18850, 20260, 17650, 54920, 49420, 44570, 81680,
40760, 59310, 69290, 9300, 11620, 12190, 10400, 28890,
19820, 68830, 65830, 52500, 47270, 45060, 20060, 90290,
12250, 92750, 82960, 12250, 92750, 82960),
umid = c(56.981, 50.997, 44.2, 42.835, 37.112, 36.112, 33.37,
28.201, 26.994, 14.642, 17.987, 16.192, 6.903, 4.112,
11.407, 9.727, 8.687, 4.676, 4.987, 7.666, 3.475, 1.925,
1.625, 2.15, 41.209, 16.878, 16.843, 17.641, 6.969, 9.219,
3.015, 3.829, 4.447, 8.929, 3.528, 8.813, 3.225, 2.35,
3.05, 0.314, 33.79, 13.655, 32.958, 23.747, 22.801, 8.306,
5.877, 15.81, 6.711, 3.782, 2.686, 4.177, 2.532, 1.1,
1.147, 3.994, 0.827, 27.964, 20.717, 25.672, 24.259, 5.409,
7.238, 3.571, 3.686, 2.141, 3.979, 1.621, 3.971, 0.938,
2.679, 1.824, 3.143, 2.679, 1.824, 3.143))
# Gráfico dos dados.
plot(umid ~ kPa, data = tb, log = "x")
Ajuste do modelo não linear
# Tabela com as variáveis.
tb_fit <- transform(tb,
x = log10(kPa),
y = umid)
head(tb_fit)
## kPa umid x y
## 1 1.01 56.981 0.004321374 56.981
## 2 2.00 50.997 0.301029996 50.997
## 3 4.00 44.200 0.602059991 44.200
## 4 6.00 42.835 0.778151250 42.835
## 5 10.00 37.112 1.000000000 37.112
## 6 33.00 36.112 1.518513940 36.112
# Expressão do modelo.
model <- y ~ res +
(pmp - res)/(1 + (atex * x)^ntex)^(1 - 1/ntex) +
(sat - pmp)/(1 + (aest * x)^nest)^(1 - 1/nest)
# Valores iniciais (passo mais crÃtico do ajuste).
start <- c(res = -15.33,
pmp = 2.73,
sat = 56.82,
atex = 0.24,
ntex = 38.93,
aest = 1.79,
nest = 1.5)
# Ajusta o modelo não linear aos dados.
fit <- nls(formula = model, data = tb_fit, start = start)
# Estimativas dos parâmetros.
summary(fit)
##
## Formula: y ~ res + (pmp - res)/(1 + (atex * x)^ntex)^(1 - 1/ntex) + (sat -
## pmp)/(1 + (aest * x)^nest)^(1 - 1/nest)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## res -15.325580 62.098705 -0.247 0.8058
## pmp 2.729099 66.017032 0.041 0.9671
## sat 56.823822 5.432517 10.460 7.07e-16 ***
## atex 0.240769 0.002964 81.225 < 2e-16 ***
## ntex 38.927190 16.900771 2.303 0.0243 *
## aest 1.788087 1.271388 1.406 0.1641
## nest 1.497447 1.088147 1.376 0.1732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.511 on 69 degrees of freedom
##
## Number of iterations to convergence: 13
## Achieved convergence tolerance: 5.902e-06
# Valores da curva estimada.
pred <- with(tb_fit,
data.frame(x = seq(min(x), max(x), length.out = 151)))
pred$y_fit <- predict(fit, newdata = pred)
# Curva sobre os dados.
plot(y ~ x, data = tb_fit,
xlab = expression("Tensão" ~ (Psi * ", " * log[10])),
ylab = "Conteúdo de água (%)")
lines(y_fit ~ x, data = pred, col = "red")
# Breve análise dos resÃduos.
plot(residuals(fit) ~ fitted(fit))
qqnorm(residuals(fit))
#