Microwaves and Ultrasonic energy for assessment physical attributes in highly oxidic soils (Brazil)
Samara Barbosa, Carla E. Carducci & Walmes M. Zeviani
Importação e análise exploratória
# Carrega os pacotes necessários.
library(nlme)
library(tidyverse)txt <- "nome soil time rep umid umrel
CXa-A 1 4 1 4.675505932 12.022705
CXa-A 1 4 2 3.806228374 9.787424411
CXa-A 1 8 1 6.837606838 17.58238169
CXa-A 1 8 2 14.06844106 36.17591746
CXa-A 1 12 1 15.5624037 40.01752782
CXa-A 1 12 2 11.35857461 29.20770366
CXa-A 1 16 1 17.83189317 45.85334597
CXa-A 1 16 2 22.74959083 58.49882845
CXa-A 1 20 1 22.24938875 57.21259715
CXa-A 1 20 2 20.38523274 52.41906291
CXa-A 1 24 1 31.92612137 82.09557309
CXa-A 1 24 2 29.53367876 75.94359035
CXa-A 1 28 1 29.31034483 75.36930428
CXa-A 1 28 2 29.98266898 77.09813428
CXa-A 1 32 1 32.74336283 84.19704684
CXa-A 1 32 2 31.34851138 80.61029329
CXa-A 1 36 1 36.36363636 93.50630263
CXa-A 1 36 2 38.12154696 98.02663494
CXa-A 1 40 1 35.74660633 91.91965723
CXa-A 1 40 2 38.76040703 99.66941462
CXbdf-A 2 4 1 16.18900077 40.767766
CXbdf-A 2 4 2 9.011627907 22.69342889
CXbdf-A 2 8 1 20.28869286 51.0917687
CXbdf-A 2 8 2 15.03067485 37.85082498
CXbdf-A 2 12 1 21.1631664 53.29390168
CXbdf-A 2 12 2 27.44265081 69.10714146
CXbdf-A 2 16 1 25.7334451 64.80295373
CXbdf-A 2 16 2 29.64563526 74.65478186
CXbdf-A 2 20 1 34.64991023 87.25674006
CXbdf-A 2 20 2 30.77593723 77.50115185
CXbdf-A 2 24 1 37.99448022 95.67916517
CXbdf-A 2 24 2 34.28827216 86.34604911
CXbdf-A 2 28 1 39.92537313 100.5416142
CXbdf-A 2 28 2 38.76040703 97.60795159
CXbdf-A 2 32 1 39.40520446 99.23170534
CXbdf-A 2 32 2 39.27576602 98.90574846
CXbdf-A 2 36 1 38.50415512 96.96264816
CXbdf-A 2 36 2 40.58106842 102.1928113
CXbdf-A 2 40 1 38.76040703 97.60795159
CXbdf-A 2 40 2 40.97744361 103.1909786
CXbdf-C 3 4 1 9.80966325 24.3832721
CXbdf-C 3 4 2 5.708245243 14.1886315
CXbdf-C 3 8 1 12.10762332 30.09516905
CXbdf-C 3 8 2 12.61261261 31.35038965
CXbdf-C 3 12 1 15.11895625 37.58025274
CXbdf-C 3 12 2 15.11895625 37.58025274
CXbdf-C 3 16 1 26.15643398 65.01542719
CXbdf-C 3 16 2 19.61722488 48.76132033
CXbdf-C 3 20 1 29.19896641 72.57806153
CXbdf-C 3 20 2 27.01100762 67.13958795
CXbdf-C 3 24 1 30.77593723 76.49784018
CXbdf-C 3 24 2 24.58471761 61.1087092
CXbdf-C 3 28 1 35.74660633 88.85312438
CXbdf-C 3 28 2 36.98630137 91.93455751
CXbdf-C 3 32 1 36.98630137 91.93455751
CXbdf-C 3 32 2 36.98630137 91.93455751
CXbdf-C 3 36 1 39.79496738 98.91588458
CXbdf-C 3 36 2 36.73655424 91.31377655
CXbdf-C 3 40 1 39.27576602 97.62533792
CXbdf-C 3 40 2 38.24884793 95.07279126
LVAd-A 4 4 1 9.489051095 24.31957264
LVAd-A 4 4 2 16.91348402 43.3477172
LVAd-A 4 8 1 16.82242991 43.1143538
LVAd-A 4 8 2 12.78195489 32.75898478
LVAd-A 4 12 1 16.18900077 41.4909327
LVAd-A 4 12 2 19.80830671 50.76688377
LVAd-A 4 16 1 25 64.07272024
LVAd-A 4 16 2 24.17218543 61.95110698
LVAd-A 4 20 1 27.11864407 69.50261178
LVAd-A 4 20 2 27.76831346 71.16765519
LVAd-A 4 24 1 33.33333333 85.43029365
LVAd-A 4 24 2 33.80909902 86.64963771
LVAd-A 4 28 1 33.68983957 86.34398663
LVAd-A 4 28 2 34.77088949 89.11461898
LVAd-A 4 32 1 38.24884793 98.0283093
LVAd-A 4 32 2 37.48854262 96.07971614
LVAd-A 4 36 1 37.11151737 95.11343479
LVAd-A 4 36 2 38.37638376 98.35517202
LVAd-A 4 40 1 38.37638376 98.35517202
LVAd-A 4 40 2 40.05602241 102.6599327
LVAd-Bw 5 4 1 10.29411765 27.48935307
LVAd-Bw 5 4 2 8.538350217 22.80076173
LVAd-Bw 5 8 1 6.157112527 16.44191818
LVAd-Bw 5 8 2 9.249817917 24.70066102
LVAd-Bw 5 12 1 12.02389843 32.10854981
LVAd-Bw 5 12 2 14.94252874 39.90244353
LVAd-Bw 5 16 1 19.71268955 52.64065376
LVAd-Bw 5 16 2 23.05168171 61.55707939
LVAd-Bw 5 20 1 24.27506214 64.82398751
LVAd-Bw 5 20 2 30.3214596 80.97025282
LVAd-Bw 5 24 1 31.00436681 82.79388435
LVAd-Bw 5 24 2 27.98634812 74.73458445
LVAd-Bw 5 28 1 35.6238698 95.12977879
LVAd-Bw 5 28 2 32.04225352 85.5654511
LVAd-Bw 5 32 1 34.52914798 92.20643988
LVAd-Bw 5 32 2 35.50135501 94.802616
LVAd-Bw 5 36 1 36.86131387 98.43424237
LVAd-Bw 5 36 2 37.23696249 99.4373723
LVAd-Bw 5 40 1 36.73655424 98.10108496
LVAd-Bw 5 40 2 37.36263736 99.77297377
LVdf-A 6 4 1 4.968509447 13.17933796
LVdf-A 6 4 2 2.040816327 5.41341591
LVdf-A 6 8 1 5.857445307 15.53730593
LVdf-A 6 8 2 5.782792666 15.33928429
LVdf-A 6 12 1 10.61946903 28.16892526
LVdf-A 6 12 2 11.69024572 31.00923946
LVdf-A 6 16 1 15.29592621 40.57357305
LVdf-A 6 16 2 18.5770751 49.27706261
LVdf-A 6 20 1 22.44897959 59.54757501
LVdf-A 6 20 2 19.90407674 52.79703239
LVdf-A 6 24 1 28.31479897 75.10709379
LVdf-A 6 24 2 25.20868114 66.86788701
LVdf-A 6 28 1 30.20833333 80.12983341
LVdf-A 6 28 2 30.66202091 81.33327318
LVdf-A 6 32 1 35.74660633 94.82051125
LVdf-A 6 32 2 34.28827216 90.95217223
LVdf-A 6 36 1 36.36363636 96.45722894
LVdf-A 6 36 2 35.99274705 95.47341767
LVdf-A 6 40 1 36.48771611 96.78635961
LVdf-A 6 40 2 36.23978202 96.12869614
LVdf-B 7 4 1 5.411103303 13.97568758
LVdf-B 7 4 2 2.459016393 6.35109754
LVdf-B 7 8 1 4.748603352 12.26459618
LVdf-B 7 8 2 6.307583274 16.29109782
LVdf-B 7 12 1 12.61261261 32.57559939
LVdf-B 7 12 2 11.77347243 30.40828519
LVdf-B 7 16 1 13.12217195 33.89167888
LVdf-B 7 16 2 22.95081967 59.27691037
LVdf-B 7 20 1 23.45679012 60.58372057
LVdf-B 7 20 2 20.87026591 53.90329843
LVdf-B 7 24 1 29.08777969 75.12732592
LVdf-B 7 24 2 28.31479897 73.13088704
LVdf-B 7 28 1 32.50883392 83.96315523
LVdf-B 7 28 2 29.42191544 75.99032495
LVdf-B 7 32 1 35.99274705 92.96133522
LVdf-B 7 32 2 32.39187996 83.66108893
LVdf-B 7 36 1 37.6146789 97.15042782
LVdf-B 7 36 2 34.89208633 90.11857109
LVdf-B 7 40 1 36.23978202 93.5993721
LVdf-B 7 40 2 35.74660633 92.32560798
LVu-Bw 9 4 1 17.64705882 51.95064232
LVu-Bw 9 4 2 13.29305136 39.13301154
LVu-Bw 9 8 1 20.96774194 61.72630082
LVu-Bw 9 8 2 23.35526316 68.75485228
LVu-Bw 9 12 1 27.87723785 82.06695671
LVu-Bw 9 12 2 23.15270936 68.15856029
LVu-Bw 9 16 1 29.19896641 85.9579534
LVu-Bw 9 16 2 31.2335958 91.94763729
LVu-Bw 9 20 1 32.62599469 96.04667825
LVu-Bw 9 20 2 32.27513228 95.01378499
LVu-Bw 9 24 1 32.9787234 97.08506562
LVu-Bw 9 24 2 33.80909902 99.52958325
LVu-Bw 9 28 1 34.16815742 100.5866044
LVu-Bw 9 28 2 32.27513228 95.01378499
LVu-Bw 9 32 1 32.50883392 95.70177219
LVu-Bw 9 32 2 32.9787234 97.08506562
LVu-Bw 9 36 1 33.92857143 99.88129447
LVu-Bw 9 36 2 32.39187996 95.35747498
LVu-Bw 9 40 1 33.4519573 98.47820455
LVu-Bw 9 40 2 32.9787234 97.08506562"
tb <- read.table(textConnection(txt), header = TRUE, sep = "\t")
str(tb)## 'data.frame': 160 obs. of 6 variables:
## $ nome : Factor w/ 8 levels "CXa-A","CXbdf-A",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ soil : int 1 1 1 1 1 1 1 1 1 1 ...
## $ time : int 4 4 8 8 12 12 16 16 20 20 ...
## $ rep : int 1 2 1 2 1 2 1 2 1 2 ...
## $ umid : num 4.68 3.81 6.84 14.07 15.56 ...
## $ umrel: num 12.02 9.79 17.58 36.18 40.02 ...ggplot(tb,
aes(x = time, y = umrel)) +
facet_wrap(facets = ~nome) +
geom_point()
Ajuste do modelo aos dados
Gompertz model was assumed to describe the relative water content as a function of time in microwave, \[ f(t; A, b, c) = A \cdot \exp\{-\exp\{-c \cdot (t - b)\}\}, \] where:
- \(t\) is the period in microwave (minutes).
- \(A\) is the asymptote parameter and represents the limit value of the function when \(t\) goes to infinity.
- \(c\) is a shape parameter related to the rate at the inflextion point.
- \(b\) is a shape parameter related to the position of inflexion point.
The pamater of interest is the relative water content for 40 minutes, \(y_{40} = f(t = 40; A, b, c)\). To improve discussion, the confidence interval for the predicted value at \(t = 40\) was determined, so \(\text{IC}(y_{40})\) is the confidence interval for the predicted value at 40 minutes. The null hypotesis is that 40 minutes is enought to achieve 100% of relative water content in soil samples,
\[ H_0: f(t = 40, A, b, c) = 100 \quad vs \quad H_a: f(t = 40, A, b, c) \neq 100. \]
A confidence interval with 95% of coverage is obtained at 40 minutes, \(\text{IC}(y_{40})\). When this interval contains the value 100 is, in terms of hypothesis test, accept the null hypohtesis at 5% of significance.
# Modelo para a curva de secagem.
model <- umrel ~ A * exp(-exp(-C * (time - B)))
# Ajuste do modelo para todos os solos conjuntamente.
m0 <- gnls(model = model,
data = tb,
params = A + C + B ~ nome,
start = c(100, rep(0, times = 7),
0.05, rep(0, times = 7),
10, rep(0, times = 7)))
# Resultados.
summary(m0)## Generalized nonlinear least squares fit
## Model: model
## Data: tb
## AIC BIC logLik
## 1034.202 1111.081 -492.101
##
## Coefficients:
## Value Std.Error t-value p-value
## A.(Intercept) 107.27077 7.499912 14.302937 0.0000
## A.nomeCXbdf-A -1.16148 8.546840 -0.135896 0.8921
## A.nomeCXbdf-C 1.14570 10.124944 0.113156 0.9101
## A.nomeLVAd-A 8.84669 11.858789 0.746003 0.4570
## A.nomeLVAd-Bw 2.65073 9.524089 0.278318 0.7812
## A.nomeLVdf-A 2.76048 10.161453 0.271662 0.7863
## A.nomeLVdf-B -5.06544 9.273899 -0.546204 0.5858
## A.nomeLVu-Bw -8.10827 7.833214 -1.035114 0.3025
## C.(Intercept) 0.08126 0.012548 6.475638 0.0000
## C.nomeCXbdf-A 0.01953 0.019441 1.004681 0.3168
## C.nomeCXbdf-C 0.00091 0.017749 0.051354 0.9591
## C.nomeLVAd-A -0.01424 0.017787 -0.800332 0.4249
## C.nomeLVAd-Bw 0.00912 0.017597 0.518200 0.6052
## C.nomeLVdf-A 0.01019 0.017747 0.574286 0.5667
## C.nomeLVdf-B 0.02006 0.018538 1.082224 0.2811
## C.nomeLVu-Bw 0.07494 0.029650 2.527341 0.0126
## B.(Intercept) 13.00004 1.158629 11.220192 0.0000
## B.nomeCXbdf-A -6.71557 1.351628 -4.968507 0.0000
## B.nomeCXbdf-C -1.77971 1.548607 -1.149231 0.2525
## B.nomeLVAd-A -4.09063 1.757753 -2.327190 0.0214
## B.nomeLVAd-Bw -1.37688 1.445025 -0.952842 0.3424
## B.nomeLVdf-A 1.99342 1.520718 1.310843 0.1921
## B.nomeLVdf-B 0.77965 1.427654 0.546105 0.5859
## B.nomeLVu-Bw -10.48393 1.403063 -7.472174 0.0000
##
## Correlation:
## A.(In) A.CX-A A.CX-C A.LVA-A A.LVA-B A.LV-A A.nmLVd-B A.nmLV-Bw
## A.nomeCXbdf-A -0.878
## A.nomeCXbdf-C -0.741 0.650
## A.nomeLVAd-A -0.632 0.555 0.468
## A.nomeLVAd-Bw -0.787 0.691 0.583 0.498
## A.nomeLVdf-A -0.738 0.648 0.547 0.467 0.581
## A.nomeLVdf-B -0.809 0.710 0.599 0.511 0.637 0.597
## A.nomeLVu-Bw -0.957 0.840 0.709 0.606 0.754 0.707 0.774
## C.(Intercept) -0.911 0.800 0.675 0.576 0.717 0.672 0.737 0.872
## C.nomeCXbdf-A 0.588 -0.826 -0.436 -0.372 -0.463 -0.434 -0.476 -0.563
## C.nomeCXbdf-C 0.644 -0.565 -0.908 -0.407 -0.507 -0.475 -0.521 -0.617
## C.nomeLVAd-A 0.643 -0.564 -0.476 -0.921 -0.506 -0.474 -0.520 -0.615
## C.nomeLVAd-Bw 0.650 -0.570 -0.481 -0.411 -0.895 -0.480 -0.525 -0.622
## C.nomeLVdf-A 0.644 -0.565 -0.477 -0.407 -0.507 -0.903 -0.521 -0.617
## C.nomeLVdf-B 0.617 -0.541 -0.457 -0.390 -0.486 -0.455 -0.875 -0.590
## C.nomeLVu-Bw 0.386 -0.338 -0.286 -0.244 -0.304 -0.285 -0.312 -0.546
## B.(Intercept) 0.862 -0.756 -0.639 -0.545 -0.679 -0.636 -0.697 -0.825
## B.nomeCXbdf-A -0.739 0.732 0.547 0.467 0.582 0.545 0.598 0.708
## B.nomeCXbdf-C -0.645 0.566 0.839 0.408 0.508 0.476 0.522 0.618
## B.nomeLVAd-A -0.568 0.499 0.421 0.845 0.447 0.419 0.460 0.544
## B.nomeLVAd-Bw -0.691 0.607 0.512 0.437 0.830 0.510 0.559 0.662
## B.nomeLVdf-A -0.657 0.576 0.487 0.415 0.517 0.857 0.531 0.629
## B.nomeLVdf-B -0.700 0.614 0.518 0.442 0.551 0.516 0.832 0.670
## B.nomeLVu-Bw -0.712 0.625 0.527 0.450 0.561 0.525 0.576 0.654
## C.(In) C.CX-A C.CX-C C.LVA-A C.LVA-B C.LV-A C.nmLVd-B C.nmLV-Bw
## A.nomeCXbdf-A
## A.nomeCXbdf-C
## A.nomeLVAd-A
## A.nomeLVAd-Bw
## A.nomeLVdf-A
## A.nomeLVdf-B
## A.nomeLVu-Bw
## C.(Intercept)
## C.nomeCXbdf-A -0.645
## C.nomeCXbdf-C -0.707 0.456
## C.nomeLVAd-A -0.705 0.455 0.499
## C.nomeLVAd-Bw -0.713 0.460 0.504 0.503
## C.nomeLVdf-A -0.707 0.456 0.500 0.499 0.504
## C.nomeLVdf-B -0.677 0.437 0.479 0.477 0.483 0.479
## C.nomeLVu-Bw -0.423 0.273 0.299 0.299 0.302 0.299 0.286
## B.(Intercept) -0.705 0.455 0.498 0.497 0.503 0.498 0.477 0.298
## B.nomeCXbdf-A 0.604 -0.386 -0.427 -0.426 -0.431 -0.427 -0.409 -0.256
## B.nomeCXbdf-C 0.527 -0.340 -0.666 -0.372 -0.376 -0.373 -0.357 -0.223
## B.nomeLVAd-A 0.465 -0.300 -0.329 -0.684 -0.331 -0.329 -0.315 -0.197
## B.nomeLVAd-Bw 0.565 -0.365 -0.400 -0.399 -0.642 -0.400 -0.383 -0.239
## B.nomeLVdf-A 0.537 -0.347 -0.380 -0.379 -0.383 -0.694 -0.364 -0.227
## B.nomeLVdf-B 0.572 -0.369 -0.404 -0.404 -0.408 -0.404 -0.631 -0.242
## B.nomeLVu-Bw 0.582 -0.376 -0.412 -0.411 -0.415 -0.412 -0.394 0.091
## B.(In) B.CX-A B.CX-C B.LVA-A B.LVA-B B.LV-A B.nmLVd-B
## A.nomeCXbdf-A
## A.nomeCXbdf-C
## A.nomeLVAd-A
## A.nomeLVAd-Bw
## A.nomeLVdf-A
## A.nomeLVdf-B
## A.nomeLVu-Bw
## C.(Intercept)
## C.nomeCXbdf-A
## C.nomeCXbdf-C
## C.nomeLVAd-A
## C.nomeLVAd-Bw
## C.nomeLVdf-A
## C.nomeLVdf-B
## C.nomeLVu-Bw
## B.(Intercept)
## B.nomeCXbdf-A -0.857
## B.nomeCXbdf-C -0.748 0.641
## B.nomeLVAd-A -0.659 0.565 0.493
## B.nomeLVAd-Bw -0.802 0.687 0.600 0.529
## B.nomeLVdf-A -0.762 0.653 0.570 0.502 0.611
## B.nomeLVdf-B -0.812 0.696 0.607 0.535 0.651 0.618
## B.nomeLVu-Bw -0.826 0.708 0.618 0.544 0.662 0.629 0.670
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.69052326 -0.62932292 -0.06360098 0.62660571 2.53523662
##
## Residual standard error: 5.685687
## Degrees of freedom: 160 total; 136 residual# Quadro de testes para o termos de efeito fixo.
anova(m0)## Denom. DF: 136
## numDF F-value p-value
## A.(Intercept) 1 32143.29 <.0001
## A.nome 7 45.30 <.0001
## C.(Intercept) 1 2554.18 <.0001
## C.nome 7 20.49 <.0001
## B.(Intercept) 1 1050.84 <.0001
## B.nome 7 29.37 <.0001# Reajuste com parâmetros por solo.
m1 <- gnls(model = model,
data = tb,
params = A + C + B ~ 0 + nome,
start = c(rep(100, times = 8),
rep(0.05, times = 8),
rep(10, times = 8)))
# Tabela com as estimativas dos parâmetros para cada solo.
summary(m1)$tTable## Value Std.Error t-value p-value
## A.nomeCXa-A 107.27076636 7.49991202 14.302937 6.811608e-29
## A.nomeCXbdf-A 106.10928886 4.09875455 25.888178 2.072395e-54
## A.nomeCXbdf-C 108.41646596 6.80189745 15.939150 6.453204e-33
## A.nomeLVAd-A 116.11745646 9.18597779 12.640729 1.045723e-24
## A.nomeLVAd-Bw 109.92149493 5.87022988 18.725245 1.799550e-39
## A.nomeLVdf-A 110.03124558 6.85612504 16.048605 3.506418e-33
## A.nomeLVdf-B 102.20532169 5.45495491 18.736236 1.698709e-39
## A.nomeLVu-Bw 99.16249669 2.26065418 43.864514 3.848443e-82
## C.nomeCXa-A 0.08125507 0.01254781 6.475638 1.590747e-09
## C.nomeCXbdf-A 0.10078658 0.01484876 6.787544 3.222896e-10
## C.nomeCXbdf-C 0.08216653 0.01255238 6.545892 1.113731e-09
## C.nomeLVAd-A 0.06701948 0.01260689 5.316101 4.234846e-07
## C.nomeLVAd-Bw 0.09037381 0.01233713 7.325352 1.901350e-11
## C.nomeLVdf-A 0.09144703 0.01255048 7.286337 2.341795e-11
## C.nomeLVdf-B 0.10131747 0.01364604 7.424678 1.116428e-11
## C.nomeLVu-Bw 0.15619137 0.02686428 5.814091 4.142388e-08
## B.nomeCXa-A 13.00003592 1.15862871 11.220191 4.348630e-21
## B.nomeCXbdf-A 6.28446129 0.69604509 9.028813 1.489910e-15
## B.nomeCXbdf-C 11.22032948 1.02750273 10.920000 2.532226e-20
## B.nomeLVAd-A 8.90940970 1.32184595 6.740127 4.117485e-10
## B.nomeLVAd-Bw 11.62315550 0.86352567 13.460116 8.830491e-27
## B.nomeLVdf-A 14.99345852 0.98496802 15.222279 3.618027e-31
## B.nomeLVdf-B 13.77968511 0.83413113 16.519807 2.576412e-34
## B.nomeLVu-Bw 2.51610722 0.79130578 3.179690 1.826021e-03# Resíduos e valores ajustados para inspecionar.
tb$res <- residuals(m0)
tb$fit <- fitted(m0)
# Resíduos contra o tempo.
ggplot(tb,
aes(x = time, y = res)) +
facet_wrap(facets = ~nome) +
geom_point()
# Resíduos contra os valores ajustados.
ggplot(tb,
aes(x = fit, y = res)) +
facet_wrap(facets = ~nome) +
geom_point()
# Normalidade dos resíduos.
ggplot(tb,
aes(sample = res)) +
facet_wrap(facets = ~nome) +
geom_qq()
# Grid para a predição com bandas.
time_seq <- seq(0, 42, by = 0.5)
grid <- with(tb,
expand.grid(time = time_seq,
nome = levels(nome)))
grid$fit <- predict(m1, newdata = grid)
# Para retornar a matriz gradiente para um vetor `time`.
partials <- deriv(model[-2],
namevec = c("A", "C", "B"),
function.arg = function(time, A, C, B) {
NULL
})
# Matriz com os parâmetros estimados por solo.
params <- matrix(coef(m1),
ncol = 3,
dimnames = list(levels(tb$nome),
c("A", "C", "B")))
# Lista de matrizes gradiente.
grad <- sapply(levels(tb$nome),
simplify = FALSE,
FUN = function(nome) {
u <- do.call(partials,
args = c(list(time = time_seq),
as.list(params[nome, ])))
u <- attr(u, "gradient")
colnames(u) <- sprintf("%s.nome%s",
colnames(u),
nome)
return(u)
})
# Cria matriz bloco diagonal e atribui nomes.
X <- as.matrix(Matrix::bdiag(grad))
colnames(X) <- c(sapply(grad, colnames))
# Troca colunas de lugar para corresponder com `vcov()`.
X <- X[, colnames(vcov(m1))]
# Cria erro padrão do valor predito e calcula o intervalo de confiança.
qtl <- qt(0.975, df = nrow(tb) - length(coef(m1)))
grid$se <- sqrt(diag(X %*% vcov(m1) %*% t(X)))
grid <- grid %>%
mutate(lwr = fit - qtl * se,
upr = fit + qtl * se)
tb_params <- params %>%
as.data.frame() %>%
rownames_to_column("nome")
tb_params <- grid %>%
filter(time == 40) %>%
select(nome, lwr, upr, fit) %>%
inner_join(tb_params) %>%
select(nome, A, B, C, fit, lwr, upr)cap <- "Observed (open circles) and fitted values (solid line) of soil water content relative to oven drying and drying period in microwave oven. The filled region arround the fitted curve is the 95% confidence band for the fitted values. Point estimates for $A$, $b$ and $c$ parameters are annotated in each soil plot. IC40 is the confidence interval (95%) for the predicted relative water content at 40 minutes. Dotted horizontal line is at 100%."
# Gráfico com bandas de confiança.
fmt <- "c: %0.2f\nb: %0.1f\nA: %0.1f\nIC40: [%0.1f; %0.1f]"
ggplot(tb,
aes(x = time, y = umrel)) +
facet_wrap(facets = ~nome, nrow = 2) +
geom_line(data = grid,
mapping = aes(x = time, y = fit)) +
geom_ribbon(data = grid,
inherit.aes = FALSE,
mapping = aes(x = time, ymin = lwr, ymax = upr),
fill = "darkorange",
alpha = 0.4) +
geom_point(pch = 1) +
geom_text(data = tb_params,
mapping = aes(x = 42,
y = 0,
label = sprintf(fmt,
C,
B,
A,
lwr,
upr)),
parse = FALSE,
size = 3.5,
hjust = 1,
vjust = 0) +
geom_hline(yintercept = 100,
lty = 3,
size = 0.5,
col = "darkorange") +
scale_y_continuous(breaks = seq(0, 100, by = 20)) +
scale_x_continuous(breaks = seq(0, 40, by = 10)) +
labs(x = "Drying period in microwave oven (min)",
y = "Relative water content (%)")
Observed (open circles) and fitted values (solid line) of soil water content relative to oven drying and drying period in microwave oven. The filled region arround the fitted curve is the 95% confidence band for the fitted values. Point estimates for \(A\), \(b\) and \(c\) parameters are annotated in each soil plot. IC40 is the confidence interval (95%) for the predicted relative water content at 40 minutes. Dotted horizontal line is at 100%.
tb_r2 <- tb %>%
group_by(nome) %>%
summarise(R2 = cor(umrel, fit)^2)
# Tabelas com estimativas pontuais para os parâmetros e valore predito
# para o time = 40 acompanhado do IC95%.
tb_params <- tb_params %>%
inner_join(tb_r2) %>%
rename_at(c("fit", "lwr", "upr"),
~paste(., "40", sep = ""))
# Tabela a ser incluida com as estimativas dos parâmetros de interesse.
tb_params %>%
mutate_if(is.numeric, round, digits = 3)## nome A B C fit40 lwr40 upr40 R2
## 1 CXa-A 107.271 13.000 0.081 95.954 90.241 101.668 0.955
## 2 CXbdf-A 106.109 6.284 0.101 102.620 97.563 107.677 0.947
## 3 CXbdf-C 108.416 11.220 0.082 98.692 93.094 104.290 0.958
## 4 LVAd-A 116.117 8.909 0.067 102.527 96.815 108.240 0.957
## 5 LVAd-Bw 109.921 11.623 0.090 101.780 96.269 107.290 0.956
## 6 LVdf-A 110.031 14.993 0.091 99.402 93.652 105.152 0.990
## 7 LVdf-B 102.205 13.780 0.101 95.278 89.732 100.823 0.970
## 8 LVu-Bw 99.162 2.516 0.156 98.879 94.774 102.983 0.947# Tabela das estimativas dos parâmetros com erros padrões.
summary(m1)$tTable## Value Std.Error t-value p-value
## A.nomeCXa-A 107.27076636 7.49991202 14.302937 6.811608e-29
## A.nomeCXbdf-A 106.10928886 4.09875455 25.888178 2.072395e-54
## A.nomeCXbdf-C 108.41646596 6.80189745 15.939150 6.453204e-33
## A.nomeLVAd-A 116.11745646 9.18597779 12.640729 1.045723e-24
## A.nomeLVAd-Bw 109.92149493 5.87022988 18.725245 1.799550e-39
## A.nomeLVdf-A 110.03124558 6.85612504 16.048605 3.506418e-33
## A.nomeLVdf-B 102.20532169 5.45495491 18.736236 1.698709e-39
## A.nomeLVu-Bw 99.16249669 2.26065418 43.864514 3.848443e-82
## C.nomeCXa-A 0.08125507 0.01254781 6.475638 1.590747e-09
## C.nomeCXbdf-A 0.10078658 0.01484876 6.787544 3.222896e-10
## C.nomeCXbdf-C 0.08216653 0.01255238 6.545892 1.113731e-09
## C.nomeLVAd-A 0.06701948 0.01260689 5.316101 4.234846e-07
## C.nomeLVAd-Bw 0.09037381 0.01233713 7.325352 1.901350e-11
## C.nomeLVdf-A 0.09144703 0.01255048 7.286337 2.341795e-11
## C.nomeLVdf-B 0.10131747 0.01364604 7.424678 1.116428e-11
## C.nomeLVu-Bw 0.15619137 0.02686428 5.814091 4.142388e-08
## B.nomeCXa-A 13.00003592 1.15862871 11.220191 4.348630e-21
## B.nomeCXbdf-A 6.28446129 0.69604509 9.028813 1.489910e-15
## B.nomeCXbdf-C 11.22032948 1.02750273 10.920000 2.532226e-20
## B.nomeLVAd-A 8.90940970 1.32184595 6.740127 4.117485e-10
## B.nomeLVAd-Bw 11.62315550 0.86352567 13.460116 8.830491e-27
## B.nomeLVdf-A 14.99345852 0.98496802 15.222279 3.618027e-31
## B.nomeLVdf-B 13.77968511 0.83413113 16.519807 2.576412e-34
## B.nomeLVu-Bw 2.51610722 0.79130578 3.179690 1.826021e-03#