Home
Ricardo de Oliveira Perdiz 2021-02-17 (Atualizado por último: 2021-07-01)
Primeiro, é importante que você tenha instalado os pacotes listados abaixo para poder executar este tutorial. Executa os comandos abaixo para instalar os pacotes:
libs <- c("data.table", "dplyr", "purrr", "tidyr", "caret", "tidymodels", "pls")
install.packages(libs)Após instalar todos os pacotes acima, carregue-os executando os comandos abaixo:
library("NIRtools")
library("caret")
library("tidymodels")
library("purrr")
library("tidyr")
library("MASS")
library("pls")Vou utilizar o conjunto de dados nir_data, um subconjunto proveniente
de minha tese de doutorado (R. O. Perdiz 2020), e que acompanha o pacote
NIRtools (Ricardo Oliveira Perdiz 2021) para servir como exemplo para
as análises neste tutorial. Vale ressaltar que esse conjunto de dados é
extremamente pequeno, por isso desconsidere os resultados das
análises. Ele servirá apenas de guia para abordar os procedimentos aqui
abordados.
Para eu poder ler o conjunto de dados que acompanha o pacote, eu faço
uso do comando abaixo, system.file() para poder descobrir onde está
este arquivo no meu computador. Você pode executá-lo no seu computador e
buscar o arquivo para poder copiá-lo e colocá-lo na sua pasta de
trabalho para facilitar o entendimento dos próximos passos.
system.file("extdata", "nir_raw.csv", package = "NIRtools")## [1] "/Library/Frameworks/R.framework/Versions/4.1/Resources/library/NIRtools/extdata/nir_raw.csv"
Seguindo então com a leitura dos dados, uso o system.file() para
providenciar o local do arquivo para utilizá-lo como primeiro argumento
da função read.table(), que é a função utilizada para ler o arquivo
.csv:
caminho <- system.file("extdata", "nir_data.csv", package = "NIRtools")
dado1 <- read.table(caminho, sep = "\t", header = TRUE)
head(dado1, 6)[,1:10]## especimenid SP1 face coletor number X10001.03 X3999.64
## 1 10194 P. aracouchini abaxial Perdiz, R.O. 2856 0.2776265 0.5961316
## 2 10194 P. aracouchini adaxial Perdiz, R.O. 2856 0.2542301 0.5672824
## 3 10194 P. aracouchini abaxial Perdiz, R.O. 2856 0.2730012 0.5813973
## 4 10194 P. aracouchini adaxial Perdiz, R.O. 2856 0.2582242 0.5764732
## 5 10194 P. aracouchini abaxial Perdiz, R.O. 2856 0.2776265 0.5961316
## 6 10194 P. aracouchini adaxial Perdiz, R.O. 2856 0.2542301 0.5672824
## X4003.497 X4007.354 X4011.211
## 1 0.5964082 0.5962480 0.5961008
## 2 0.5674334 0.5673867 0.5674006
## 3 0.5817484 0.5816703 0.5815431
## 4 0.5767084 0.5766730 0.5765890
## 5 0.5964082 0.5962480 0.5961008
## 6 0.5674334 0.5673867 0.5674006
Dimensões do objeto dado1:
dim(dado1)## [1] 48 1562
Variável identificadora de cada amostra neste objeto se encontra na
variável especimenid. Em seus dados, ela pode ser outra variável, pode
ser um nome, pode ser um código etc. No meu caso, são números. Números
iguais demonstram que se trata do mesmo espécime, podendo ser leituras
adicionais de uma mesma face de um folíolo (meus dados são de plantas!),
ou simplesmente leituras de duas faces diferentes do folíolo:
table(dado1$especimenid)##
## 10194 10196 10197 12647 12648 12650
## 8 8 8 8 8 8
A variável identificadora do táxon está presente na coluna SP1:
table(dado1$SP1)##
## P. aracouchini P. calanense
## 24 24
Vamos utilizar aqui as funcionalidades do pacote tidymodels (Kuhn e
Wickham 2020b). Vou estabelecer uma seed para que vocês cheguem ao
mesmo resultado, caso utilizem estas “receitas” durante o aprendizado do
tutorial.
set.seed(221015)Para separar os dados em conjunto treino e teste, utilizamos três funções:
-
initial_split()do pacotersample(Kuhn, Chow, e Wickham 2020), que é carregado automaticamente ao se chamar o pacotetidymodels(Kuhn e Wickham 2020b) - esta função cria índices para separar o conjunto de dados emtreinoeteste. Por padrão, a função divide o conjunto de dados em 75% paratreinoe 25% parateste. Caso você queira outra proporção, altere o argumentoprop. Em caso de dúvida, veja o?desta função executando o comando?initial_split.
dado1_split <- initial_split(dado1, strata = "SP1")- funções
training()etesting(), ambas do mesmo pacotersample- ambas recebem o objeto criado com a funçãoinitial_split(),dado1_split, e dividem o treino em conjuntostreinoeteste, respectivamente:
treino <- training(dado1_split)
table(treino$SP1)##
## P. aracouchini P. calanense
## 18 18
teste <- testing(dado1_split)
table(teste$SP1)##
## P. aracouchini P. calanense
## 6 6
Vou utilizar abaixo o ambiente de pré-processamento do pacote recipes
(Kuhn e Wickham 2020a). Ele facilita a criação de “receitas” que
permitem a criação e documentação de conjuntos de dados para serem
utilizados posteriormente em análises do tidymodels, permitindo ao
usuário um fluxo contínuo de análise de dados e a facilidade na
reutilização do mesmo conjunto de dados para diversas análises.
Recomendo o aprendizado destas novas (algumas não tão novas, como é o
caso do pacote caret (Kuhn 2020)) ferramentas de pré-processamento e
análise de dados em R.
Separei o tópico de análises em Sem receitas e Com receitas. Caso
você prossiga nesta seção, você poderá utilizar os objetos aqui
construídos na seção Com receitas de todas as análises.
A primeira função desta abordagem é a função recipe(), em que usamos
uma fórmula SP1 ~ . para dizer que SP1 é nossa variável resposta e
as demais são as variáveis preditoras, representadas pelo . na
fórmula. Porém, sabemos que nosso conjunto de dados dado1 possui
outras variáveis que não são preditoras, como face e coletor.
Aí entra o papel da função update_role() em que informamos o papel de
cada variável extra em nosso conjunto de dados. Por exemplo, as
variáveis coletore number são identificadores da coleta. Então,
iniciamos a função com o objeto contendo a receita, dado1_receita, e
colocamos estas duas variáveis, sem aspas, em seguida, seguido do
argumento new_role, em que colocamos um texto informando que papel é
esse. Ao utilziarmos esta receita nos passos de análises, essas
variáveis não serão analisadas pois terão a informação desses papéis
que elas têm no conjunto de dados.
update_role(dado1_receita, coletor, number, new_role = "dados do espécime")# receita para analise
dado1_receita <- recipes::recipe(SP1 ~ ., data = treino)
dado1_receita <- update_role(dado1_receita, especimenid, new_role = "id")
dado1_receita <- update_role(dado1_receita, coletor, number, new_role = "dados do espécime")
dado1_receita <- update_role(dado1_receita, face, new_role = "identificador da face do foliolo")
dado1_receita <- prep(dado1_receita)
dado1_receita## Data Recipe
##
## Inputs:
##
## role #variables
## dados do espécime 2
## id 1
## identificador da face do foliolo 1
## outcome 1
## predictor 1557
##
## Training data contained 36 data points and no missing data.
Vejam que os comandos acima ficaram muito repetitivos, pois a receita
dado1_receita aparece múltiplas vezes. Caso utilizássemos o
encadeamento de ações proporcionado pelo operador %>% do pacote
magrittr (Bache e Wickham 2020), o texto ficaria mais limpo. Porém,
devo admitir que isto é uma questão de preferência. Aos interessados,
sugiro ler o texto Base R
vs. Tidyverse
de nosso Curso básico de Introdução ao R
para aprender um pouquinho mais sobre este operador, presente no pacote
magrittr.
dado1_receita <-
recipes::recipe(SP1 ~ ., data = treino) %>%
update_role(., especimenid, new_role = "id") %>%
update_role(., coletor, number, new_role = "dados do espécime") %>%
update_role(., face, new_role = "identificador da face do foliolo") %>%
prep()
dado1_receitateste_processado <- bake(dado1_receita, new_data = teste)
teste_processado## # A tibble: 12 x 1,562
## especimenid face coletor number X10001.03 X3999.64 X4003.497 X4007.354
## <int> <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 10194 abaxial Perdiz, R.… 2856 0.278 0.596 0.596 0.596
## 2 10194 abaxial Perdiz, R.… 2856 0.278 0.596 0.596 0.596
## 3 10197 abaxial Perdiz, R.… 2859 0.271 0.681 0.681 0.681
## 4 10197 adaxial Perdiz, R.… 2859 0.274 0.662 0.662 0.661
## 5 10197 adaxial Perdiz, R.… 2859 0.274 0.662 0.662 0.661
## 6 10197 adaxial Perdiz, R.… 2859 0.275 0.710 0.710 0.709
## 7 12647 adaxial Perdiz, R.… 3236 0.329 0.714 0.715 0.715
## 8 12647 adaxial Perdiz, R.… 3236 0.320 0.735 0.736 0.736
## 9 12647 abaxial Perdiz, R.… 3236 0.340 0.724 0.724 0.724
## 10 12647 adaxial Perdiz, R.… 3236 0.320 0.735 0.736 0.736
## 11 12650 adaxial Perdiz, R.… 3239 0.342 0.744 0.744 0.744
## 12 12650 abaxial Perdiz, R.… 3239 0.342 0.713 0.714 0.714
## # … with 1,554 more variables: X4011.211 <dbl>, X4015.068 <dbl>,
## # X4018.925 <dbl>, X4022.781 <dbl>, X4026.638 <dbl>, X4030.495 <dbl>,
## # X4034.352 <dbl>, X4038.209 <dbl>, X4042.066 <dbl>, X4045.923 <dbl>,
## # X4049.78 <dbl>, X4053.637 <dbl>, X4057.494 <dbl>, X4061.351 <dbl>,
## # X4065.208 <dbl>, X4069.065 <dbl>, X4072.922 <dbl>, X4076.779 <dbl>,
## # X4080.635 <dbl>, X4084.492 <dbl>, X4088.349 <dbl>, X4092.206 <dbl>,
## # X4096.063 <dbl>, X4099.92 <dbl>, X4103.777 <dbl>, X4107.634 <dbl>,
## # X4111.491 <dbl>, X4115.348 <dbl>, X4119.205 <dbl>, X4123.062 <dbl>,
## # X4126.918 <dbl>, X4130.775 <dbl>, X4134.632 <dbl>, X4138.489 <dbl>,
## # X4142.346 <dbl>, X4146.203 <dbl>, X4150.06 <dbl>, X4153.917 <dbl>,
## # X4157.774 <dbl>, X4161.631 <dbl>, X4165.488 <dbl>, X4169.345 <dbl>,
## # X4173.202 <dbl>, X4177.059 <dbl>, X4180.916 <dbl>, X4184.772 <dbl>,
## # X4188.629 <dbl>, X4192.486 <dbl>, X4196.343 <dbl>, X4200.2 <dbl>,
## # X4204.057 <dbl>, X4207.914 <dbl>, X4211.771 <dbl>, X4215.628 <dbl>,
## # X4219.485 <dbl>, X4223.342 <dbl>, X4227.199 <dbl>, X4231.056 <dbl>,
## # X4234.913 <dbl>, X4238.77 <dbl>, X4242.626 <dbl>, X4246.483 <dbl>,
## # X4250.34 <dbl>, X4254.197 <dbl>, X4258.054 <dbl>, X4261.911 <dbl>,
## # X4265.768 <dbl>, X4269.625 <dbl>, X4273.482 <dbl>, X4277.339 <dbl>,
## # X4281.196 <dbl>, X4285.053 <dbl>, X4288.91 <dbl>, X4292.767 <dbl>,
## # X4296.624 <dbl>, X4300.48 <dbl>, X4304.337 <dbl>, X4308.194 <dbl>,
## # X4312.051 <dbl>, X4315.908 <dbl>, X4319.765 <dbl>, X4323.622 <dbl>,
## # X4327.479 <dbl>, X4331.336 <dbl>, X4335.193 <dbl>, X4339.05 <dbl>,
## # X4342.907 <dbl>, X4346.764 <dbl>, X4350.621 <dbl>, X4354.478 <dbl>,
## # X4358.334 <dbl>, X4362.191 <dbl>, X4366.048 <dbl>, X4369.905 <dbl>,
## # X4373.762 <dbl>, X4377.619 <dbl>, X4381.476 <dbl>, X4385.333 <dbl>,
## # X4389.19 <dbl>, X4393.047 <dbl>, …
juice(dado1_receita)## # A tibble: 36 x 1,562
## especimenid face coletor number X10001.03 X3999.64 X4003.497 X4007.354
## <int> <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 10194 adaxial Perdiz, R.… 2856 0.254 0.567 0.567 0.567
## 2 10194 abaxial Perdiz, R.… 2856 0.273 0.581 0.582 0.582
## 3 10194 adaxial Perdiz, R.… 2856 0.258 0.576 0.577 0.577
## 4 10194 adaxial Perdiz, R.… 2856 0.254 0.567 0.567 0.567
## 5 10194 abaxial Perdiz, R.… 2856 0.273 0.581 0.582 0.582
## 6 10194 adaxial Perdiz, R.… 2856 0.258 0.576 0.577 0.577
## 7 10196 abaxial Perdiz, R.… 2858 0.274 0.681 0.681 0.681
## 8 10196 adaxial Perdiz, R.… 2858 0.278 0.673 0.672 0.672
## 9 10196 abaxial Perdiz, R.… 2858 0.276 0.703 0.703 0.703
## 10 10196 adaxial Perdiz, R.… 2858 0.267 0.712 0.712 0.711
## # … with 26 more rows, and 1,554 more variables: X4011.211 <dbl>,
## # X4015.068 <dbl>, X4018.925 <dbl>, X4022.781 <dbl>, X4026.638 <dbl>,
## # X4030.495 <dbl>, X4034.352 <dbl>, X4038.209 <dbl>, X4042.066 <dbl>,
## # X4045.923 <dbl>, X4049.78 <dbl>, X4053.637 <dbl>, X4057.494 <dbl>,
## # X4061.351 <dbl>, X4065.208 <dbl>, X4069.065 <dbl>, X4072.922 <dbl>,
## # X4076.779 <dbl>, X4080.635 <dbl>, X4084.492 <dbl>, X4088.349 <dbl>,
## # X4092.206 <dbl>, X4096.063 <dbl>, X4099.92 <dbl>, X4103.777 <dbl>,
## # X4107.634 <dbl>, X4111.491 <dbl>, X4115.348 <dbl>, X4119.205 <dbl>,
## # X4123.062 <dbl>, X4126.918 <dbl>, X4130.775 <dbl>, X4134.632 <dbl>,
## # X4138.489 <dbl>, X4142.346 <dbl>, X4146.203 <dbl>, X4150.06 <dbl>,
## # X4153.917 <dbl>, X4157.774 <dbl>, X4161.631 <dbl>, X4165.488 <dbl>,
## # X4169.345 <dbl>, X4173.202 <dbl>, X4177.059 <dbl>, X4180.916 <dbl>,
## # X4184.772 <dbl>, X4188.629 <dbl>, X4192.486 <dbl>, X4196.343 <dbl>,
## # X4200.2 <dbl>, X4204.057 <dbl>, X4207.914 <dbl>, X4211.771 <dbl>,
## # X4215.628 <dbl>, X4219.485 <dbl>, X4223.342 <dbl>, X4227.199 <dbl>,
## # X4231.056 <dbl>, X4234.913 <dbl>, X4238.77 <dbl>, X4242.626 <dbl>,
## # X4246.483 <dbl>, X4250.34 <dbl>, X4254.197 <dbl>, X4258.054 <dbl>,
## # X4261.911 <dbl>, X4265.768 <dbl>, X4269.625 <dbl>, X4273.482 <dbl>,
## # X4277.339 <dbl>, X4281.196 <dbl>, X4285.053 <dbl>, X4288.91 <dbl>,
## # X4292.767 <dbl>, X4296.624 <dbl>, X4300.48 <dbl>, X4304.337 <dbl>,
## # X4308.194 <dbl>, X4312.051 <dbl>, X4315.908 <dbl>, X4319.765 <dbl>,
## # X4323.622 <dbl>, X4327.479 <dbl>, X4331.336 <dbl>, X4335.193 <dbl>,
## # X4339.05 <dbl>, X4342.907 <dbl>, X4346.764 <dbl>, X4350.621 <dbl>,
## # X4354.478 <dbl>, X4358.334 <dbl>, X4362.191 <dbl>, X4366.048 <dbl>,
## # X4369.905 <dbl>, X4373.762 <dbl>, X4377.619 <dbl>, X4381.476 <dbl>,
## # X4385.333 <dbl>, X4389.19 <dbl>, X4393.047 <dbl>, …
Para derivar espectros, precisaremos do pacote prospectr (Stevens,
Ramirez-Lopez, e Hans 2020). Primeiro, instale-o:
install.packages("prospectr")Agora, lembre-se de carregá-lo:
library("prospectr")## prospectr version 0.2.1 -- 'seville'
## check the github repository at: http://github.com/l-ramirez-lopez/prospectr
##
## Attaching package: 'prospectr'
## The following object is masked from 'package:pls':
##
## msc
Tomando o conjunto dado1 como exemplo, precisamos selecionar apenas as
variáveis NIR do conjunto de dados. Em nosso caso, as variáveis NIR
possuem a letra X precedendo o nome da variável. Podemos utilizar
tanto a função select() do pacote dplyr (Wickham et al. 2021):
variaveis_nir <- dplyr::select(dado1, starts_with("X"))… quanto mesclar a função grep() com o indexador [ no pacote base do
R:
variaveis_nir <- dado1[, grep("^X", names(dado1))]Para derivar as variáveis NIR, utilizamos a função savitzkyGolay() do
pacote prospectr. A função pede a matriz, ou data.frame, contendo
apenas as variáveis de interesse no argumento X; ordem da derivada no
argumento m; ordem do polinômio no argumento p; e a janela no
argumento w. Recomendo que vejam o ? da função para entender como
ela funciona:
?savitzkyGolayderivada_m1_p3_w11 <- savitzkyGolay(
X = variaveis_nir,
m = 1,
p = 3,
w = 11
)
derivada_m1_p3_w11[1:5, 1:5]## X4015.068 X4018.925 X4022.781 X4026.638 X4030.495
## [1,] -0.01928127 -0.0010994357 -0.001444311 -0.001804412 -0.002127605
## [2,] -0.01883717 -0.0009765530 -0.001368026 -0.001742961 -0.002063381
## [3,] -0.01861778 -0.0009799171 -0.001270707 -0.001560611 -0.001826906
## [4,] -0.01915667 -0.0009731902 -0.001330068 -0.001678564 -0.001959641
## [5,] -0.01928127 -0.0010994357 -0.001444311 -0.001804412 -0.002127605
Uma maneira de derivar em uma escala maior é criar uma tabela com
combinações de valores e rodar um for loop para todas essas
combinações. Temos então:
# Combinacoes
### tamanho da janela sob o argumento `w`.
## valores para tamanho da janela
w_vals <- seq(3, 101, by = 2) # valores de w devem ser ímpares
### ordem do polinômio sob o argumento `p`
p_vals <- 2:3
### ordem de diferenciação sob o argumento `m`
m_vals <- 1:2
# Create a list with these values
values_list <- list(w_vals, p_vals, m_vals)
# Create a data.frame containing all possible combinations
combinations <-
do.call(expand.grid, values_list)
combinations$w <- combinations$Var1
combinations$p <- combinations$Var2
combinations$m <- combinations$Var3
combinations$Var1 <- NULL
combinations$Var2 <- NULL
combinations$Var3 <- NULL
attr(combinations,"out.attrs") <- NULL
combinations <- combinations %>%
dplyr::filter(w > p) %>%
dplyr::arrange(p)
combinations## w p m
## 1 3 2 1
## 2 5 2 1
## 3 7 2 1
## 4 9 2 1
## 5 11 2 1
## 6 13 2 1
## 7 15 2 1
## 8 17 2 1
## 9 19 2 1
## 10 21 2 1
## 11 23 2 1
## 12 25 2 1
## 13 27 2 1
## 14 29 2 1
## 15 31 2 1
## 16 33 2 1
## 17 35 2 1
## 18 37 2 1
## 19 39 2 1
## 20 41 2 1
## 21 43 2 1
## 22 45 2 1
## 23 47 2 1
## 24 49 2 1
## 25 51 2 1
## 26 53 2 1
## 27 55 2 1
## 28 57 2 1
## 29 59 2 1
## 30 61 2 1
## 31 63 2 1
## 32 65 2 1
## 33 67 2 1
## 34 69 2 1
## 35 71 2 1
## 36 73 2 1
## 37 75 2 1
## 38 77 2 1
## 39 79 2 1
## 40 81 2 1
## 41 83 2 1
## 42 85 2 1
## 43 87 2 1
## 44 89 2 1
## 45 91 2 1
## 46 93 2 1
## 47 95 2 1
## 48 97 2 1
## 49 99 2 1
## 50 101 2 1
## 51 3 2 2
## 52 5 2 2
## 53 7 2 2
## 54 9 2 2
## 55 11 2 2
## 56 13 2 2
## 57 15 2 2
## 58 17 2 2
## 59 19 2 2
## 60 21 2 2
## 61 23 2 2
## 62 25 2 2
## 63 27 2 2
## 64 29 2 2
## 65 31 2 2
## 66 33 2 2
## 67 35 2 2
## 68 37 2 2
## 69 39 2 2
## 70 41 2 2
## 71 43 2 2
## 72 45 2 2
## 73 47 2 2
## 74 49 2 2
## 75 51 2 2
## 76 53 2 2
## 77 55 2 2
## 78 57 2 2
## 79 59 2 2
## 80 61 2 2
## 81 63 2 2
## 82 65 2 2
## 83 67 2 2
## 84 69 2 2
## 85 71 2 2
## 86 73 2 2
## 87 75 2 2
## 88 77 2 2
## 89 79 2 2
## 90 81 2 2
## 91 83 2 2
## 92 85 2 2
## 93 87 2 2
## 94 89 2 2
## 95 91 2 2
## 96 93 2 2
## 97 95 2 2
## 98 97 2 2
## 99 99 2 2
## 100 101 2 2
## 101 5 3 1
## 102 7 3 1
## 103 9 3 1
## 104 11 3 1
## 105 13 3 1
## 106 15 3 1
## 107 17 3 1
## 108 19 3 1
## 109 21 3 1
## 110 23 3 1
## 111 25 3 1
## 112 27 3 1
## 113 29 3 1
## 114 31 3 1
## 115 33 3 1
## 116 35 3 1
## 117 37 3 1
## 118 39 3 1
## 119 41 3 1
## 120 43 3 1
## 121 45 3 1
## 122 47 3 1
## 123 49 3 1
## 124 51 3 1
## 125 53 3 1
## 126 55 3 1
## 127 57 3 1
## 128 59 3 1
## 129 61 3 1
## 130 63 3 1
## 131 65 3 1
## 132 67 3 1
## 133 69 3 1
## 134 71 3 1
## 135 73 3 1
## 136 75 3 1
## 137 77 3 1
## 138 79 3 1
## 139 81 3 1
## 140 83 3 1
## 141 85 3 1
## 142 87 3 1
## 143 89 3 1
## 144 91 3 1
## 145 93 3 1
## 146 95 3 1
## 147 97 3 1
## 148 99 3 1
## 149 101 3 1
## 150 5 3 2
## 151 7 3 2
## 152 9 3 2
## 153 11 3 2
## 154 13 3 2
## 155 15 3 2
## 156 17 3 2
## 157 19 3 2
## 158 21 3 2
## 159 23 3 2
## 160 25 3 2
## 161 27 3 2
## 162 29 3 2
## 163 31 3 2
## 164 33 3 2
## 165 35 3 2
## 166 37 3 2
## 167 39 3 2
## 168 41 3 2
## 169 43 3 2
## 170 45 3 2
## 171 47 3 2
## 172 49 3 2
## 173 51 3 2
## 174 53 3 2
## 175 55 3 2
## 176 57 3 2
## 177 59 3 2
## 178 61 3 2
## 179 63 3 2
## 180 65 3 2
## 181 67 3 2
## 182 69 3 2
## 183 71 3 2
## 184 73 3 2
## 185 75 3 2
## 186 77 3 2
## 187 79 3 2
## 188 81 3 2
## 189 83 3 2
## 190 85 3 2
## 191 87 3 2
## 192 89 3 2
## 193 91 3 2
## 194 93 3 2
## 195 95 3 2
## 196 97 3 2
## 197 99 3 2
## 198 101 3 2
Depois executamos o for loop e estocamos os resultados em uma lista,
por exemplo:
combinations$resultado <- as.list(NA)
for (i in seq_along(combinations$p)) {
# i = 1
sg_combs <- combinations[i,]
spec_sg <- savitzkyGolay(
X = variaveis_nir,
m = sg_combs$m,
p = sg_combs$p,
w = sg_combs$w
)
# Converte spec_sg para data.frame
spec_sg <- as.data.frame(spec_sg)
combinations$resultado[[i]] <- spec_sg
}# dado1 e o conjunto com teus dados e que voce usou para poder derivar
variaveis_nir <- dplyr::select(dado1, starts_with("X")) %>% as.matrix(.)
combinations$acuracia <- NA
for (i in 1:100) {
# i = 1
sg_combs <- combinations[i,]
spec_sg <- savitzkyGolay(
X = variaveis_nir,
m = sg_combs$m,
p = sg_combs$p,
w = sg_combs$w
)
# Converte spec_sg para data.frame
spec_sg <- as.data.frame(spec_sg)
spec_sg$SP1 <- dado1$SP1 # coluna com categoria !!!!
# Divide dataset em 70% treino e 30% de teste
# set.seed(15121890) # se voces quiserem reproduzir os resultados para cada uma das combinacoes, se voces quiserem repetir cada uma das particoes 70/30, usem o SEED; caso contrario
training_samples <- spec_sg$SP1 %>%
createDataPartition(p = 0.7, list = FALSE)
train_data <- spec_sg[training_samples,]
test_data <- spec_sg[-training_samples,]
preproc_param <-
train_data %>%
preProcess(method = c("center", "scale"))
# Transform the data using the estimated parameters
train_transformed <- preproc_param %>% predict(train_data)
test_transformed <- preproc_param %>% predict(test_data)
# Treina o modelo
model <-
lda(SP1 ~ ., data = train_transformed)
# Faz predições
predictions <- model %>% predict(test_transformed)
## Acurácia do modelo
acuracia <-
mean(predictions$class == test_transformed$SP1)
combinations$acuracia[i] <- acuracia
# Remove objetos desnecessarios
rm(
sg_combs,
model,
predictions,
acuracia,
test_transformed,
train_transformed,
preproc_param,
training_samples,
spec_sg,
test_data,
train_data
)
}
combinationscombinations[1, ]Cinco primeiras linhas e cinco primeiras colunas da combinação 1 da
tabela combinations:
combinations$resultado[[1]][1:5, 1:5]Em construção.
ldactrl <-
trainControl(
method = "LOOCV",
verboseIter = TRUE,
savePredictions = "final")lda_loocv <-
train(
dado1_receita, treino,
method = "lda",
trControl = ldactrl,
verbose = TRUE
)# Check CV profile
print(lda_loocv)## Linear Discriminant Analysis
##
## 36 samples
## 1561 predictors
## 2 classes: 'P. aracouchini', 'P. calanense'
##
## Recipe steps:
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 35, 35, 35, 35, 35, 35, ...
## Resampling results:
##
## Accuracy Kappa
## 1 1
# vamos predizer as amostras
lda_preditos <-
predict(lda_loocv, newdata = teste_processado)
confusionMatrix(data = lda_preditos, reference = teste_processado$SP1)## Confusion Matrix and Statistics
##
## Reference
## Prediction P. aracouchini P. calanense
## P. aracouchini 6 0
## P. calanense 0 6
##
## Accuracy : 1
## 95% CI : (0.7354, 1)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : 0.0002441
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0
## Specificity : 1.0
## Pos Pred Value : 1.0
## Neg Pred Value : 1.0
## Prevalence : 0.5
## Detection Rate : 0.5
## Detection Prevalence : 0.5
## Balanced Accuracy : 1.0
##
## 'Positive' Class : P. aracouchini
##
Em construção.
Vamos usar a função trainControl() para estocar os parâmetros de nossa
análise. Em seguida, utilizaremos a função train() para treinar nosso
conjunto de dados treino.
ctrl <-
trainControl(
method = "cv",
savePredictions = "final")plsfitcv <-
train(
dado1_receita, treino,
method = "pls",
metric = "Accuracy",
trControl = ctrl,
tuneLength = 10,
verbose = TRUE
)# Check CV profile
print(plsfitcv)## Partial Least Squares
##
## 36 samples
## 1561 predictors
## 2 classes: 'P. aracouchini', 'P. calanense'
##
## Recipe steps:
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 32, 33, 33, 32, 32, 32, ...
## Resampling results across tuning parameters:
##
## ncomp Accuracy Kappa
## 1 1 1
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 6 1 1
## 7 1 1
## 8 1 1
## 9 1 1
## 10 1 1
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was ncomp = 1.
plot(plsfitcv)
# vamos predizer as amostras
preditos <- predict(plsfitcv, newdata = teste_processado)
confusionMatrix(data = preditos, reference = teste_processado$SP1)## Confusion Matrix and Statistics
##
## Reference
## Prediction P. aracouchini P. calanense
## P. aracouchini 6 0
## P. calanense 0 6
##
## Accuracy : 1
## 95% CI : (0.7354, 1)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : 0.0002441
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0
## Specificity : 1.0
## Pos Pred Value : 1.0
## Neg Pred Value : 1.0
## Prevalence : 0.5
## Detection Rate : 0.5
## Detection Prevalence : 0.5
## Balanced Accuracy : 1.0
##
## 'Positive' Class : P. aracouchini
##
ctrlrepeatcv <-
trainControl(
method = "repeatedcv",
number = 10,
repeats = 10,
savePredictions = "final")plsfitcv100rptcv <-
train(
dado1_receita, treino,
method = "pls",
metric = "Accuracy",
trControl = ctrlrepeatcv,
tuneLength = 10,
verbose = TRUE
)# Check CV profile
print(plsfitcv100rptcv)## Partial Least Squares
##
## 36 samples
## 1561 predictors
## 2 classes: 'P. aracouchini', 'P. calanense'
##
## Recipe steps:
## Resampling: Cross-Validated (10 fold, repeated 10 times)
## Summary of sample sizes: 32, 33, 33, 32, 32, 32, ...
## Resampling results across tuning parameters:
##
## ncomp Accuracy Kappa
## 1 1 1
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 6 1 1
## 7 1 1
## 8 1 1
## 9 1 1
## 10 1 1
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was ncomp = 1.
plot(plsfitcv100rptcv)
# vamos predizer as amostras
preditos100rptcv <- predict(plsfitcv100rptcv, newdata = teste_processado)
confusionMatrix(data = preditos100rptcv, reference = teste_processado$SP1)## Confusion Matrix and Statistics
##
## Reference
## Prediction P. aracouchini P. calanense
## P. aracouchini 6 0
## P. calanense 0 6
##
## Accuracy : 1
## 95% CI : (0.7354, 1)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : 0.0002441
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0
## Specificity : 1.0
## Pos Pred Value : 1.0
## Neg Pred Value : 1.0
## Prevalence : 0.5
## Detection Rate : 0.5
## Detection Prevalence : 0.5
## Balanced Accuracy : 1.0
##
## 'Positive' Class : P. aracouchini
##
Bache, Stefan Milton, e Hadley Wickham. 2020. magrittr: A Forward-Pipe Operator for R. https://CRAN.R-project.org/package=magrittr.
Kuhn, Max. 2020. caret: Classification and Regression Training. https://github.com/topepo/caret/.
Kuhn, Max, Fanny Chow, e Hadley Wickham. 2020. rsample: General Resampling Infrastructure. https://CRAN.R-project.org/package=rsample.
Kuhn, Max, e Hadley Wickham. 2020a. recipes: Preprocessing Tools to Create Design Matrices. https://CRAN.R-project.org/package=recipes.
———. 2020b. tidymodels: Easily Install and Load the Tidymodels Packages. https://CRAN.R-project.org/package=tidymodels.
Perdiz, R. O. 2020. “Delimitação específica e filogeografia do complexo Protium aracouchini (Aubl.) Marchand (Burseraceae)”. Tese de doutorado, Manaus, Amazonas, Brasil: Programa de pós-graduação em Ciências Biológicas (Botânica), Instituto Nacional de Pesquisas da Amazônia. https://repositorio.inpa.gov.br/handle/1/36948.
Perdiz, Ricardo Oliveira. 2021. NIRtools: Tools to deal with near infrared (NIR) spectrocospy data. https://github.com/ricoperdiz/NIRtools.
Stevens, Antoine, Leonardo Ramirez-Lopez, e Guillaume Hans. 2020. prospectr: Miscellaneous Functions for Processing and Sample Selection of Spectroscopic Data. https://github.com/l-ramirez-lopez/prospectr.
Wickham, Hadley, Romain François, Lionel Henry, e Kirill Müller. 2021. dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.