0

Estou começando a usar o R e tenho um banco de dados que precisa ser dividido em sub-bases para que se possa determinar os coeficientes de uma regressão não linear em cada uma delas.

O critério de separação é definido pelo argumento da variável 'amo'. Para dividir os dados usei a função split.data.frame e foram criados os 87 objetos que precisava em uma lista, mas agora não sei como aplicar a função nls em todos eles, sem precisar separar manualmente cada um.

O que fiz, foi:

dados #banco de dados com 1539 linhas e 18 colunas
subbases <- split.data.frame(dados, f = dados$amo)

O que pensei em fazer foi:

Primeiramente separar cada um dos objetos

amo1 <- subbases$`1`
amo2 <- subbases$`2`
amo3 <- subbases$`3`
... 

Depois aplicar a função nls em cada um deles:

regressao1 <- nls(variavel18 ~ k1*(variavel16^k2)*(variavel17^k3), 
                  data = amo1, 
                  start= list(k1=2000, k2=0.4, k3=-0.2))

regressao2 <- nls(variavel18 ~ k1*(variavel16^k2)*(variavel17^k3), 
                  data = amo2, 
                  start= list(k1=2000, k2=0.4, k3=-0.2))

...

O problema é que são 87 sub-bases e acredito que deve haver alguma forma mais elegante de fazer esse procedimento.

2
  • Olá @StephannyEgito, por favor, disponibilize uma amostra da base.
    – itamar
    Commented 6/02/2023 às 20:55
  • Oi @AugustoVasques, muito obrigada pela disponibilidade. Sou completamente iniciante nessa área. Não consegui anexar o arquivo aqui, então fiz um link para o GitHub. Espero que tenha feito certo e você consiga visualizar. link. Commented 7/02/2023 às 11:26

1 Resposta 1

3

A forma como eu gosto de fazer isso é criando uma lista de modelos no próprio objeto para fácil controle e organização com o comando nest do pacote tidyr.

Primeiro eu crio uma função que descreva o modelo padrão que irei aplicar. Na lista que criamos através do nest usamos a função tidyr::map para aplicar o modelo em cada elemento da lista separadamente.

library(dplyr)
library(tidyr)
library(purrr)

urlfile <- 'https://raw.githubusercontent.com/StephannyEgito/RNA/main/Amostra_Git_csv.csv'
dt <-read.csv2(urlfile)

modelo_padrao <- function(dados)
  nls(MR ~ k1*(confinante^k2)*(desvio^k3), 
      data = dados, 
      start= list(k1=2000, k2=0.4, k3=-0.2))

dt_new <- 
dt %>% 
  nest(data = -amo) %>% 
  mutate(teste = map(data, modelo_padrao),
         tidy_mod = map(teste, broom::tidy)) 

Dessa forma podemos acessar tanto os teste individualmente em cada elemento da lista

head(dt_new$teste, 3)  
#> [[1]]
#> Nonlinear regression model
#>   model: MR ~ k1 * (confinante^k2) * (desvio^k3)
#>    data: dados
#>        k1        k2        k3 
#> 1525.8538    0.5458   -0.1743 
#>  residual sum-of-squares: 32784
#> 
#> Number of iterations to convergence: 6 
#> Achieved convergence tolerance: 7.96e-06
#> 
#> [[2]]
#> Nonlinear regression model
#>   model: MR ~ k1 * (confinante^k2) * (desvio^k3)
#>    data: dados
#>        k1        k2        k3 
#> 8.115e+03 7.308e-01 2.655e-02 
#>  residual sum-of-squares: 408876
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 1.048e-06
#> 
#> [[3]]
#> Nonlinear regression model
#>   model: MR ~ k1 * (confinante^k2) * (desvio^k3)
#>    data: dados
#>        k1        k2        k3 
#> 8863.3252    0.6966    0.0106 
#>  residual sum-of-squares: 64422
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 4.794e-07

Quanto em um formato limpo para comparação direta através de broom::tidy e tidyr::unnest

unnest(dt_new, tidy_mod)
#> # A tibble: 30 × 8
#>      amo data               teste  term   estimate std.error statistic  p.value
#>    <int> <list>             <list> <chr>     <dbl>     <dbl>     <dbl>    <dbl>
#>  1     1 <tibble [18 × 18]> <nls>  k1    1526.      156.         9.76  6.84e- 8
#>  2     1 <tibble [18 × 18]> <nls>  k2       0.546     0.0605     9.02  1.90e- 7
#>  3     1 <tibble [18 × 18]> <nls>  k3      -0.174     0.0471    -3.70  2.15e- 3
#>  4     2 <tibble [18 × 18]> <nls>  k1    8115.     1431.         5.67  4.45e- 5
#>  5     2 <tibble [18 × 18]> <nls>  k2       0.731     0.102      7.15  3.31e- 6
#>  6     2 <tibble [18 × 18]> <nls>  k3       0.0265    0.0732     0.363 7.22e- 1
#>  7     3 <tibble [18 × 18]> <nls>  k1    8863.      500.        17.7   1.78e-11
#>  8     3 <tibble [18 × 18]> <nls>  k2       0.697     0.0328    21.3   1.30e-12
#>  9     3 <tibble [18 × 18]> <nls>  k3       0.0106    0.0238     0.446 6.62e- 1
#> 10     4 <tibble [18 × 18]> <nls>  k1    1313.       40.8       32.2   2.90e-15
#> # … with 20 more rows

Created on 2023-02-07 with reprex v2.0.2

Edit

Além da função tidy podemos utilizar broom::glance() para obter estatística relacionadas somente ao modelo (loglik, AIC, BIC, etc...)

dt_new %>% 
  mutate(glance_mod = map(teste, broom::glance)) %>% 
  unnest(glance_mod)
#> # A tibble: 10 × 13
#>      amo data     teste tidy_mod sigma isConv  finTol logLik   AIC   BIC devia…¹
#>    <int> <list>   <lis> <list>   <dbl> <lgl>    <dbl>  <dbl> <dbl> <dbl>   <dbl>
#>  1     1 <tibble> <nls> <tibble>  46.8 TRUE   7.96e-6  -93.1  194.  198.  32784.
#>  2     2 <tibble> <nls> <tibble> 165.  TRUE   1.05e-6 -116.   240.  243. 408876.
#>  3     3 <tibble> <nls> <tibble>  65.5 TRUE   4.79e-7  -99.2  206.  210.  64422.
#>  4     4 <tibble> <nls> <tibble>  21.9 TRUE   5.01e-6  -79.5  167.  171.   7207.
#>  5     5 <tibble> <nls> <tibble>  47.3 TRUE   7.93e-7  -93.3  195.  198.  33565.
#>  6     6 <tibble> <nls> <tibble>  56.2 TRUE   1.94e-6  -96.4  201.  204.  47372.
#>  7     7 <tibble> <nls> <tibble>  71.1 TRUE   2.48e-6 -101.   209.  213.  75778.
#>  8     8 <tibble> <nls> <tibble>  70.5 TRUE   1.83e-7 -100.   209.  213.  74469.
#>  9     9 <tibble> <nls> <tibble>  70.4 TRUE   2.24e-6 -100.   209.  212.  74260.
#> 10    10 <tibble> <nls> <tibble>  42.4 TRUE   2.90e-6  -91.4  191.  194.  26981.
#> # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
#> #   variable name ¹​deviance
3
  • Muitíssimo obrigada pelo socorro que me prestou. Ainda não conhecia nenhuma das funções que você mencionou. Testei e deu tudo certo. Gostaria de tirar mais um dúvida, após realizar todo esse processo, é possível extrair os valores dos resíduos, para que eu possa calcular o erro residual e o R² de cada uma das amostras? Fazendo da forma arcaica que comecei, estava usando o 'fitted' e depois o 'cor' para cada amostra. Commented 7/02/2023 às 19:26
  • @StephannyEgito Você pode criar uma função para isso, da mesma forma que foi criado a função modelo_padrao e depois aplicar ela dentro do map para retornar o R² de cada elemento da lista. Não coloco um exemplo pois o cálculo e interpretação do R² é controverso para modelos não lineares. De forma alternativa você pode usar a função broom::glance para retornar estatística relacionadas diretamente ao modelo, vou editar a resposta para exemplificar isso.
    – M Aurélio
    Commented 7/02/2023 às 19:35
  • @StephannyEgito se a resposta foi satisfatória e atendeu a sua dúvida inicial gostaria que a marcasse como aceita.
    – M Aurélio
    Commented 7/02/2023 às 19:35

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .