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 Análise Gráfica do Modelo
par(mfrow=c(2,2))
plot(modelo, pch=20)
2.4.5 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.494
## 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.6 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 Densidade e Propriedades Mecânicas do Combustível
A densidade se refere a quantidade de combustível por unidade de volume. O objetivo principal dessa modelagem é verificar se a Densidade tem correlação com as propriedades mecânicas do combustível (resistência e rigidez). Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência de cada variável no valor da densidade, por meio da comparação da magnitude dos betas da equação de regressão ajustada. 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.
3.1 Banco de Dados
tabela1 <- read_excel("Dados.xlsx",
sheet = "Densidade")tabela1 %>%
kbl(caption = "Tabela dos Dados", align = 'cc') %>%
kable_classic(full_width = F, html_font = "Times", font_size = 12)%>%
scroll_box(width = "400px", height = "300px")| Densidade | Resistencia | Rigidez |
|---|---|---|
| 688 | 50.5 | 12876 |
| 1170 | 79.5 | 20827 |
| 694 | 59.8 | 12912 |
| 1170 | 76.7 | 16694 |
| 803 | 48.1 | 13481 |
| 677 | 59.1 | 14098 |
| 871 | 52.0 | 14613 |
| 801 | 56.0 | 16224 |
| 759 | 54.8 | 11105 |
| 504 | 39.0 | 9839 |
| 500 | 31.5 | 8058 |
| 1090 | 93.2 | 23002 |
| 838 | 54.4 | 13627 |
| 1221 | 83.8 | 19426 |
| 705 | 47.3 | 13409 |
| 899 | 48.0 | 13286 |
| 999 | 62.0 | 18421 |
| 822 | 51.8 | 13963 |
| 690 | 48.9 | 18029 |
| 640 | 40.3 | 12813 |
| 931 | 63.5 | 18099 |
| 924 | 48.3 | 14431 |
| 929 | 54.9 | 16782 |
| 1087 | 72.7 | 19881 |
| 952 | 51.6 | 15561 |
| 948 | 78.5 | 19360 |
| 731 | 46.8 | 14933 |
| 899 | 57.7 | 17198 |
| 755 | 53.9 | 14617 |
| 889 | 42.7 | 14577 |
| 739 | 46.0 | 13166 |
| 892 | 78.4 | 18359 |
| 825 | 71.4 | 14624 |
| 919 | 62.4 | 17212 |
| 1068 | 76.0 | 18011 |
| 1074 | 93.3 | 23607 |
| 684 | 56.5 | 14185 |
| 1143 | 82.9 | 22733 |
| 856 | 71.4 | 18971 |
| 756 | 69.9 | 14719 |
| 544 | 37.8 | 9067 |
| 1106 | 95.2 | 21724 |
| 940 | 79.5 | 19583 |
| 580 | 40.9 | 15225 |
| 579 | 35.4 | 8431 |
| 537 | 32.6 | 7110 |
| 535 | 42.3 | 9868 |
| 560 | 40.4 | 11889 |
| 538 | 43.6 | 10904 |
| 645 | 44.4 | 13304 |
3.2 Análise Descritiva
# Densidade
y <- tabela1$Densidade
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.0f $\\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)
Densidade<-data.frame(Med_e_sd, Coef_de_Variacao)
# Resistencia
y <- tabela1$Resistencia
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.0f $\\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)
Resistencia<-data.frame(Med_e_sd, Coef_de_Variacao)
#Rigidez
y <- tabela1$Rigidez
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.0f $\\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)
Rigidez<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(Densidade,Resistencia)
dados<-rbind(dados, Rigidez)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Densidade (g/cm$^{3}$)", "Resistência (MPa)", "Rigidez (MPa)"))
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(%) |
|---|---|---|
| Densidade (g/cm\(^{3}\)) | 822 \(\pm\) 197.37 | 24.01 |
| Resistência (MPa) | 58 \(\pm\) 16.77 | 28.83 |
| Rigidez (MPa) | 15297 \(\pm\) 3930.17 | 25.69 |
3.2.1 Densidade - outlier
y <- tabela1$Densidade
scale(y) # padronização da variável ## [,1]
## [1,] -0.6795386825
## [2,] 1.7625851242
## [3,] -0.6491388010
## [4,] 1.7625851242
## [5,] -0.0968742888
## [6,] -0.7352717984
## [7,] 0.2476577006
## [8,] -0.1070075826
## [9,] -0.3198067524
## [10,] -1.6118017124
## [11,] -1.6320683000
## [12,] 1.3572533720
## [13,] 0.0804583528
## [14,] 2.0209841161
## [15,] -0.5934056851
## [16,] 0.3895238138
## [17,] 0.8961885040
## [18,] -0.0006079976
## [19,] -0.6694053886
## [20,] -0.9227377337
## [21,] 0.5516565147
## [22,] 0.5161899863
## [23,] 0.5415232208
## [24,] 1.3420534313
## [25,] 0.6580560996
## [26,] 0.6377895120
## [27,] -0.4616728657
## [28,] 0.3895238138
## [29,] -0.3400733400
## [30,] 0.3388573448
## [31,] -0.4211396905
## [32,] 0.3540572855
## [33,] 0.0145919431
## [34,] 0.4908567518
## [35,] 1.2457871402
## [36,] 1.2761870216
## [37,] -0.6998052701
## [38,] 1.6257856578
## [39,] 0.1716579970
## [40,] -0.3350066931
## [41,] -1.4091358363
## [42,] 1.4383197224
## [43,] 0.5972563368
## [44,] -1.2267365478
## [45,] -1.2318031947
## [46,] -1.4446023646
## [47,] -1.4547356584
## [48,] -1.3280694859
## [49,] -1.4395357177
## [50,] -0.8974044992
## attr(,"scaled:center")
## [1] 822.12
## attr(,"scaled:scale")
## [1] 197.3692score_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
## outlier tipical
## 1 493.2.2 Resistência - outlier
y <- tabela1$Resistencia
scale(y) # padronização da variável ## [,1]
## [1,] -0.45634808
## [2,] 1.27314673
## [3,] 0.09828302
## [4,] 1.10616103
## [5,] -0.59947868
## [6,] 0.05653659
## [7,] -0.36689145
## [8,] -0.12834044
## [9,] -0.19990575
## [10,] -1.14218223
## [11,] -1.58946537
## [12,] 2.09018394
## [13,] -0.22376085
## [14,] 1.52958907
## [15,] -0.64718889
## [16,] -0.60544246
## [17,] 0.22948607
## [18,] -0.37881900
## [19,] -0.55176848
## [20,] -1.06465315
## [21,] 0.31894270
## [22,] -0.58755113
## [23,] -0.19394197
## [24,] 0.86761002
## [25,] -0.39074655
## [26,] 1.21350898
## [27,] -0.67700776
## [28,] -0.02695626
## [29,] -0.25357972
## [30,] -0.92152255
## [31,] -0.72471796
## [32,] 1.20754521
## [33,] 0.79008094
## [34,] 0.25334117
## [35,] 1.06441460
## [36,] 2.09614771
## [37,] -0.09852157
## [38,] 1.47591509
## [39,] 0.79008094
## [40,] 0.70062431
## [41,] -1.21374753
## [42,] 2.20945944
## [43,] 1.27314673
## [44,] -1.02887050
## [45,] -1.35687814
## [46,] -1.52386384
## [47,] -0.94537765
## [48,] -1.05868938
## [49,] -0.86784857
## [50,] -0.82013837
## attr(,"scaled:center")
## [1] 58.152
## attr(,"scaled:scale")
## [1] 16.7679score_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
## outlier tipical
## 3 473.2.3 Rigidez - outlier
y <- tabela1$Rigidez
scale(y) # padronização da variável ## [,1]
## [1,] -0.61592270
## [2,] 1.40714577
## [3,] -0.60676279
## [4,] 0.35553692
## [5,] -0.46198528
## [6,] -0.30499456
## [7,] -0.17395692
## [8,] 0.23594917
## [9,] -1.06653951
## [10,] -1.38866311
## [11,] -1.84182434
## [12,] 1.96055716
## [13,] -0.42483675
## [14,] 1.05067251
## [15,] -0.48030511
## [16,] -0.51160147
## [17,] 0.79495828
## [18,] -0.33934423
## [19,] 0.69521701
## [20,] -0.63195255
## [21,] 0.71302795
## [22,] -0.22026537
## [23,] 0.37792782
## [24,] 1.16644363
## [25,] 0.06725411
## [26,] 1.03387934
## [27,] -0.09253547
## [28,] 0.48377570
## [29,] -0.17293915
## [30,] -0.18311683
## [31,] -0.54213452
## [32,] 0.77918288
## [33,] -0.17115806
## [34,] 0.48733788
## [35,] 0.69063705
## [36,] 2.11449458
## [37,] -0.28285810
## [38,] 1.89211226
## [39,] 0.93490139
## [40,] -0.14698606
## [41,] -1.58509234
## [42,] 1.63538026
## [43,] 1.09061990
## [44,] -0.01823840
## [45,] -1.74691747
## [46,] -2.08303537
## [47,] -1.38128429
## [48,] -0.86705697
## [49,] -1.11768236
## [50,] -0.50702152
## attr(,"scaled:center")
## [1] 15296.68
## attr(,"scaled:scale")
## [1] 3930.169score_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
## outlier tipical
## 2 483.3 Gráfico Box Plot
par(mfrow=c(1,3))
boxplot(tabela1$Densidade, col="blue", main="Box Plot",
ylab= "Densidade (g/cm³)", adj=0.5, ylim= c(400,1200))
boxplot(tabela1$Resistencia, col="purple",main="Box Plot",
ylab= "Resistência (MPa)", adj=0.5, ylim= c(30,100))
boxplot(tabela1$Rigidez,col="green2",main="Box Plot", ylab= "Rigidez (MPa)",
adj=0.5, ylim= c(1000,30000)) 
3.4 Dispersão, Histograma e Correlação
chart.Correlation(tabela1, histogram =TRUE)
3.5 Regressão Linear Múltipla
# Ajuste do modelo
modelo1<-lm(formula =Densidade~Resistencia+Rigidez,
data=tabela1); modelo1##
## Call:
## lm(formula = Densidade ~ Resistencia + Rigidez, data = tabela1)
##
## Coefficients:
## (Intercept) Resistencia Rigidez
## 164.91421 4.26645 0.02674# Estatísticas das estimativas dos paramentros do modelo
summary(modelo1)##
## Call:
## lm(formula = Densidade ~ Resistencia + Rigidez, data = tabela1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -166.60 -84.37 -22.51 75.49 231.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.649e+02 5.712e+01 2.887 0.00586 **
## Resistencia 4.266e+00 1.809e+00 2.359 0.02253 *
## Rigidez 2.675e-02 7.716e-03 3.466 0.00114 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 99.55 on 47 degrees of freedom
## Multiple R-squared: 0.756, Adjusted R-squared: 0.7456
## F-statistic: 72.8 on 2 and 47 DF, p-value: 4.022e-15# Anova da regressão do modelo
anova(modelo1)## Analysis of Variance Table
##
## Response: Densidade
## Df Sum Sq Mean Sq F value Pr(>F)
## Resistencia 1 1323934 1323934 133.593 2.442e-15 ***
## Rigidez 1 119063 119063 12.014 0.001139 **
## Residuals 47 465778 9910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.5.1 Medidas de Precisão para Análise do Modelo
Densidade_est<-predict.lm(modelo1)
newdata1 = data.frame(tabela1, Densidade_est)
Res<-(((newdata1$Densidade-newdata1$Densidade_est)/(newdata1$Densidade))*100)
newdata1 = data.frame(tabela1, Densidade_est,Res)
resid<-(newdata1$Densidade-newdata1$Densidade_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
SQTotal<-var(newdata1$Densidade)*((length(newdata1$Densidade))-1)
R2<-(1-(SQResiduo/SQTotal))
a<-(length(newdata1$Densidade)-1)/(length(newdata1$Densidade)-2)
R2adj<-(1-(a*(1-R2)))
Syx<-sqrt(SQResiduo/((length(newdata1$Densidade))-2))
Syx_Porc<-(Syx/(mean(newdata1$Densidade)))*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}\)% |
|---|---|---|---|
| 0.76 | 0.75 | 98.51 | 11.98 |
3.5.2 Predição da Densidade aplicando o modelo
newdatas<-data.frame(Resistencia= c(50),Rigidez= c(12876))
Densidade_estimado1 = predict(modelo1, newdatas, interval = "confidence");
Densidade_estimado1 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 722.5998 689.0572 756.14253.5.3 Análise grafica do modelo
par(mfrow=c(2,2))
plot(modelo1, pch=20)
3.5.4 Testes de hipoteses
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(modelo1$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo1$residuals
## W = 0.96394, p-value = 0.13ks.test(modelo1$residuals, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelo1$residuals
## D = 0.54, p-value = 3.642e-14
## 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(modelo1)## lag Autocorrelation D-W Statistic p-value
## 1 0.1013323 1.785325 0.426
## 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(modelo1)##
## studentized Breusch-Pagan test
##
## data: modelo1
## BP = 1.6976, df = 2, p-value = 0.4279ncvTest(modelo1)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.9172137, Df = 1, p = 0.338213.5.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=Densidade, y=Densidade_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
geom_point(size=2, aes(x=Densidade_est, y=Densidade), color="blue2")+
theme_classic()+ # aparencia clean do gráfico
labs(x = "Densidade Observada (g/cm³)",
y = "Densidade Estimada (g/cm³)",
title = "Gráfico 1:1")+ # legendas
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=Densidade, y=Res)
,color="black")+ # gráfico de pontos
theme_classic()+ # aparencia clean do gráfico
labs(x = "Densidade Observada (g/cm³)",
y = "Resíduo (%)",
title = "Gráfico de Resíduos")+ # legendas
scale_x_continuous(breaks = seq(450, 1250, 250), limits= c(450,1250)) + # eixo x
scale_y_continuous(breaks = seq(-30, 30, 10), limits= c(-30,30)) + # eixo y
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) 
3.5.6 Gráfico 3D de Dispersão em função do modelo
graph<-scatterplot3d(tabela1$Densidade~tabela1$Resistencia+tabela1$Rigidez,
pch=16, angle=45,color="blue", box=FALSE,
xlab="Resistência (MPa)", ylab="Rigidez (Mpa)", zlab = "Densidade (g/cm³)")
# previsão do modelo: Quando o modelo é apropriado todos os pontos ficam dentro
#do plano.
graph$plane3d(modelo1, col="black", draw_polygon=TRUE) 
3.5.7 Dados Padronizados
tabela_padronizada1<- scale(tabela1) # padrozinação do banco de dados
tabela_padronizada1<- data.frame(tabela_padronizada1) # tornar os dados padronizados em `data.frame`
lm(formula =tabela_padronizada1$Densidade~tabela_padronizada1$Resistencia+
tabela_padronizada1$Rigidez,
data=tabela_padronizada1) # coeficientes com dados padronizados##
## Call:
## lm(formula = tabela_padronizada1$Densidade ~ tabela_padronizada1$Resistencia +
## tabela_padronizada1$Rigidez, data = tabela_padronizada1)
##
## Coefficients:
## (Intercept) tabela_padronizada1$Resistencia
## -2.002e-17 3.625e-01
## tabela_padronizada1$Rigidez
## 5.326e-013.5.8 Método de Seleção de modelos utilizando critérios matemáticos.
# O modelo 1 é o que nós criamos com todas as variaveis, o mod.simples é o modelo contendo apenas o intercepto. E com a função StepAIC o R roda qual seria o melhor ajuste e quais variáveis são relevantes para o modelo.
modelo1<-lm(formula =tabela1$Densidade~tabela1$Resistencia+tabela1$Rigidez, data=tabela1)
mod.simples<-lm(formula=tabela1$Densidade~1, data=tabela1)
stepAIC(modelo1, scope=list(upper=modelo1,
lower=mod.simples), direction = "backward")## Start: AIC=462.97
## tabela1$Densidade ~ tabela1$Resistencia + tabela1$Rigidez
##
## Df Sum of Sq RSS AIC
## <none> 465778 462.97
## - tabela1$Resistencia 1 55154 520932 466.57
## - tabela1$Rigidez 1 119063 584841 472.35##
## Call:
## lm(formula = tabela1$Densidade ~ tabela1$Resistencia + tabela1$Rigidez,
## data = tabela1)
##
## Coefficients:
## (Intercept) tabela1$Resistencia tabela1$Rigidez
## 164.91421 4.26645 0.026744 Poder Calorífico Superior e a Composição Elementar
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 (PC) tem correlação com a Composição Elementar da biomassa (teor de carbono, hidrogênio e oxigênio). Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência de cada variável no valor do PC, por meio da comparação da magnitude dos betas da equação de regressão ajustada. 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.
4.1 Banco de Dados
tabela2 <- read_excel("Dados.xlsx",
sheet = "Analise_elementar")tabela2 %>%
kbl(caption = "Tabela dos Dados", align = 'cc') %>%
kable_classic(full_width = F, html_font = "Times", font_size = 12)| C | H | O | PCS |
|---|---|---|---|
| 83.67 | 3.56 | 2.84 | 32.856 |
| 82.62 | 3.02 | 3.66 | 33.000 |
| 63.89 | 4.97 | 24.54 | 25.100 |
| 92.04 | 2.45 | 2.96 | 34.388 |
| 89.13 | 0.43 | 0.98 | 31.124 |
| 46.90 | 6.07 | 43.99 | 18.261 |
| 49.14 | 6.34 | 43.52 | 19.423 |
| 54.41 | 4.99 | 39.69 | 21.010 |
| 48.79 | 5.91 | 43.41 | 19.260 |
| 46.58 | 5.87 | 45.46 | 18.770 |
| 48.10 | 5.99 | 45.74 | 19.916 |
| 47.84 | 5.80 | 45.76 | 18.981 |
| 50.64 | 5.98 | 42.88 | 20.720 |
| 48.15 | 5.87 | 44.75 | 19.777 |
| 46.04 | 5.82 | 44.49 | 18.640 |
4.2 Análise Descritiva
# Poder Calorífico Superior
y <- tabela2$PCS
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)
PCS<-data.frame(Med_e_sd, Coef_de_Variacao)
# Carbono
y <- tabela2$C
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)
Cb<-data.frame(Med_e_sd, Coef_de_Variacao)
#Hidrogenio
y <- tabela2$H
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)
Hd<-data.frame(Med_e_sd, Coef_de_Variacao)
#Oxigenio
y <- tabela2$O
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)
Ox<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(PCS,Cb)
dados<-rbind(dados, Hd)
dados<-rbind(dados, Ox)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Poder Calorífico Superior (MJ/kg)",
"Carbono (%)", "Hidrogênio (%)", "Oxigênio (%)"))
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(%) |
|---|---|---|
| Poder Calorífico Superior (MJ/kg) | 23.42 \(\pm\) 6.13 | 26.18 |
| Carbono (%) | 59.86 \(\pm\) 17.53 | 29.28 |
| Hidrogênio (%) | 4.87 \(\pm\) 1.73 | 35.42 |
| Oxigênio (%) | 31.64 \(\pm\) 18.85 | 59.57 |
4.2.1 Poder Calorífico Superior - outlier
y <- tabela2$PCS
scale(y) # padronização da variável ## [,1]
## [1,] 1.5402241
## [2,] 1.5637167
## [3,] 0.2748854
## [4,] 1.7901595
## [5,] 1.2576601
## [6,] -0.8408509
## [7,] -0.6512785
## [8,] -0.3923703
## [9,] -0.6778709
## [10,] -0.7578110
## [11,] -0.5708489
## [12,] -0.7233878
## [13,] -0.4396818
## [14,] -0.5935258
## [15,] -0.7790197
## attr(,"scaled:center")
## [1] 23.41507
## attr(,"scaled:scale")
## [1] 6.129584score_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
## 154.2.2 Carbono - outlier
y <- tabela2$C
scale(y) # padronização da variável ## [,1]
## [1,] 1.3583062
## [2,] 1.2983994
## [3,] 0.2297759
## [4,] 1.8358492
## [5,] 1.6698217
## [6,] -0.7395734
## [7,] -0.6117722
## [8,] -0.3110971
## [9,] -0.6317412
## [10,] -0.7578307
## [11,] -0.6711085
## [12,] -0.6859426
## [13,] -0.5261911
## [14,] -0.6682558
## [15,] -0.7886400
## attr(,"scaled:center")
## [1] 59.86267
## attr(,"scaled:scale")
## [1] 17.52722score_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
## 154.2.3 Hidrogênio - outlier
y <- tabela2$H
scale(y) # padronização da variável ## [,1]
## [1,] -0.76009633
## [2,] -1.07309990
## [3,] 0.05719078
## [4,] -1.40349256
## [5,] -2.57435777
## [6,] 0.69479064
## [7,] 0.85129243
## [8,] 0.06878350
## [9,] 0.60204884
## [10,] 0.57886340
## [11,] 0.64841974
## [12,] 0.53828886
## [13,] 0.64262338
## [14,] 0.57886340
## [15,] 0.54988158
## attr(,"scaled:center")
## [1] 4.871333
## attr(,"scaled:scale")
## [1] 1.72522score_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
## outlier tipical
## 1 144.2.4 Oxigênio - outlier
y <- tabela2$O
scale(y) # padronização da variável ## [,1]
## [1,] -1.5280538
## [2,] -1.4845538
## [3,] -0.3768942
## [4,] -1.5216880
## [5,] -1.6267246
## [6,] 0.6549055
## [7,] 0.6299725
## [8,] 0.4267955
## [9,] 0.6241371
## [10,] 0.7328872
## [11,] 0.7477409
## [12,] 0.7488019
## [13,] 0.5960213
## [14,] 0.6952226
## [15,] 0.6814299
## attr(,"scaled:center")
## [1] 31.64467
## attr(,"scaled:scale")
## [1] 18.85056score_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
## 154.3 Gráfico Box Plot
par(mfrow=c(1,4))
boxplot(tabela2$C, col="blue", main="Box Plot",
ylab= "Carbono (%)", adj=0.5, ylim=c(40,100))
boxplot(tabela2$H, col="purple",main="Box Plot",
ylab= "Hidrogênio (%)", adj=0.5, ylim=c(2,7))
boxplot(tabela2$O,col="green2",main="Box Plot", ylab= "Oxigênio (%)",
adj=0.5, ylim=c(0,50))
boxplot(tabela2$PCS,col="gray",main="Box Plot", ylab= "PCS (MJ/kg)",
adj=0.5, ylim=c(15,35)) 
4.4 Dispersão, Histograma e Correlação
chart.Correlation(tabela2, histogram =TRUE)
4.5 Regressão Linear Múltipla
# Ajuste do modelo
modelo2<-lm(formula =PCS~-1+C+H+O, data=tabela2); modelo2##
## Call:
## lm(formula = PCS ~ -1 + C + H + O, data = tabela2)
##
## Coefficients:
## C H O
## 0.34876 1.14845 -0.09662# Estatísticas das estimativas dos paramentros do modelo
summary(modelo2)##
## Call:
## lm(formula = PCS ~ -1 + C + H + O, data = tabela2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81639 -0.35435 0.05688 0.26616 1.07100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## C 0.348757 0.006411 54.401 9.8e-16 ***
## H 1.148451 0.224550 5.114 0.000256 ***
## O -0.096622 0.024853 -3.888 0.002157 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5792 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 8690 on 3 and 12 DF, p-value: < 2.2e-16# Anova da regressão do modelo
anova(modelo2)## Analysis of Variance Table
##
## Response: PCS
## Df Sum Sq Mean Sq F value Pr(>F)
## C 1 8729.7 8729.7 26022.698 < 2.2e-16 ***
## H 1 11.2 11.2 33.345 8.818e-05 ***
## O 1 5.1 5.1 15.115 0.002157 **
## Residuals 12 4.0 0.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.5.1 Predição do PCS aplicando o modelo
newdatas<-data.frame(C= c(83),H= c(3), O= c(3))
PCS_estimado2 = predict(modelo2, newdatas, interval = "confidence");
PCS_estimado2 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 32.10232 31.43632 32.768334.5.2 Medidas de Precisão para Análise do Modelo
PCS2_est<-predict.lm(modelo2)
newdata1 = data.frame(tabela2, PCS2_est)
Res<-(((newdata1$PCS-newdata1$PCS2_est)/(newdata1$PCS))*100)
newdata1 = data.frame(tabela2, PCS2_est,Res)
resid<-(newdata1$PCS-newdata1$PCS2_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
SQTotal<-var(newdata1$PCS)*((length(newdata1$PCS))-1)
R2<-(1-(SQResiduo/SQTotal))
a<-(length(newdata1$PCS)-1)/(length(newdata1$PCS)-2)
R2adj<-(1-(a*(1-R2)))
Syx<-sqrt(SQResiduo/((length(newdata1$PCS))-2))
Syx_Porc<-(Syx/(mean(newdata1$PCS)))*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}\)% |
|---|---|---|---|
| 0.99 | 0.99 | 0.56 | 2.38 |
4.5.3 Análise grafica do modelo
par(mfrow=c(2,2))
plot(modelo2, pch=20)
4.5.4 Testes de hipoteses
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(modelo2$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.97638, p-value = 0.9387ks.test(modelo2$residuals, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelo2$residuals
## D = 0.20714, p-value = 0.4781
## 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(modelo2)## lag Autocorrelation D-W Statistic p-value
## 1 0.1754137 1.634667 0.29
## 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(modelo2)##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 1.6078, df = 2, p-value = 0.4476ncvTest(modelo2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.1071283, Df = 1, p = 0.743444.5.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=PCS, y=PCS2_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
geom_point(size=2, aes(x=PCS2_est, y=PCS), color="blue2")+
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observada (MJ/kg)",
y = "PCS Estimada (Mj/kg)",
title = "Gráfico 1:1")+ # legendas
scale_x_continuous(breaks = seq(10, 40, 5), limits= c(10,40))+
scale_y_continuous(breaks = seq(10, 40, 5), limits= c(10,40))+
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=PCS, y=Res)
,color="black")+ # gráfico de pontos
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observada (g/cm³)",
y = "Resíduo (%)",
title = "Gráfico de Resíduos")+ # legendas
scale_x_continuous(breaks = seq(10, 40, 5), limits= c(10,40)) +
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) 
4.5.6 Dados Padronizados
tabela_padronizada2<- scale(tabela2) # padrozinação do banco de dados
tabela_padronizada2<- data.frame(tabela_padronizada2) # tornar os dados padronizados em `data.frame`
lm(formula =PCS~-1+C+H+O, data=tabela_padronizada2) # coeficientes com dados padronizados##
## Call:
## lm(formula = PCS ~ -1 + C + H + O, data = tabela_padronizada2)
##
## Coefficients:
## C H O
## 0.9539 0.3158 -0.33354.5.7 Método de Seleção de modelos utilizando critérios matemáticos.
# O modelo 2 é o que nós criamos com todas as variaveis, o mod.simples é o modelo contendo apenas o intercepto. E com a função StepAIC o R roda qual seria o melhor ajuste e quais variáveis são relevantes para o modelo.
modelo2<-lm(formula =PCS~-1+C+H+O, data=tabela2)
mod.simples2<-lm(formula=PCS~1, data=tabela2)
stepAIC(modelo2, scope=list(upper=modelo2,
lower=mod.simples2), direction = "backward")## Start: AIC=-13.73
## PCS ~ -1 + C + H + O
##
## Df Sum of Sq RSS AIC
## <none> 4.03 -13.731
## - O 1 5.07 9.10 -3.503
## - H 1 8.77 12.80 1.622
## - C 1 992.80 996.82 66.948##
## Call:
## lm(formula = PCS ~ -1 + C + H + O, data = tabela2)
##
## Coefficients:
## C H O
## 0.34876 1.14845 -0.096625 Poder Calorífico Superior e a Composição Química
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 (PC) tem correlação com a Composição Química da biomassa (teor de lignina, polissacarídeos e extrativos). Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência de cada variável no valor do PC, por meio da comparação da magnitude dos betas da equação de regressão ajustada. 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.
5.1 Banco de Dados
tabela3 <- read_excel("Dados.xlsx",
sheet = "Composicao_quimica")tabela3 %>%
kbl(caption = "Tabela dos Dados", align = 'cc') %>%
kable_classic(full_width = F, html_font = "Times", font_size = 12)| extrativos | lignina | holocelulose | PCS |
|---|---|---|---|
| 2.9434407 | 29.86403 | 67.19253 | 17.52720 |
| 2.7417815 | 29.96502 | 67.29320 | 17.64854 |
| 2.6133189 | 30.74463 | 66.64205 | 17.67364 |
| 1.0326479 | 30.52474 | 68.44261 | 17.48536 |
| 1.2695355 | 29.73810 | 68.99236 | 17.49372 |
| 0.8540925 | 29.79510 | 69.35080 | 17.58996 |
| 1.2518575 | 29.12709 | 69.62106 | 17.62762 |
| 1.3169764 | 29.30707 | 69.37595 | 17.78661 |
| 1.0532025 | 28.41579 | 70.53101 | 17.59414 |
| 1.6052158 | 29.68110 | 68.71369 | 17.07531 |
| 1.7641220 | 30.03399 | 68.20188 | 17.32218 |
| 1.8730302 | 29.15542 | 68.97155 | 17.25105 |
| 2.5166644 | 27.63171 | 69.85163 | 17.34310 |
| 1.9856993 | 27.77445 | 70.23986 | 17.29707 |
| 1.9593548 | 28.24153 | 69.79912 | 17.36820 |
| 1.8144348 | 29.25707 | 68.92849 | 17.43933 |
| 1.9256605 | 29.18541 | 68.88893 | 17.38494 |
| 1.2597507 | 30.25092 | 68.48932 | 17.42678 |
| 1.8393206 | 28.46292 | 69.69776 | 17.31381 |
| 1.9124465 | 28.74425 | 69.34330 | 17.46025 |
| 1.9287968 | 28.08158 | 69.98963 | 17.54812 |
| 2.7028249 | 29.19124 | 68.10593 | 17.80335 |
| 3.7316788 | 31.36373 | 64.90459 | 17.76569 |
| 3.3813145 | 30.27092 | 66.34777 | 17.74059 |
| 1.9874762 | 30.61163 | 67.40089 | 17.48954 |
| 2.1662125 | 31.01380 | 66.81999 | 17.41423 |
| 2.3980380 | 30.81384 | 66.78812 | 17.58577 |
5.2 Análise Descritiva
# Poder Calorífico Superior
y <- tabela3$PCS
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)
PCS<-data.frame(Med_e_sd, Coef_de_Variacao)
# Extrativos
y <- tabela3$extrativos
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)
extrativos<-data.frame(Med_e_sd, Coef_de_Variacao)
# Lignina
y <- tabela3$lignina
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)
lignina<-data.frame(Med_e_sd, Coef_de_Variacao)
# Holocelulose
y <- tabela3$holocelulose
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)
holocelulose<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(PCS,extrativos)
dados<-rbind(dados, lignina)
dados<-rbind(dados, holocelulose)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Poder Calorífico Superior (MJ/kg)",
"Extrativos", "Lignina", "Holocelulose"))
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(%) |
|---|---|---|
| Poder Calorífico Superior (MJ/kg) | 17.50 \(\pm\) 0.18 | 1.02 |
| Extrativos | 1.99 \(\pm\) 0.71 | 35.78 |
| Lignina | 29.53 \(\pm\) 1.01 | 3.43 |
| Holocelulose | 68.48 \(\pm\) 1.39 | 2.03 |
5.2.1 Poder Calorífico Superior - outlier
y <- tabela3$PCS
scale(y) # padronização da variável ## [,1]
## [1,] 0.16212574
## [2,] 0.84462282
## [3,] 0.98582911
## [4,] -0.07321808
## [5,] -0.02614931
## [6,] 0.51514147
## [7,] 0.72695091
## [8,] 1.62125742
## [9,] 0.53867585
## [10,] -2.37958750
## [11,] -0.99105897
## [12,] -1.39114346
## [13,] -0.87338706
## [14,] -1.13226526
## [15,] -0.73218077
## [16,] -0.33209628
## [17,] -0.63804324
## [18,] -0.40269942
## [19,] -1.03812773
## [20,] -0.21442437
## [21,] 0.27979765
## [22,] 1.71539494
## [23,] 1.50358551
## [24,] 1.36237922
## [25,] -0.04968370
## [26,] -0.47330257
## [27,] 0.49160709
## attr(,"scaled:center")
## [1] 17.49837
## attr(,"scaled:scale")
## [1] 0.1777867score_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
## outlier tipical
## 1 265.2.2 Extrativos - outlier
y <- tabela3$extrativos
scale(y) # padronização da variável ## [,1]
## [1,] 1.331445495
## [2,] 1.048749709
## [3,] 0.868664554
## [4,] -1.347198040
## [5,] -1.015117292
## [6,] -1.597505782
## [7,] -1.039899188
## [8,] -0.948612380
## [9,] -1.318383529
## [10,] -0.544544085
## [11,] -0.321781574
## [12,] -0.169108733
## [13,] 0.733169486
## [14,] -0.011163553
## [15,] -0.048094631
## [16,] -0.251250556
## [17,] -0.095328870
## [18,] -1.028834087
## [19,] -0.216364485
## [20,] -0.113852965
## [21,] -0.090932378
## [22,] 0.994138401
## [23,] 2.436436619
## [24,] 1.945278655
## [25,] -0.008672669
## [26,] 0.241888773
## [27,] 0.566873104
## attr(,"scaled:center")
## [1] 1.993663
## attr(,"scaled:scale")
## [1] 0.7133434score_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
## outlier tipical
## 1 265.2.3 Lignina - outlier
y <- tabela3$lignina
scale(y) # padronização da variável ## [,1]
## [1,] 0.3318672
## [2,] 0.4315091
## [3,] 1.2007107
## [4,] 0.9837564
## [5,] 0.2076259
## [6,] 0.2638626
## [7,] -0.3952338
## [8,] -0.2176548
## [9,] -1.0970325
## [10,] 0.1513780
## [11,] 0.4995639
## [12,] -0.3672771
## [13,] -1.8706458
## [14,] -1.7298169
## [15,] -1.2689705
## [16,] -0.2669823
## [17,] -0.3376925
## [18,] 0.7135992
## [19,] -1.0505316
## [20,] -0.7729587
## [21,] -1.4267870
## [22,] -0.3319350
## [23,] 1.8115447
## [24,] 0.7333263
## [25,] 1.0694917
## [26,] 1.4662866
## [27,] 1.2689963
## attr(,"scaled:center")
## [1] 29.52767
## attr(,"scaled:scale")
## [1] 1.013532score_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
## 275.2.4 Holocelulose - outlier
y <- tabela3$holocelulose
scale(y) # padronização da variável ## [,1]
## [1,] -0.927465940
## [2,] -0.854870991
## [3,] -1.324430479
## [4,] -0.025999083
## [5,] 0.370436514
## [6,] 0.628920745
## [7,] 0.823805163
## [8,] 0.647056570
## [9,] 1.479993882
## [10,] 0.169479455
## [11,] -0.199595501
## [12,] 0.355428315
## [13,] 0.990075132
## [14,] 1.270038207
## [15,] 0.952210656
## [16,] 0.324379028
## [17,] 0.295852212
## [18,] 0.007684768
## [19,] 0.879116955
## [20,] 0.623510368
## [21,] 1.089592398
## [22,] -0.268789189
## [23,] -2.577357669
## [24,] -1.536648217
## [25,] -0.777213098
## [26,] -1.196115607
## [27,] -1.219094595
## attr(,"scaled:center")
## [1] 68.47867
## attr(,"scaled:scale")
## [1] 1.38672score_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
## outlier tipical
## 1 265.3 Gráfico Box Plot
par(mfrow=c(1,4))
boxplot(tabela3$extrativos, col="blue", main="Box Plot",
ylab= "Extrativos", adj=0.5, ylim=c(0,4))
boxplot(tabela3$lignina, col="purple",main="Box Plot",
ylab= "Lignina", adj=0.5, ylim=c(26,32))
boxplot(tabela3$holocelulose,col="green2",main="Box Plot", ylab= "Holocelulose",
adj=0.5, ylim=c(64,72))
boxplot(tabela3$PCS,col="gray",main="Box Plot", ylab= "PCS (MJ/kg)",
adj=0.5, ylim=c(17.0,17.8)) 
5.4 Dispersão, Histograma e Correlação
chart.Correlation(tabela3, histogram =TRUE)
5.5 Regressão Linear Múltipla
# Ajuste do modelo
modelo3<-lm(formula =PCS~-1+extrativos+lignina+holocelulose, data=tabela3)
modelo3##
## Call:
## lm(formula = PCS ~ -1 + extrativos + lignina + holocelulose,
## data = tabela3)
##
## Coefficients:
## extrativos lignina holocelulose
## 0.2244 0.2057 0.1603# Estatísticas das estimativas dos paramentros do modelo
summary(modelo3)##
## Call:
## lm(formula = PCS ~ -1 + extrativos + lignina + holocelulose,
## data = tabela3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40514 -0.09816 0.00261 0.10095 0.34158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## extrativos 0.224352 0.050789 4.417 0.000183 ***
## lignina 0.205672 0.024158 8.514 1.03e-08 ***
## holocelulose 0.160314 0.009798 16.361 1.61e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1687 on 24 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 9.678e+04 on 3 and 24 DF, p-value: < 2.2e-16# Anova da regressão do modelo
anova(modelo3)## Analysis of Variance Table
##
## Response: PCS
## Df Sum Sq Mean Sq F value Pr(>F)
## extrativos 1 7376.7 7376.7 259062.87 < 2.2e-16 ***
## lignina 1 883.0 883.0 31011.91 < 2.2e-16 ***
## holocelulose 1 7.6 7.6 267.69 1.614e-14 ***
## Residuals 24 0.7 0.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.5.1 Medidas de Precisão para Análise do Modelo
PCS3_est<- predict(modelo3)
resid<-(tabela3$PCS-PCS3_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
Syx<-sqrt(SQResiduo/((length(tabela3$PCS))-3))
Syx_Porc<-(Syx/(mean(tabela3$PCS)))*100
Syx<-round(Syx,2)
Syx_Porc<-round(Syx_Porc, 2)
newdata33 = data.frame(Syx,Syx_Porc)
# Criar da Tabela Descritiva
knitr::kable(newdata33, align = 'cc',
caption = "Medidas de Precisão do Modelo",
col.names = c("S$_{yx}$", "S$_{yx}$%")) %>%
kable_classic(full_width = F, html_font = "Times", font_size = 16)| S\(_{yx}\) | S\(_{yx}\)% |
|---|---|
| 0.17 | 0.96 |
5.5.2 Predição do PCS aplicando o modelo
newdatas<-data.frame(extrativos= c(1.25),lignina= c(30), holocelulose= c(68))
PCS_estimado1 = predict(modelo3, newdatas, interval = "confidence");
PCS_estimado1 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 17.35194 17.23403 17.469855.5.3 Análise grafica do modelo
par(mfrow=c(2,2))
plot(modelo3, pch=20)
5.5.4 Testes de hipoteses
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(modelo3$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo3$residuals
## W = 0.98761, p-value = 0.981ks.test(modelo3$residuals, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelo3$residuals
## D = 0.37484, p-value = 0.0006472
## 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(modelo3)## lag Autocorrelation D-W Statistic p-value
## 1 0.4867611 1.023195 0
## 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(modelo3)##
## studentized Breusch-Pagan test
##
## data: modelo3
## BP = 1.326, df = 2, p-value = 0.5153ncvTest(modelo3)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.173309, Df = 1, p = 0.278725.5.5 Gráfico 1:1 e de Resíduos
Valores em vermelho estimados; valores em azul observados.
PCS3_est<-predict.lm(modelo3)
newdata3 = data.frame(tabela3, PCS3_est)
Res<-(((newdata3$PCS-newdata3$PCS3_est)/(newdata3$PCS))*100)
newdata3 = data.frame(tabela3, PCS3_est,Res)
G1<-ggplot(newdata3) +
geom_point(size=2, aes(x=PCS, y=PCS3_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
geom_point(size=2, aes(x=PCS3_est, y=PCS), color="blue2")+
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observada (MJ/kg)",
y = "PCS Estimada (Mj/kg)",
title = "Gráfico 1:1")+ # legendas
scale_x_continuous(breaks = seq(17.0, 18.0, 0.20), limits= c(17.0, 18.0))+
scale_y_continuous(breaks = seq(17.0, 18.0, 0.20), limits= c(17.0, 18.0))+
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(newdata3) +
geom_point(size=2, aes(x=PCS, y=Res)
,color="black")+ # gráfico de pontos
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observada (g/cm³)",
y = "Resíduo (%)",
title = "Gráfico de Resíduos")+ # legendas
scale_x_continuous(breaks = seq(17.0, 18.0, 0.20), limits= c(17.0, 18.0)) +
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) 
5.5.6 Dados Padronizados
tabela_padronizada3<- scale(tabela3) # padrozinação do banco de dados
tabela_padronizada3<- data.frame(tabela_padronizada3) # tornar os dados padronizados em `data.frame`
lm(formula =PCS~-1+extrativos+lignina+holocelulose, data=tabela_padronizada3) # coeficientes com dados padronizados##
## Call:
## lm(formula = PCS ~ -1 + extrativos + lignina + holocelulose,
## data = tabela_padronizada3)
##
## Coefficients:
## extrativos lignina holocelulose
## 0.2569 0.2586 NA5.5.7 Método de Seleção de modelos utilizando critérios matemáticos.
modelo3<-lm(formula =PCS~-1+extrativos+lignina+holocelulose, data=tabela3)
mod.simples3<-lm(formula=PCS~1, data=tabela3)
stepAIC(modelo3, scope=list(upper=modelo3,
lower=mod.simples3), direction = "backward")## Start: AIC=-93.27
## PCS ~ -1 + extrativos + lignina + holocelulose
##
## Df Sum of Sq RSS AIC
## <none> 0.6834 -93.266
## - extrativos 1 0.5556 1.2390 -79.201
## - lignina 1 2.0639 2.7473 -57.701
## - holocelulose 1 7.6224 8.3057 -27.830##
## Call:
## lm(formula = PCS ~ -1 + extrativos + lignina + holocelulose,
## data = tabela3)
##
## Coefficients:
## extrativos lignina holocelulose
## 0.2244 0.2057 0.16036 Poder Calorífico Superior e a Composição Imediata
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 (PC) tem correlação com a Composição Imediata da biomassa (teor de carbono fixo, voláteis e cinzas). Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência de cada variável no valor do PC, por meio da comparação da magnitude dos betas da equação de regressão ajustada. 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.
6.1 Banco de Dados
tabela4 <- read_excel("Dados.xlsx",
sheet = "Analise_imediata")tabela4 %>%
kbl(caption = "Tabela dos Dados") %>%
kable_classic(full_width = F, html_font = "Times", font_size = 12)%>%
scroll_box(height = "300px")| CF | TV | CZ | PCS |
|---|---|---|---|
| 20.9 | 78.7 | 0.4 | 20.70 |
| 10.1 | 89.9 | 0.0 | 19.76 |
| 17.2 | 81.6 | 1.3 | 19.84 |
| 20.3 | 79.6 | 0.0 | 20.44 |
| 17.8 | 81.3 | 0.9 | 18.98 |
| 14.2 | 85.6 | 0.2 | 19.91 |
| 16.6 | 80.6 | 2.9 | 19.35 |
| 19.5 | 79.4 | 1.1 | 19.16 |
| 20.7 | 77.7 | 1.6 | 22.22 |
| 18.0 | 76.1 | 5.9 | 16.18 |
| 27.1 | 71.0 | 1.9 | 20.28 |
| 22.5 | 75.9 | 1.7 | 21.10 |
| 15.3 | 84.6 | 0.1 | 20.21 |
| 9.1 | 90.9 | 0.1 | 19.91 |
| 13.3 | 86.7 | 0.1 | 20.14 |
| 21.4 | 78.4 | 0.2 | 22.01 |
| 24.7 | 70.6 | 4.7 | 18.67 |
| 19.6 | 76.2 | 4.2 | 16.55 |
| 18.3 | 78.5 | 3.2 | 18.69 |
| 20.0 | 79.6 | 0.4 | 19.41 |
| 17.8 | 82.0 | 0.2 | 19.72 |
| 17.4 | 82.4 | 0.2 | 20.11 |
| 16.0 | 83.2 | 0.7 | 19.70 |
| 14.8 | 84.9 | 0.3 | 19.84 |
| 5.4 | 93.9 | 0.8 | 19.97 |
| 19.7 | 78.4 | 1.8 | 19.85 |
| 17.6 | 76.2 | 6.1 | 16.65 |
| 15.3 | 84.1 | 0.7 | 19.31 |
| 20.0 | 79.4 | 0.6 | 22.22 |
| 15.6 | 82.6 | 1.8 | 19.88 |
| 18.2 | 80.9 | 0.9 | 20.49 |
| 17.1 | 82.0 | 1.0 | 18.92 |
| 19.3 | 79.7 | 1.0 | 20.55 |
| 19.1 | 76.9 | 4.1 | 19.99 |
| 19.6 | 80.1 | 0.3 | 19.85 |
| 16.7 | 82.8 | 1.7 | 20.20 |
| 22.3 | 74.8 | 2.9 | 19.53 |
| 19.8 | 78.0 | 2.3 | 19.59 |
| 17.6 | 81.8 | 0.6 | 19.14 |
| 16.8 | 82.6 | 0.7 | 19.87 |
6.2 Análise Descritiva
# Poder Calorífico Superior
y <- tabela4$PCS
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)
PCS<-data.frame(Med_e_sd, Coef_de_Variacao)
# Carbono Fixo
y <- tabela4$CF
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)
CF<-data.frame(Med_e_sd, Coef_de_Variacao)
# Teor de Voláteis
y <- tabela4$TV
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)
TV<-data.frame(Med_e_sd, Coef_de_Variacao)
# Teor de Cinzas
y <- tabela4$CZ
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)
CZ<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(PCS,CF)
dados<-rbind(dados, TV)
dados<-rbind(dados, CZ)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Poder Calorífico Superior (MJ/kg)",
"Carbono Fixo", "Teor de Voláteis", "Teor de Cinzas"))
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(%) |
|---|---|---|
| Poder Calorífico Superior (MJ/kg) | 19.72 \(\pm\) 1.24 | 6.30 |
| Carbono Fixo | 17.82 \(\pm\) 3.93 | 22.03 |
| Teor de Voláteis | 80.74 \(\pm\) 4.71 | 5.84 |
| Teor de Cinzas | 1.49 \(\pm\) 1.61 | 107.78 |
6.2.1 Poder Calorífico Superior - outlier
y <- tabela4$PCS
scale(y) # padronização da variável ## [,1]
## [1,] 0.786308913
## [2,] 0.030358641
## [3,] 0.094694835
## [4,] 0.577216285
## [5,] -0.596919244
## [6,] 0.150989004
## [7,] -0.299364350
## [8,] -0.452162809
## [9,] 2.008696587
## [10,] -2.848686012
## [11,] 0.448543898
## [12,] 1.107989880
## [13,] 0.392249729
## [14,] 0.150989004
## [15,] 0.335955560
## [16,] 1.839814080
## [17,] -0.846221993
## [18,] -2.551131117
## [19,] -0.830137945
## [20,] -0.251112205
## [21,] -0.001809455
## [22,] 0.311829487
## [23,] -0.017893504
## [24,] 0.094694835
## [25,] 0.199241149
## [26,] 0.102736859
## [27,] -2.470710876
## [28,] -0.331532446
## [29,] 2.008696587
## [30,] 0.126862931
## [31,] 0.617426406
## [32,] -0.645171389
## [33,] 0.665678551
## [34,] 0.215325197
## [35,] 0.102736859
## [36,] 0.384207705
## [37,] -0.154607915
## [38,] -0.106355770
## [39,] -0.468246857
## [40,] 0.118820907
## attr(,"scaled:center")
## [1] 19.72225
## attr(,"scaled:scale")
## [1] 1.243468score_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
## outlier tipical
## 5 356.2.2 Carbono Fixo - outlier
y <- tabela4$CF
scale(y) # padronização da variável ## [,1]
## [1,] 0.785246830
## [2,] -1.965982940
## [3,] -0.157304110
## [4,] 0.632400732
## [5,] -0.004458011
## [6,] -0.921534601
## [7,] -0.310150208
## [8,] 0.428605934
## [9,] 0.734298131
## [10,] 0.046490688
## [11,] 2.364656514
## [12,] 1.192836426
## [13,] -0.641316754
## [14,] -2.220726438
## [15,] -1.150803749
## [16,] 0.912618579
## [17,] 1.753272120
## [18,] 0.454080284
## [19,] 0.122913737
## [20,] 0.555977683
## [21,] -0.004458011
## [22,] -0.106355410
## [23,] -0.462996306
## [24,] -0.768688503
## [25,] -3.163277378
## [26,] 0.479554634
## [27,] -0.055406711
## [28,] -0.641316754
## [29,] 0.555977683
## [30,] -0.564893705
## [31,] 0.097439388
## [32,] -0.182778459
## [33,] 0.377657235
## [34,] 0.326708535
## [35,] 0.454080284
## [36,] -0.284675858
## [37,] 1.141887727
## [38,] 0.505028983
## [39,] -0.055406711
## [40,] -0.259201508
## attr(,"scaled:center")
## [1] 17.8175
## attr(,"scaled:scale")
## [1] 3.925517score_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
## outlier tipical
## 3 376.2.3 Teor de Voláteis - outlier
y <- tabela4$TV
scale(y) # padronização da variável ## [,1]
## [1,] -0.43272206
## [2,] 1.94300688
## [3,] 0.18242204
## [4,] -0.24181527
## [5,] 0.11878645
## [6,] 1.03089666
## [7,] -0.02969661
## [8,] -0.28423900
## [9,] -0.64484071
## [10,] -0.98423056
## [11,] -2.06603570
## [12,] -1.02665429
## [13,] 0.81877801
## [14,] 2.15512553
## [15,] 1.26422718
## [16,] -0.49635765
## [17,] -2.15088316
## [18,] -0.96301869
## [19,] -0.47514579
## [20,] -0.24181527
## [21,] 0.26726951
## [22,] 0.35211697
## [23,] 0.52181189
## [24,] 0.88241360
## [25,] 2.79148150
## [26,] -0.49635765
## [27,] -0.96301869
## [28,] 0.71271868
## [29,] -0.28423900
## [30,] 0.39454070
## [31,] 0.03393898
## [32,] 0.26726951
## [33,] -0.22060340
## [34,] -0.81453564
## [35,] -0.13575594
## [36,] 0.43696443
## [37,] -1.25998481
## [38,] -0.58120511
## [39,] 0.22484577
## [40,] 0.39454070
## attr(,"scaled:center")
## [1] 80.74
## attr(,"scaled:scale")
## [1] 4.714343score_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
## outlier tipical
## 4 366.2.4 Teor de Cinzas - outlier
y <- tabela4$CZ
scale(y) # padronização da variável ## [,1]
## [1,] -0.67875278
## [2,] -0.92783637
## [3,] -0.11831470
## [4,] -0.92783637
## [5,] -0.36739829
## [6,] -0.80329458
## [7,] 0.87801965
## [8,] -0.24285650
## [9,] 0.06849799
## [10,] 2.74614657
## [11,] 0.25531068
## [12,] 0.13076888
## [13,] -0.86556547
## [14,] -0.86556547
## [15,] -0.86556547
## [16,] -0.80329458
## [17,] 1.99889580
## [18,] 1.68754132
## [19,] 1.06483234
## [20,] -0.67875278
## [21,] -0.80329458
## [22,] -0.80329458
## [23,] -0.49194009
## [24,] -0.74102368
## [25,] -0.42966919
## [26,] 0.19303978
## [27,] 2.87068837
## [28,] -0.49194009
## [29,] -0.55421099
## [30,] 0.19303978
## [31,] -0.36739829
## [32,] -0.30512740
## [33,] -0.30512740
## [34,] 1.62527042
## [35,] -0.74102368
## [36,] 0.13076888
## [37,] 0.87801965
## [38,] 0.50439427
## [39,] -0.55421099
## [40,] -0.49194009
## attr(,"scaled:center")
## [1] 1.49
## attr(,"scaled:scale")
## [1] 1.605887score_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
## outlier tipical
## 2 386.3 Gráfico Box Plot
par(mfrow=c(1,4))
boxplot(tabela4$CF, col="blue", main="Box Plot",
ylab= "Carbono Fixo (%)", adj=0.5, ylim=c(5,30))
boxplot(tabela4$TV, col="purple",main="Box Plot",
ylab= "Teor de Voláteis (%)", adj=0.5, ylim=c(70, 95))
boxplot(tabela4$CZ,col="green2",main="Box Plot", ylab= "Teor de Cinzas (%)",
adj=0.5, ylim=c(0,6))
boxplot(tabela4$PCS,col="gray",main="Box Plot", ylab= "PCS (MJ/kg)",
adj=0.5) 
6.4 Matriz de correlação (valor de r e o pvalor)
chart.Correlation(tabela4, histogram =TRUE)
6.5 Regressão Linear Múltipla
# Ajuste do modelo
modelo4<-lm(formula =PCS~-1+CF+TV+CZ, data=tabela4); modelo4##
## Call:
## lm(formula = PCS ~ -1 + CF + TV + CZ, data = tabela4)
##
## Coefficients:
## CF TV CZ
## 0.2907 0.1876 -0.4051# Estatísticas das estimativas dos paramentros do modelo
summary(modelo4)##
## Call:
## lm(formula = PCS ~ -1 + CF + TV + CZ, data = tabela4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.74134 -0.45320 -0.05022 0.38873 2.27426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## CF 0.290744 0.031081 9.354 2.74e-11 ***
## TV 0.187586 0.006478 28.958 < 2e-16 ***
## CZ -0.405080 0.092066 -4.400 8.84e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8667 on 37 degrees of freedom
## Multiple R-squared: 0.9982, Adjusted R-squared: 0.9981
## F-statistic: 6919 on 3 and 37 DF, p-value: < 2.2e-16# Anova da regressão do modelo
anova(modelo4)## Analysis of Variance Table
##
## Response: PCS
## Df Sum Sq Mean Sq F value Pr(>F)
## CF 1 14888.0 14888.0 19821.305 < 2.2e-16 ***
## TV 1 688.7 688.7 916.870 < 2.2e-16 ***
## CZ 1 14.5 14.5 19.359 8.843e-05 ***
## Residuals 37 27.8 0.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.5.1 Medidas de Precisão para Análise do Modelo
PCS4_est<- predict(modelo4)
resid<-(tabela4$PCS-PCS4_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
Syx<-sqrt(SQResiduo/((length(tabela4$PCS))-3))
Syx_Porc<-(Syx/(mean(tabela4$PCS)))*100
Syx<-round(Syx,2)
Syx_Porc<-round(Syx_Porc, 2)
newdata33 = data.frame(Syx,Syx_Porc)
# Criar da Tabela Descritiva
knitr::kable(newdata33, align = 'cc',
caption = "Medidas de Precisão do Modelo",
col.names = c("S$_{yx}$", "S$_{yx}$%")) %>%
kable_classic(full_width = F, html_font = "Times", font_size = 16)| S\(_{yx}\) | S\(_{yx}\)% |
|---|---|
| 0.87 | 4.39 |
6.5.2 Predição do PCS aplicando o modelo
newdatas<-data.frame(CF= c(20),TV= c(78), CZ= c(0.4))
PCS_estimado4 = predict(modelo4, newdatas, interval = "confidence");
PCS_estimado4 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 20.2846 19.87058 20.698616.5.3 Análise grafica do modelo
par(mfrow=c(2,2))
plot(modelo4, pch=20)
6.5.4 Testes de hipoteses
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(modelo4$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo4$residuals
## W = 0.96481, p-value = 0.2436ks.test(modelo4$residuals, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelo4$residuals
## D = 0.13132, p-value = 0.4568
## 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(modelo4)## lag Autocorrelation D-W Statistic p-value
## 1 -0.1595594 2.31727 0.318
## 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(modelo4)##
## studentized Breusch-Pagan test
##
## data: modelo4
## BP = 1.3741, df = 2, p-value = 0.5031ncvTest(modelo4)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.2500128, Df = 1, p = 0.617076.5.5 Gráfico 1:1 e de Resíduos
Valores em vermelho estimados; valores em azul observados.
PCS4_est<-predict.lm(modelo4)
newdata4 = data.frame(tabela4, PCS4_est)
Res<-(((newdata4$PCS-newdata4$PCS4_est)/(newdata4$PCS))*100)
newdata4 = data.frame(tabela4, PCS4_est,Res)
G1<-ggplot(newdata4) +
geom_point(size=2, aes(x=PCS, y=PCS4_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
geom_point(size=2, aes(x=PCS4_est, y=PCS), color="blue2")+
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observada (MJ/kg)",
y = "PCS Estimada (Mj/kg)",
title = "Gráfico 1:1")+ # legendas
scale_x_continuous(breaks = seq(15, 25, 5), limits= c(15, 25))+
scale_y_continuous(breaks = seq(15, 25, 5), limits= c(15, 25))+
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(newdata4) +
geom_point(size=2, aes(x=PCS, y=Res)
,color="black")+ # gráfico de pontos
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observada (g/cm³)",
y = "Resíduo (%)",
title = "Gráfico de Resíduos")+ # legendas
scale_x_continuous(breaks = seq(15, 25, 5), limits= c(15, 25)) +
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) 
6.5.6 Dados Padronizados
tabela_padronizada4<- scale(tabela4) # padrozinação do banco de dados
tabela_padronizada4<- data.frame(tabela_padronizada4) # tornar os dados padronizados em `data.frame`
lm(lm(formula =PCS~-1+CF+TV+CZ, data=tabela_padronizada4)) # coeficientes com dados padronizados##
## Call:
## lm(formula = lm(formula = PCS ~ -1 + CF + TV + CZ, data = tabela_padronizada4))
##
## Coefficients:
## CF TV CZ
## 3.0413 3.2419 0.33476.5.7 Método de Seleção de modelos utilizando critérios matemáticos.
modelo4<-lm(formula =PCS~-1+CF+TV+CZ, data=tabela4)
mod.simples4<-lm(formula=PCS~1, data=tabela4)
stepAIC(modelo4, scope=list(upper=modelo4,
lower=mod.simples4), direction = "backward")## Start: AIC=-8.57
## PCS ~ -1 + CF + TV + CZ
##
## Df Sum of Sq RSS AIC
## <none> 27.79 -8.567
## - CZ 1 14.54 42.33 6.266
## - CF 1 65.73 93.52 37.971
## - TV 1 629.84 657.63 115.990##
## Call:
## lm(formula = PCS ~ -1 + CF + TV + CZ, data = tabela4)
##
## Coefficients:
## CF TV CZ
## 0.2907 0.1876 -0.40517 Poder Calorífico Superior e a Densidade Básica
A densidade básica se refere a quantidade de combustível seco (lenha) por unidade de volume saturado, sendo uma propriedade física importante na caracterização de um combustível a base de biomassa. O objetivo principal dessa modelagem é verificar se as Propriedades Energéticas do combustível (Poder Calorífico) tem correlação com a Densidade Básica. Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência da Densidade Básica no valor do Poder Calorífico, através da magnitude do beta 1 da equação de regressão ajustada. 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.
7.1 Banco de Dados
tabela5 <- read_excel("Dados.xlsx",
sheet = "PCxDb")tabela5 %>%
kbl(caption = "Tabela dos Dados") %>%
kable_classic(full_width = F, html_font = "Times", font_size = 12)%>%
scroll_box(width = "300px", height = "300px")| Nome comum | PCS | Db |
|---|---|---|
| Amargozinho | 4989 | 0.74 |
| Macucu de paca | 5075 | 0.73 |
| Melancieira | 4927 | 0.53 |
| Cajuaçu | 4456 | 0.52 |
| Cajuaçu, Cajuí | 4411 | 0.42 |
| Sucupira vermelha | 4876 | 0.67 |
| Bolsinha | 4827 | 0.61 |
| Guatambu | 4863 | 0.58 |
| Piquiá-marfim | 4742 | 0.86 |
| Maria preta | 4516 | 0.46 |
| Amapá amargoso | 4685 | 0.73 |
| Pau-marfim | 4798 | 0.91 |
| Tanimbuca | 4685 | 0.72 |
| Murici vermelho | 4844 | 0.59 |
| Murici | 4781 | 0.56 |
| Murici | 4771 | 0.48 |
| Andiroba | 4633 | 0.43 |
| Tauari da amazônia | 4721 | 0.49 |
| Pequi | 4839 | 0.61 |
| Castanha de paca | 4714 | 0.61 |
| Cedro | 4707 | 0.38 |
| Cedrorana | 4746 | 0.46 |
| Huimba negra | 4625 | 0.57 |
| Guariúba | 4848 | 0.59 |
| Coração de negro | 4813 | 0.52 |
| Castanha jacaré | 4748 | 0.84 |
| Tauarí | 4735 | 0.60 |
| Jacarandá do cerrado | 4896 | 0.77 |
| Arapari Branco | 4663 | 0.73 |
| Faveira | 4940 | 0.70 |
| Cumaru | 4866 | 0.97 |
| Cumaru | 4828 | 1.08 |
| Cumarurana | 4907 | 0.83 |
| Louro preto | 4920 | 0.48 |
|
|
4737 | 0.62 |
| Sucupira amarela | 4772 | 0.68 |
|
|
4738 | 0.57 |
| Punga colorada | 3888 | 0.39 |
| Paineira | 4565 | 0.36 |
| Paineira do cerrado | 4565 | 0.38 |
| Quarubarana | 4523 | 0.55 |
| Fruto de passarinho | 4638 | 0.52 |
| Muchiba | 4549 | 0.62 |
| Muchiba comprida | 4931 | 0.54 |
| Matá-matá | 4747 | 0.81 |
| Casca doce | 4676 | 0.73 |
| Cupiúba | 4654 | 0.69 |
|
|
4622 | 0.47 |
| Gitó | 4828 | 0.66 |
| Jaruta | 4653 | 0.71 |
| Seringueira | 4485 | 0.51 |
| Jatobá do cerrado | 4851 | 0.78 |
| Jatobá | 4792 | 0.88 |
| Jutaí | 4743 | 0.78 |
| Angelim da mata | 4828 | 0.66 |
| Angelim rajado | 4837 | 0.67 |
| Ucuuba puna | 4645 | 0.69 |
| Ucuubarana | 4792 | 0.64 |
| Caroba | 4696 | 0.35 |
| Pau santo | 4747 | 0.46 |
| Pau santo | 4882 | 0.58 |
| Mangaba brava | 4788 | 0.74 |
| Macucu fofo | 4761 | 0.88 |
| Louro aritu | 4770 | 0.79 |
| Louro chumbo | 4889 | 1.04 |
| Ingá-cumaru | 4680 | 0.68 |
| Maçranduba | 4793 | 0.92 |
| Machin sapote | 4110 | 0.48 |
| Sapote | 4062 | 0.42 |
| Itaúba | 5263 | 0.70 |
| Lacre | 4777 | 0.65 |
| Lacre | 4626 | 0.57 |
|
|
4700 | 0.52 |
| Abiurana | 4564 | 0.88 |
| Louro inhamui | 5150 | 0.47 |
| Ucuubarana | 4827 | 0.42 |
| Cabelo de negro | 4926 | 0.50 |
| Bate caixa | 4695 | 0.43 |
| Faveira-folha-fina | 4647 | 0.72 |
| Coração de negro | 4744 | 0.42 |
| Angelim | 4861 | 0.81 |
| Macacaúba | 4987 | 0.75 |
| Abiurana | 4878 | 0.90 |
| Grão de galo | 4779 | 0.70 |
|
|
4752 | 0.20 |
| Sucupira branca | 4953 | 0.73 |
| Mandioqueira | 4398 | 0.63 |
| Pau terra folha grande | 4736 | 0.69 |
| Pau terra liso | 4725 | 0.66 |
| Mandioqueira | 4626 | 0.66 |
| Pau terra roxo | 4710 | 0.69 |
| Sapotilho | 4334 | 0.46 |
|
|
4667 | 0.47 |
| Mandiocão do cerrado | 4740 | 0.68 |
| Morototó | 4556 | 0.40 |
| Carvoeiro | 4849 | 0.72 |
| Cardeiro | 4709 | 0.59 |
| Marupá | 4627 | 0.35 |
| Quina do cerrado | 4756 | 0.72 |
| Barbatimão | 4816 | 0.55 |
| Laranjeira do cerrado | 4755 | 0.49 |
| Coração de negro | 4904 | 0.97 |
| Ipê | 4760 | 0.62 |
| Pau d’arco | 4882 | 0.87 |
| Ipê | 4823 | 0.69 |
| Ipê | 4957 | 1.05 |
| Tachi preto | 4667 | 0.51 |
| Breu sucuruba | 4606 | 0.44 |
| Ucuúba grande | 4574 | 0.50 |
| Pau doce | 4736 | 0.57 |
|
|
4680 | 0.40 |
| Gomeira | 4713 | 0.49 |
7.2 Análise Descritiva
# Poder Calorífico
y <- tabela5$PCS
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.0f $\\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)
PCS2<-data.frame(Med_e_sd, Coef_de_Variacao)
# Densidade Básica
y <- tabela5$Db
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)
Db<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(Db,PCS2)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Densidade (g/cm$^{3}$)", "Poder Calorífico (MJ/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(%) |
|---|---|---|
| Densidade (g/cm\(^{3}\)) | 0.63 \(\pm\) 0.17 | 27.12 |
| Poder Calorífico (MJ/kg) | 4733 \(\pm\) 187.45 | 3.96 |
7.2.1 Poder Calorífico Superior - outlier
y <- tabela5$PCS
scale(y) # padronização da variável ## [,1]
## [1,] 1.36704612
## [2,] 1.82582415
## [3,] 1.03629917
## [4,] -1.47631074
## [5,] -1.71636901
## [6,] 0.76423313
## [7,] 0.50283634
## [8,] 0.69488296
## [9,] 0.04939294
## [10,] -1.15623304
## [11,] -0.25468087
## [12,] 0.34813212
## [13,] -0.25468087
## [14,] 0.59352502
## [15,] 0.25744344
## [16,] 0.20409716
## [17,] -0.53208154
## [18,] -0.06263425
## [19,] 0.56685188
## [20,] -0.09997665
## [21,] -0.13731905
## [22,] 0.07073145
## [23,] -0.57475856
## [24,] 0.61486354
## [25,] 0.42815155
## [26,] 0.08140071
## [27,] 0.01205054
## [28,] 0.87092569
## [29,] -0.37204269
## [30,] 1.10564933
## [31,] 0.71088684
## [32,] 0.50817097
## [33,] 0.92960660
## [34,] 0.99895677
## [35,] 0.02271980
## [36,] 0.20943179
## [37,] 0.02805443
## [38,] -4.50637958
## [39,] -0.89483626
## [40,] -0.89483626
## [41,] -1.11889064
## [42,] -0.50540840
## [43,] -0.98019031
## [44,] 1.05763768
## [45,] 0.07606608
## [46,] -0.30269252
## [47,] -0.42005434
## [48,] -0.59076245
## [49,] 0.50817097
## [50,] -0.42538897
## [51,] -1.32160652
## [52,] 0.63086742
## [53,] 0.31612435
## [54,] 0.05472757
## [55,] 0.50817097
## [56,] 0.55618263
## [57,] -0.46806600
## [58,] 0.31612435
## [59,] -0.19599996
## [60,] 0.07606608
## [61,] 0.79624090
## [62,] 0.29478584
## [63,] 0.15075088
## [64,] 0.19876253
## [65,] 0.83358329
## [66,] -0.28135401
## [67,] 0.32145898
## [68,] -3.32209211
## [69,] -3.57815427
## [70,] 2.82873426
## [71,] 0.23610493
## [72,] -0.56942393
## [73,] -0.17466144
## [74,] -0.90017089
## [75,] 2.22592127
## [76,] 0.50283634
## [77,] 1.03096454
## [78,] -0.20133459
## [79,] -0.45739674
## [80,] 0.06006220
## [81,] 0.68421370
## [82,] 1.35637686
## [83,] 0.77490238
## [84,] 0.24677419
## [85,] 0.10273922
## [86,] 1.17499950
## [87,] -1.78571917
## [88,] 0.01738517
## [89,] -0.04129574
## [90,] -0.56942393
## [91,] -0.12131516
## [92,] -2.12713538
## [93,] -0.35070418
## [94,] 0.03872369
## [95,] -0.94284791
## [96,] 0.62019816
## [97,] -0.12664979
## [98,] -0.56408931
## [99,] 0.12407774
## [100,] 0.44415543
## [101,] 0.11874311
## [102,] 0.91360272
## [103,] 0.14541625
## [104,] 0.79624090
## [105,] 0.48149783
## [106,] 1.19633801
## [107,] -0.35070418
## [108,] -0.67611650
## [109,] -0.84682460
## [110,] 0.01738517
## [111,] -0.28135401
## [112,] -0.10531128
## attr(,"scaled:center")
## [1] 4732.741
## attr(,"scaled:scale")
## [1] 187.4545score_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
## outlier tipical
## 6 1067.2.2 Densidade Básica - outlier
y <- tabela5$Db
scale(y) # padronização da variável ## [,1]
## [1,] 0.66235109
## [2,] 0.60356874
## [3,] -0.57207820
## [4,] -0.63086055
## [5,] -1.21868402
## [6,] 0.25087466
## [7,] -0.10181942
## [8,] -0.27816646
## [9,] 1.36773925
## [10,] -0.98355463
## [11,] 0.60356874
## [12,] 1.66165099
## [13,] 0.54478639
## [14,] -0.21938412
## [15,] -0.39573116
## [16,] -0.86598993
## [17,] -1.15990167
## [18,] -0.80720759
## [19,] -0.10181942
## [20,] -0.10181942
## [21,] -1.45381340
## [22,] -0.98355463
## [23,] -0.33694881
## [24,] -0.21938412
## [25,] -0.63086055
## [26,] 1.25017456
## [27,] -0.16060177
## [28,] 0.83869813
## [29,] 0.60356874
## [30,] 0.42722170
## [31,] 2.01434507
## [32,] 2.66095089
## [33,] 1.19139221
## [34,] -0.86598993
## [35,] -0.04303708
## [36,] 0.30965701
## [37,] -0.33694881
## [38,] -1.39503106
## [39,] -1.57137810
## [40,] -1.45381340
## [41,] -0.45451350
## [42,] -0.63086055
## [43,] -0.04303708
## [44,] -0.51329585
## [45,] 1.07382752
## [46,] 0.60356874
## [47,] 0.36843935
## [48,] -0.92477228
## [49,] 0.19209231
## [50,] 0.48600405
## [51,] -0.68964289
## [52,] 0.89748048
## [53,] 1.48530395
## [54,] 0.89748048
## [55,] 0.19209231
## [56,] 0.25087466
## [57,] 0.36843935
## [58,] 0.07452762
## [59,] -1.63016044
## [60,] -0.98355463
## [61,] -0.27816646
## [62,] 0.66235109
## [63,] 1.48530395
## [64,] 0.95626282
## [65,] 2.42582150
## [66,] 0.30965701
## [67,] 1.72043333
## [68,] -0.86598993
## [69,] -1.21868402
## [70,] 0.42722170
## [71,] 0.13330997
## [72,] -0.33694881
## [73,] -0.63086055
## [74,] 1.48530395
## [75,] -0.92477228
## [76,] -1.21868402
## [77,] -0.74842524
## [78,] -1.15990167
## [79,] 0.54478639
## [80,] -1.21868402
## [81,] 1.07382752
## [82,] 0.72113344
## [83,] 1.60286864
## [84,] 0.42722170
## [85,] -2.51189565
## [86,] 0.60356874
## [87,] 0.01574527
## [88,] 0.36843935
## [89,] 0.19209231
## [90,] 0.19209231
## [91,] 0.36843935
## [92,] -0.98355463
## [93,] -0.92477228
## [94,] 0.30965701
## [95,] -1.33624871
## [96,] 0.54478639
## [97,] -0.21938412
## [98,] -1.63016044
## [99,] 0.54478639
## [100,] -0.45451350
## [101,] -0.80720759
## [102,] 2.01434507
## [103,] -0.04303708
## [104,] 1.42652160
## [105,] 0.36843935
## [106,] 2.48460385
## [107,] -0.68964289
## [108,] -1.10111932
## [109,] -0.74842524
## [110,] -0.33694881
## [111,] -1.33624871
## [112,] -0.80720759
## attr(,"scaled:center")
## [1] 0.6273214
## attr(,"scaled:scale")
## [1] 0.1701191score_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
## outlier tipical
## 6 1067.3 Gráfico Box Plot
par(mfrow=c(1,2))
boxplot(tabela5$Db, col="blue", main="Box Plot",
ylab= "Densidade Básica (g/cm³)", adj=0.5, ylim=c(0.2,1.2))
boxplot(tabela5$PCS, col="purple",main="Box Plot",
ylab= "Poder Calorífico (MJ/kg)", adj=0.5, ylim=c(4000,5500))
7.4 Matriz de correlação (valor de r e o pvalor)
chart.Correlation(tabela5[-1], histogram =TRUE)
7.5 Regressão Linear Múltipla
# Ajuste do modelo
modelo5<-lm(formula =PCS~Db, data=tabela5); modelo5##
## Call:
## lm(formula = PCS ~ Db, data = tabela5)
##
## Coefficients:
## (Intercept) Db
## 4452.2 447.2# Estatísticas das estimativas dos paramentros do modelo
summary(modelo5)##
## Call:
## lm(formula = PCS ~ Db, data = tabela5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -738.61 -77.19 14.12 85.70 497.76
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4452.20 62.39 71.358 < 2e-16 ***
## Db 447.20 96.02 4.657 9.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 172.1 on 110 degrees of freedom
## Multiple R-squared: 0.1647, Adjusted R-squared: 0.1571
## F-statistic: 21.69 on 1 and 110 DF, p-value: 9.018e-06# Anova da regressão do modelo
anova(modelo5)## Analysis of Variance Table
##
## Response: PCS
## Df Sum Sq Mean Sq F value Pr(>F)
## Db 1 642439 642439 21.691 9.018e-06 ***
## Residuals 110 3258010 29618
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17.5.1 Medidas de Precisão para Análise do Modelo
PCS_est<-predict.lm(modelo5)
newdata1 = data.frame(tabela5, PCS_est)
Res<-(((newdata1$PCS-newdata1$PCS_est)/(newdata1$PCS))*100)
newdata1 = data.frame(tabela5, PCS_est,Res)
resid<-(newdata1$PCS-newdata1$PCS_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
SQTotal<-var(newdata1$PCS)*((length(newdata1$PCS))-1)
R2<-(1-(SQResiduo/SQTotal))
a<-(length(newdata1$PCS)-1)/(length(newdata1$PCS)-3)
R2adj<-(1-(a*(1-R2)))
Syx<-sqrt(SQResiduo/((length(newdata1$PCS))-2))
Syx_Porc<-(Syx/(mean(newdata1$PCS)))*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}\)% |
|---|---|---|---|
| 0.16 | 0.15 | 172.1 | 3.64 |
7.5.2 Predição do PCS aplicando o modelo
newdatas<-data.frame(Db= c(0.50))
PCS_estimado1 = predict(modelo5, newdatas, interval = "confidence");
PCS_estimado1 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 4675.803 4635.484 4716.1227.5.3 Análise grafica do modelo
par(mfrow=c(2,2))
plot(modelo5, pch=20)
7.5.4 Testes de hipoteses
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(modelo5$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo5$residuals
## W = 0.90824, p-value = 1.112e-06ks.test(modelo5$residuals, "pnorm")## Warning in ks.test.default(modelo5$residuals, "pnorm"): ties should not be
## present for the Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: modelo5$residuals
## D = 0.53571, p-value < 2.2e-16
## 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(modelo5)## lag Autocorrelation D-W Statistic p-value
## 1 0.09057317 1.805312 0.292
## 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(modelo5)##
## studentized Breusch-Pagan test
##
## data: modelo5
## BP = 5.6185, df = 1, p-value = 0.01777ncvTest(modelo5)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 17.91577, Df = 1, p = 2.309e-057.5.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=PCS, y=PCS_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
geom_point(size=2, aes(x=PCS_est, y=PCS), color="blue2")+
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observado (MJ/kg)",
y = "PCS Estimado (MJ/kg)",
title = "Gráfico 1:1")+ # legendas
scale_x_continuous(breaks = seq(3800, 5600, 200), limits= c(3800, 5600)) + # eixo x
scale_y_continuous(breaks = seq(3800, 5600, 200), limits= c(3800, 5600)) + # eixo y
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=PCS, y=Res)
,color="black")+ # gráfico de pontos
theme_classic()+ # aparencia clean do gráfico
labs(x = "PCS Observado (MJ/kg)",
y = "Resíduo (%)",
title = "Gráfico de Resíduos")+ # legendas
scale_x_continuous(breaks = seq(3800, 5600, 400), limits= c(3800, 5600)) +
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) 
7.6 Regressão Linear
grafico<-ggplot(data = tabela5, #banco de dados
mapping = aes(x = Db, y = PCS)) + # 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.20, label.y=5300)+ # localização xy da formula
scale_x_continuous(breaks = seq(0.20, 1.20, 0.20), limits= c(0.20, 1.20)) +
scale_y_continuous(breaks = seq(3800, 5400, 400), limits= c(3800, 5400))+
labs(x = "Densidade Básica (g/cm³)", # legenda do eixo x
y = "Poder Calorífico Superior (MJ/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
8 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.
No artigo de QUIRINO et al. (2005) vemos informações equivocadas sobre a relação do poder calorífico e a densidade básica.
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.
QUIRINO, Waldir F; VALE, AT do; ANDRADE, APA de; ABREU, Vera Lúcia
Silva; AZEVEDO, AC dos S. Poder calorı́fico da madeira e de
materiais ligno-celulósicos. Revista da
madeira, vol. 89, no. 100, p. 100–106, 2005.
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.