Poder Calorífico Util em função do Teor de Umidade
Poder Calorífico Util em função do Teor de Umidade
“Elaborado por Profª Drª Isabel Homczinski”
“04 agosto, 2022”
1 Instalação dos pacotes
Instalação do pacote e verificação se todos estão instalados. O Comando abaixo ajuda a verificar se os pacotes estão instalados, ou se necessitam ser instalados.
pacotes <- c("tidyverse", "readxl", "ggpubr", "kableExtra",
"PerformanceAnalytics", "scatterplot3d", "MASS",
"gridExtra", "car", "lmtest", "plyr", "knitr", "dplyr")
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
instalador <- pacotes[!pacotes %in% installed.packages()]
for(i in 1:length(instalador)) {
install.packages(instalador, dependencies = T)
break()}
sapply(pacotes, require, character = T)
} else {
sapply(pacotes, require, character = T)
}## tidyverse readxl ggpubr
## TRUE TRUE TRUE
## kableExtra PerformanceAnalytics scatterplot3d
## TRUE TRUE TRUE
## MASS gridExtra car
## TRUE TRUE TRUE
## lmtest plyr knitr
## TRUE TRUE TRUE
## dplyr
## TRUE2 Poder Calorífico Util e Teor de Umidade
O Poder calorífico é uma medida da quantidade de energia térmica liberada, na forma de calor, na queima completa de 1kg de combustível. O objetivo principal dessa modelagem é verificar se o Poder Calorífico Util tem correlação com o Teor de Umidade. Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada. Adicionalmente, avaliar a qualidade dos dados experimentais por meio das seguintes estatísticas descritivas: média, desvio padrão, coeficiente de variação e score z. Como referência para a Área de Tecnologia da Madeira, nas disciplinas ministradas pela professora Gilmara, deve-se considerar dados de excelente qualidade se o Coeficiente de variação tiver valor de até 10%, de 10 a 20% dados bons, de 20 a 30 % dados ruins e acima de 30% dados péssimos. Também, no que se refere ao score Z, valores em módulo acima de 2 serão considerados como outliers (pontos atípicos). Os outliers também serão identificados por meio de gráficos box plot das variáveis do modelo.
2.1 Banco de Dados
tabela <- read_excel("Dados.xlsx",
sheet = "Teor_umidade")tabela %>%
kbl(caption = "Tabela dos Dados", align = 'cc') %>%
kable_classic(full_width = T, html_font = "Times", font_size = 12)| U | PCU |
|---|---|
| 0 | 4756 |
| 10 | 4221 |
| 20 | 3687 |
| 30 | 3153 |
| 40 | 2620 |
| 50 | 2085 |
| 60 | 1551 |
| 70 | 1016 |
2.2 Análise Descritiva
# Teor de Umidade
y <- tabela$U
mean_and_sd <- function(y) {
m <- mean(y) # média
s <- sd(y) # desvio padrão
# sprintf("%0.3f U00B1 %0.3f", m, s)
sprintf("%0.2f $\\pm$ %0.2f", m, s)}
Med_e_sd<- mean_and_sd(y)
cv <- function(y) {100 * sd(y)/mean(y)}
Coef_de_Variacao<-cv(y) # coeficiente de variação
Coef_de_Variacao<-round(Coef_de_Variacao, 2)
TU<-data.frame(Med_e_sd, Coef_de_Variacao)
# Poder Calorífico Util
y <- tabela$PCU
mean_and_sd <- function(y) {
m <- mean(y) # média
s <- sd(y) # desvio padrão
# sprintf("%0.3f U00B1 %0.3f", m, s)
sprintf("%0.2f $\\pm$ %0.2f", m, s)}
Med_e_sd<- mean_and_sd(y)
cv <- function(y) {100 * sd(y)/mean(y)}
Coef_de_Variacao<-cv(y) # coeficiente de variação
Coef_de_Variacao<-round(Coef_de_Variacao, 2)
PCU<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(TU,PCU)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Teor de Umidade (%)", "Poder Calorífico Util (kcal/kg)"))
dados<-data.frame(dados)
#Organizar as colunas no data.frame
date<-dplyr::select(dados, Variaveis, Med_e_sd, Coef_de_Variacao)
# Criar da Tabela Descritiva
knitr::kable(date, align = 'cc', caption = "Tabela Descritiva",
col.names = c("Variáveis", "Média $\\pm$ sd", " CV(%)")) %>%
kable_classic(full_width = F, html_font = "Times", font_size = 16)| Variáveis | Média \(\pm\) sd | CV(%) |
|---|---|---|
| Teor de Umidade (%) | 35.00 \(\pm\) 24.49 | 69.99 |
| Poder Calorífico Util (kcal/kg) | 2886.12 \(\pm\) 1308.41 | 45.33 |
2.2.1 Poder Calorífico Util - outlier
y <- tabela$PCU
scale(y) # padronização da variável ## [,1]
## [1,] 1.4291237
## [2,] 1.0202294
## [3,] 0.6120994
## [4,] 0.2039695
## [5,] -0.2033962
## [6,] -0.6122905
## [7,] -1.0204205
## [8,] -1.4293147
## attr(,"scaled:center")
## [1] 2886.125
## attr(,"scaled:scale")
## [1] 1308.407score_z <- function(y, media = NULL, despad = NULL) {if (is.null(media)) {
media <- mean(y)}
if (is.null(despad)) {despad <- sd(y)}
(y - media)/despad}
status <- ifelse(abs(score_z(y)) > 2, "outlier", "tipical") # ocorrência de outlier
table(status)## status
## tipical
## 82.2.2 Teor de Umidade - outlier
y <- tabela$U
scale(y) # padronização da variável ## [,1]
## [1,] -1.4288690
## [2,] -1.0206207
## [3,] -0.6123724
## [4,] -0.2041241
## [5,] 0.2041241
## [6,] 0.6123724
## [7,] 1.0206207
## [8,] 1.4288690
## attr(,"scaled:center")
## [1] 35
## attr(,"scaled:scale")
## [1] 24.4949score_z <- function(y, media = NULL, despad = NULL) {if (is.null(media)) {
media <- mean(y)}
if (is.null(despad)) {despad <- sd(y)}
(y - media)/despad}
status <- ifelse(abs(score_z(y)) > 2, "outlier", "tipical") # ocorrência de outlier
table(status)## status
## tipical
## 82.3 Gráfico Box Plot
par(mfrow=c(1,2))
boxplot(tabela$U,
col="red",
main="Box Plot", #título,
ylab= "Teor de Umidade (%)", # texto do eixo y,
adj=0.5, #alinhamento dos textos,
ylim= c(0,70)) # limite do eixo y,)
boxplot(tabela$PCU, col="blue", main="Box Plot", ylab= "Poder Calorífico Util (kcal/kg)",
adj=0.5, ylim= c(1000,5000))
2.4 Dispersão, Histograma e Correlação
chart.Correlation(tabela, histogram =TRUE)
2.4.1 Regressão Linear
# Ajuste do modelo
modelo<-lm(formula=PCU~U, data=tabela); modelo##
## Call:
## lm(formula = PCU ~ U, data = tabela)
##
## Coefficients:
## (Intercept) U
## 4755.67 -53.42# Estatísticas das estimativas dos paramentros do modelo
summary(modelo)##
## Call:
## lm(formula = PCU ~ U, data = tabela)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.58333 -0.39583 -0.04762 0.27976 0.95238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.756e+03 3.603e-01 13200 <2e-16 ***
## U -5.342e+01 8.612e-03 -6202 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5581 on 6 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 3.847e+07 on 1 and 6 DF, p-value: < 2.2e-16# Anova da regressão do modelo
anova(modelo)## Analysis of Variance Table
##
## Response: PCU
## Df Sum Sq Mean Sq F value Pr(>F)
## U 1 11983495 11983495 38469309 < 2.2e-16 ***
## Residuals 6 2 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.4.2 Medidas de Precisão para Análise do Modelo
PCU_est<-predict.lm(modelo)
newdata1 = data.frame(tabela, PCU_est)
Res<-(((newdata1$PCU-newdata1$PCU_est)/(newdata1$PCU))*100)
newdata1 = data.frame(tabela, PCU_est,Res)
resid<-(newdata1$PCU-newdata1$PCU_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
SQTotal<-var(newdata1$PCU)*((length(newdata1$PCU))-1)
R2<-(1-(SQResiduo/SQTotal))
a<-(length(newdata1$PCU)-1)/(length(newdata1$PCU)-2)
R2adj<-(1-(a*(1-R2)))
Syx<-sqrt(SQResiduo/((length(newdata1$PCU))-2))
Syx_Porc<-(Syx/(mean(newdata1$PCU)))*100
R2<-round(R2, 2)
R2adj<-round(R2adj,2)
Syx<-round(Syx,2)
Syx_Porc<-round(Syx_Porc, 2)
newdata11 = data.frame(R2,R2adj,Syx,Syx_Porc)
# Criar da Tabela Descritiva
knitr::kable(newdata11, align = 'cc',
caption = "Medidas de Precisão do Modelo",
col.names = c("R$^{2}$", "R$_{adj}^{2}$", "S$_{yx}$", "S$_{yx}$%")) %>%
kable_classic(full_width = F, html_font = "Times", font_size = 16)| R\(^{2}\) | R\(_{adj}^{2}\) | S\(_{yx}\) | S\(_{yx}\)% |
|---|---|---|---|
| 1 | 1 | 0.56 | 0.02 |
2.4.3 Predição do PCU aplicando o modelo
newdatas<-data.frame(U= c(50))
PCU_estimado1 = predict(modelo, newdatas, interval = "confidence");
PCU_estimado1 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 2084.893 2084.316 2085.472.4.4 Testes de Hipóteses
NORMALIDADE:
Ho: resíduos seguem distribuição normal;
Hi: resíduos não seguem distribuição normal
Resultado:
Valor de p > 0.05, não rejeita Ho para alfa=5%.
## Normalidade dos Resíduos
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.93568, p-value = 0.5692ks.test(modelo$residuals, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelo$residuals
## D = 0.27983, p-value = 0.4755
## alternative hypothesis: two-sidedINDEPENDÊNCIA DOS RESÍDUOS:
Ho: resíduos são independentes;
Hi: resíduos não são independentes
Resultado:
Valor de p > 0.05, não rejeita Ho para alfa=5%.
## Independencia dos Resíduos (Durbin-Watson)
durbinWatsonTest(modelo)## lag Autocorrelation D-W Statistic p-value
## 1 -0.07006369 1.89862 0.476
## Alternative hypothesis: rho != 0HOMOSCEDASTICIDADE DOS RESÍDUOS:
Ho: resíduos são homogênios;
Hi: resíduos não são homogênios;
Resultado:
Valor de p > 0.05, não rejeita Ho para alfa=5%.
## Homocedasticidade (Breusch-Pagan)
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 0.10485, df = 1, p-value = 0.7461ncvTest(modelo)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.07270845, Df = 1, p = 0.787432.4.5 Gráfico 1:1 e de Resíduos
Valores em vermelho estimados; valores em azul observados.
G1<-ggplot(newdata1) +
geom_point(size=2, aes(x=PCU, y=PCU_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCU Observado (kcal/kg)",
y = "PCU Estimado (kcal/kg)",
title = "Gráfico 1:1")+ # legendas
scale_x_continuous(breaks = seq(1000, 5000, 1000), limits= c(1000, 5000)) +
scale_y_continuous(breaks = seq(1000, 5000, 1000), limits= c(1000, 5000)) +
theme( plot.title = element_text(hjust = 0.5),
legend.title = element_text(hjust = 0.5))+ # centralizar títulos
theme(plot.title = element_text(size = 13, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)) # negrito nos títulos
G2<-ggplot(newdata1) +
geom_point(size=2, aes(x=PCU, y=Res)
,color="black")+ # gráfico de pontos
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCU Observado (kcal/kg)",
y = "Resíduo (%)",
title = "Gráfico de Resíduos")+ # legendas
scale_x_continuous(breaks = seq(1000, 5000, 1000), limits= c(1000, 5000))+
scale_y_continuous(breaks = seq(-30, 30, 10), limits= c(-30,30)) +
theme( plot.title = element_text(hjust = 0.5),
legend.title = element_text(hjust = 0.5))+ # centralizar títulos
theme(plot.title = element_text(size = 13, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10))+
# negrito nos títulos
geom_hline(yintercept = 0, color="black")+ # linha no 0 do eixo y
theme(axis.title.y = element_text(angle = 90))
grid.arrange(G1,G2, ncol=1) 
2.5 Regressão Linear
grafico<-ggplot(data = tabela, #banco de dados
mapping = aes(x = U, y = PCU)) + # eixo x e y
geom_point(size=2) + # gráfico de dispersão
geom_smooth(method = "lm", # método linear
formula='y ~ x', # fórmula
color="red", # cor
lwd=0.8)+ # tamanho da linha
theme_classic()+ # modo classico de gráfico (fundo branco)
stat_regline_equation(aes(label=paste(..eq.label.., # inserir a equação
..rr.label.., # inserir o R2
sep = "~~~~")), # espaço entre a formula e o R2
formula=y ~ x, # formula
label.x = 0, label.y=5000)+ # localização xy da formula
scale_x_continuous(breaks = seq(0, 70, 10), limits= c(0, 70)) +
scale_y_continuous(breaks = seq(1000, 5000, 1000), limits= c(1000, 5000))+
labs(x = "Teor de Umidade (%)", # legenda do eixo x
y = "Poder Calorífico Util (kcal/kg)")+ # legenda do eixo y
theme(plot.title = element_text(hjust = 0.5),
legend.title = element_text(hjust = 0.5))+ # Alinhamento do Título e eixos
theme(plot.title = element_text(size = 12, face = "bold",family="TT Times New Roman"),
axis.title = element_text(size = 12, face = "bold", family="TT Times New Roman"),
axis.text = element_text(size = 12, face = "bold", family="TT Times New Roman"))
# formatação dos títulos e eixos, com tamanho da fonte (size), negrido (bold) e
# tipo de fonte (TI Times New Roman)
grafico
3 Citações
De acordo com GROLEMUND (2018), a ciência de dados é uma área em crescimento.
A regressão não linear geralmente possui parâmetros interpretáveis (ROSSE; VENCOVSKY, 2000). Para CHEN et al. (2021), a análise multivariada é a melhor opção.
Referências
CHEN, Xue; LI, Xiaohui; XIE, Jinmei; YANG, Hao; LIU, Aichun.
Non-invasive discrimination of multiple myeloma using label-free serum
surface-enhanced raman scattering spectroscopy in combination with
multivariate analysis. Analytica Chimica Acta, p.
339296, Nov. 2021. DOI 10.1016/j.aca.2021.339296.
Available at: https://doi.org/10.1016/j.aca.2021.339296.
GROLEMUND, Garrett. R para data science. S.l: Alta
Books, 2018.
ROSSE, Leornado Novaes; VENCOVSKY, Roland. Modelo de
regressão não-linear aplicado ao estudo da
estabilidade fenotı́pica de genótipos de
feijão no estado do paraná.
Bragantia, vol. 59, no. 1, p. 99–107, 2000. DOI 10.1590/s0006-87052000000100016.
Available at: https://doi.org/10.1590/s0006-87052000000100016.