Análise de Dados COVID-19 em Vargem Grande Paulista SP 2020-2023

Vargem Grande Paulista Preparação do R para as questões desta atividade setwdUsersCarolineLibraryMobileDocumentscomappleCloudDocsUnifesp8ºtermoAnalise de dadosBasesHISTPAINELCOVIDBR06out2023 init Systime dados20201HISTPAINELCOVIDBR2020Parte106out2023csv dados1readcsv2dados20201header TRUE sep dados20202HISTPAINELCOVIDBR2020Parte206out2023csv dados2readcsv2dados20202header TRUE sep dados20211HISTPAINELCOVIDBR2021Parte106out2023csv dados3readcsv2dados20211header TRUE sep dados20212HISTPAINELCOVIDBR2021Parte206out2023csv dados4readcsv2dados20212header TRUE sep dados20221HISTPAINELCOVIDBR2022Parte106out2023csv dados5readcsv2dados20221header TRUE sep dados20222HISTPAINELCOVIDBR2022Parte206out2023csv dados6readcsv2dados20222header TRUE sep dados20231HISTPAINELCOVIDBR2023Parte106out2023csv dados7readcsv2dados20231header TRUE sep dados20232HISTPAINELCOVIDBR2023Parte206out2023csv dados8readcsv2dados20232header TRUE sep Final Systime init Aqqui eu vou descobrir o tempo que ficou para rodar as oito bases no meu computador Final Time difference of 1468633 secs dadosrbinddados1dados2dados3dados4dados5dados6dados7dados8 objectsize dados1000000 6669 bytes Separando a cidade escolhida dadosVargemdadoswhichdadosestadoSPdadosmunicipioVargem Grande Paulista dadosdata 20200328dadosdata20230905 head dadosVargem regiao estado municipio coduf codmun codRegiaoSaude nomeRegiaoSaude 501758 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 501759 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 501760 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 501761 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 501762 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 501763 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS data semanaEpi populacaoTCU2019 casosAcumulado casosNovos obitosAcumulado 501758 20200328 13 52597 2 2 1 501759 20200329 14 52597 2 0 1 501760 20200330 14 52597 2 0 1 501761 20200331 14 52597 2 0 1 501762 20200401 14 52597 2 0 1 501763 20200402 14 52597 3 1 1 obitosNovos Recuperadosnovos emAcompanhamentoNovos interiormetropolitana 501758 1 NA NA 1 501759 0 NA NA 1 501760 0 NA NA 1 501761 0 NA NA 1 501762 0 NA NA 1 501763 0 NA NA 1 tail dadosVargem regiao estado municipio coduf codmun codRegiaoSaude nomeRegiaoSaude 7079671 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 7079672 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 7079673 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 7079674 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 7079675 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS 7079676 Sudeste SP Vargem Grande Paulista 35 355645 35013 MANANCIAIS data semanaEpi populacaoTCU2019 casosAcumulado casosNovos obitosAcumulado 7079671 20230831 35 52597 7137 0 140 7079672 20230901 35 52597 7137 0 140 7079673 20230902 35 52597 7153 16 140 7079674 20230903 36 52597 7153 0 140 7079675 20230904 36 52597 7153 0 140 7079676 20230905 36 52597 7153 0 140 obitosNovos Recuperadosnovos emAcompanhamentoNovos interiormetropolitana 7079671 0 NA NA 1 7079672 0 NA NA 1 7079673 0 NA NA 1 7079674 0 NA NA 1 7079675 0 NA NA 1 7079676 0 NA NA 1 dim dadosVargem 1 1257 17 1 não pontua Com dados obtidos de httpscovidsaudegovbr escolha uma cidade do estado de São Paulo diferente da cidade de São Paulo a partir de 28032020 até exatamente e inclusive 05092023 Verifique se os dados de casos acumulados da cidade escolhida são diferentes de zero a partir de 21042020 Em seguida utilizando do R importe os dados e faça o data mining necessário para obter as séries de dados diários de todas as variáveis disponíveis somente para a cidade escolhida I Apresentar os códigos do data mining realizado através DO CONSOLE de modo a se verificar que foram implementados II Apresentar o head e o tail do data frame gerado para as estimativas com todas as colunas do data frame utilizado na modelagem III Apresentar o resultado DO CONSOLE do R da aplicação da função dim sobre o data frame utilizado Definindo 𝑌𝑡 como o número de óbitos acumulados coluna casosAcumulado causados por COVID19 no dia t desenvolva as seguintes questões 1 Considere o modelo de crescimento logístico yt α 1 e β γ t e que yt N y tσ Separar títulos names dadosVargem 1 regiao estado municipio 4 coduf codmun codRegiaoSaude 7 nomeRegiaoSaude data semanaEpi 10 populacaoTCU2019 casosAcumulado casosNovos 13 obitosAcumulado obitosNovos Recuperadosnovos 16 emAcompanhamentoNovos interiormetropolitana Separando os óbitos da cidade Em Yt Yt dadosVargemcasosAcumulado2lengthdadosVargemcasosAcumulado head Yt 1 2 2 2 2 3 4 tail Yt 1 7137 7137 7153 7153 7153 7153 21 Escreva textualmente a expressão matemática para a função log da Verossimilhança l θ para a estimação dos parâmetros do modelo Temos que I l θ y n ln 2π 2 n ln σ 1 2 σ 2 i1 n u i 2 y i y i u i então II u i y i y i III y i α 1 e β γ t Juntando I II e III l θ y n ln 2π 2 n ln σ 1 2 σ 2 i1 n y i α 1 e β γ t 2 22 Escreva textualmente a expressão matemática para o vetor de parâmetros θ da função l θ para esse modelo Quais seriam bons chutes iniciais para a estimativa dos parâmetros 𝛽 e 𝜎 Apresente os resultados dos cálculos para esses chutes através do console do R Os parâmetros estão na expressão de estimação para a função log θ θ α β γ σ Os bons chutes iniciais seriam σ 2 i1 n u 2 n No R Tamanho da base de Vargem Grande T length Yt T 1 1256 Criação de vetor de 1 ao último número do tamanho da base t seq 1T sigma var Yt T1T05 sigma 1 00004011782 β ln α min Yt 1 dado que α min Yt então para o melhor chute α max Yt No R ac max Yt ac 1 7153 bc log ac min Yt 1 bc 1 818186 23 Escreva implemente a função l θ no R para o modelo logístico Apresente o resultado através do console do R Implementando l n θ inicialmente crio o vetor dos parâmetros chute ac max Yt bc log ac min Yt 1 gc 0 sigma var Yt T1T05 sigma 1 00004011782 theta c acbcgcsigma theta 1 7153000e03 8181860e00 0000000e00 4011782e04 Equação do modelo u i y i α 1 e β γ t l θ y n ln 2π 2 n ln σ 1 2 σ 2 i1 n y i α 1 e β γ t 2 Criação da função do modelo logístico funcML function theta defino os coeficientes a theta 1 alpha b theta 2 beta g theta 3 gama sigma theta 4 cirase o vetor de residuos u Yt a1expbgt escrevo minha função de logartimo da verossimilhança lnL com uma distribuição normal lnL 05Tlog2piTlog abs sigma05sumu2sigma2 return lnL na frente para maximizála 24 05 Considerando diferentes chutes iniciais para α partindo de max yt acrescentando 5 10 15 e com o uso da função nlm com opção printlevel 1 estime os parâmetros da função lθ para o período de 21042020 a 01082023 Busque convergência das estimativas calibrando steptol gradtol e iterlim se necessário Para cada um dos chutes iniciais para α apresente os valores da função do log da verossimilhança e dos parâmetros estimados Apresente os códigos e resultados através do console do R Somente considerada se itens 23 e 24 estiverem corretos Usando a função NLM para todas as bases lnLmax nlm funcML theta steptol 1e8 gradtol 1e7 printlevel 1 ndigit 12 iterlim 100000 iteration 0 Step 1 0 0 0 0 Parameter 1 7153000e03 8181860e00 0000000e00 4011782e04 Function Value 1 7196218e16 Gradient 1 7630551e09 5456610e13 4824198e16 3574173e20 iteration 2 Parameter 1 7153000e03 8181860e00 1423010e04 1054687e00 Function Value 1 10410880324 Gradient 1 1252752e03 8958118e06 8009612e09 1974209e10 Successive iterates within tolerance Current iterate is probably solution Separando as datas solicitadas no t ítulo do exercício em ystar date asDate dadosVargemdata ystar Yt whichdate20200421date20230801 tstar 1length ystar Tstar length ystar Tstar 1 1198 ac1 max ystar bc1 logac1min ystar 1 gc1 0 sigma1 var ystar Tstar1 Tstar 05 sigma1 1 00004119329 theta1 cac1bc1gc1sigma1 theta1 1 7123000e03 6567094e00 0000000e00 4119329e04 funcMLstar function theta1 Definindo coeficientes a theta11 alpha b theta12 beta g theta13 gamma sigma var ystar Tstar1 Tstar 05 Criase vetor de residuos u ystar a1expbg tstar funcao de log da verossimilhança lnL 05length ystar log2pi 05length ystar logabssigma05sumu2sigma2 minlnL max lnL return lnL Usando a função NLM para ystar e adaptando o steptol para o resultado da base ser tolerável lnLmax nlm funcMLstar theta1 steptol 1e 12 gradtol 1e 9 printlevel 1 ndigit 12 iterlim 100000 iteration 0 Step 1 0 0 0 0 Parameter 1 7123000e03 6567094e00 0000000e00 4119329e04 Function Value 1 7582765 Gradient 1 9862676e04 7015348e00 5891715e03 0000000e00 iteration 9 Parameter 1 7122998e03 3355158e00 5337867e03 4119329e04 Function Value 1 5 786 126 Gradient 1 001335193 2357811708 91348117076 000000000 Successive iterates within tolerance Current iterate is probably solution Por fim irei comparar os parâmetros chutes do ystar aos parâmetros maximizados ah lnLmaxestimate 1 bh lnLmaxestimate 2 gh lnLmaxestimate 3 sigmah lnLmaxestimate 4 thetah c ahbhghsigmah theta1 chute 1 7123000e03 6567094e00 0000000e00 4119329e04 thetah 1 7122998e03 3355158e00 5337867e03 4119329e04 yh ah 1exp bhgh tstar plot ystartstartype p lines yhtstar typel lty 2colpurple Para o alfa com acréscimo de 5 ac2 max ystar 1 05 bc2 logac2min ystar 1 gc2 0 sigma1 como a fórmula do sigma não muda este permanecera o mesmo 1 00004119329 theta2 cac2bc2gc2sigma1 theta2 1 7479150e03 6615951e00 0000000e00 4119329e04 funcMLstar2 function theta2 Definindo coeficientes a theta21 alpha b theta22 beta g theta23 gamma sigma var ystar Tstar1 Tstar 05 u ystar a1expbg tstar fun çã o de log da verossimilhança lnL 05length ystar log2pi 05length ystar logabssigma05sumu2sigma2 return lnL lnLmax2 nlm funcMLstar2 theta2 steptol 1e4 gradtol 1e7 printlevel 1 ndigit 12 iterlim 100000 iteration 0 Step 1 0 0 0 0 Parameter 1 7479150e03 6615951e00 0000000e00 4119329e04 Function Value 1 7582765 Gradient 1 9393025e04 7015817e00 5892110e03 0000000e00 iteration 6 Parameter 1 7479150e03 6615882e00 9531821e03 4119329e04 Function Value 1 5851183 Gradient 1 002212979 2628031795 98129152957 000000000 Successive iterates within tolerance Current iterate is probably solution Comparação os parâmetros chutes do ystar aos parâmetros maximizados ah2 lnLmax2estimate1 bh2 lnLmax2estimate2 gh2 lnLmax2estimate3 sigmah2 lnLmax2estimate4 thetah2 cah2bh2gh2sigmah2 lnLmax2 lnLmax2estimate theta2 1 7479150e03 6615951e00 0000000e00 4119329e04 thetah2 1 7479150e03 6615882e00 9531821e03 4119329e04 Para o alfa com acréscimo de 10 ac3 max ystar 11 bc3 logac3min ystar 1 gc3 0 sigma1 como a fórmula do sigma não muda este permanecera o mesmo 1 00004119329 theta3 cac3bc3gc3sigma1 theta3 1 7835300e03 6662532e00 0000000e00 4119329e04 funcMLstar3 function theta3 Definindo coeficientes a theta31 alpha b theta32 beta g theta33 gamma sigma var ystar Tstar1 Tstar 05 Criase vetor de residuos u ystar a1expbg tstar funcao de log da verossimilhança lnL 05length ystar log2pi 05length ystar logabssigma05sumu2sigma2 minlnL max lnL return lnL lnLmax3 nlm funcMLstar3 theta3 steptol 1e 12 gradtol 1e 10 printlevel 1ndigit 12 iterlim 100 iteration 0 Step 1 0 0 0 0 Parameter 1 7835300e03 6662532e00 0000000e00 4119329e04 Function Value 1 7582765 Gradient 1 8966070e04 7016245e00 5892469e03 0000000e00 iteration 27 Parameter 1 7835296e03 3209625e00 4727869e03 4119329e04 Function Value 1 5781125 Gradient 1 0003452703 0047567883 38287837015 0000000000 Successive iterates within tolerance Current iterate is probably solution Comparação os parâmetros chutes do ystar aos parâmetros maximizados ah3 lnLmax 3 estimate1 bh3 lnLmax 3 estimate2 gh3 lnLmax 3 estimate3 sigmah3 lnLmax3estimate4 thetah3 cah3bh3gh3sigmah3 lnLmax3v lnLmax3minimum theta3 1 7835300e03 6662532e00 0000000e00 4119329e04 thetah3 1 7835300e03 6662532e00 0000000e00 4119329e04 lnLmax3v 1 5 78 1925 Para o alfa com acréscimo de 1 5 ac4 max ystar 115 bc4 logac4min ystar 1 gc4 0 sigma1 como a fórmula do sigma não muda este permanecera o mesmo 1 00004119329 theta4 cac4bc4gc4sigma1 theta4 1 8191450e03 6707040e00 0000000e00 4119329e04 funcMLstar4 function theta4 Definindo coeficientes a theta41 alpha b theta42 beta g theta43 gamma sigma var ystar Tstar1 Tstar 05 Criase vetor de residuos u ystar a1expbg tstar funcao de log da verossimilhança lnL 05length ystar log2pi 05length ystar logabssigma05sumu2sigma2 minlnL max lnL return lnL lnLmax4 nlm funcMLstar4 theta4 steptol 1e4 gradtol 1e 12 printlevel 1 ndigit 12 iterlim 100000 iteration 0 Step 1 0 0 0 0 Parameter 1 8191450e03 6707040e00 0000000e00 4119329e04 Function Value 1 7582765 Gradient 1 8576240e04 7016634e00 5892796e03 0000000e00 iteration 8 Parameter 1 8191450e03 6706971e00 9089734e03 4119329e04 Function Value 1 5875018 Gradient 1 003664179 2942906023 124355718799 000000000 Successive iterates within tolerance Current iterate is probably solution Comparação os parâmetros chutes do ystar aos parâmetros maximizados ah4 lnLmax4estimate1 bh4 lnLmax4estimate2 gh4 lnLmax4estimate3 sigmah4 lnLmax4estimate4 thetah4 cah4bh4gh4sigmah4 lnLmax4v lnLmax4minimum theta4 1 8191450e03 6707040e00 0000000e00 4119329e04 thetah4 1 8191450e03 6706971e00 9089734e03 4119329e04 lnLmax4v 1 5875018 25 025 Qual dos chutes para α maximizou lθ Você considera que o algoritmo de otimização foi bemsucedido nas estimativas dos parâmetros para cada um dos chutes de α utilizados Justifique devidamente sua resposta Somente considerada se itens 23 24 e 25 estiverem corretos Foi criada a seguinte tabela no R com os chutes do alfa em relação ao fa t or da l θ Notase que entre os chutes que de com alfa maximizado ao alfa mais 15 que o alfa maximizado apresentou um valor de função menor ou seja entre estes o alfa que mais maximizou a l θ Chutes para o Y max Valor da função maximizada Max y 5786 Max y 105 5851 Max y 110 5781 Max y 115 5875 Após isso foi criado um loop para tentar encontrar se há outro alfa que maximizaria melhor esta l θ Segue o código utilizado acloop max ystar seq092005 bcloop c funclog c for i in 1length acloop bcloop i log acloop i min ystar 1 gcloop 0 chutesloop c acloop i bcloop i gcloop sigma1 cat Valor para i acloop i acloop 1n funclog i nlm funcMLstar chutesloopsteptol 1e4 gradtol 1e8 printlevel 1 iterlim 100000 ndigit 12 fscale 10000 minimum Criar a tabela dados dataframe Valoresalfa max ystar seq 092005 Valordaagravo seq 092005 Valordafunção funclog 123 dados Foi criada uma tabela com esses valores podese concluir que o 135 do alfa seria o valor que mais maximizaria esta função Valores de alfa Aumento Valor da Função 641070 090 7250350 676685 095 7476481 712300 100 5843334 747915 105 5851183 783530 110 5861925 819145 115 5875018 854760 120 5889885 890375 125 5906036 925990 130 5922992 961605 135 5781840 997220 140 5782651 1032835 145 5783532 1068450 150 5784411 1104065 155 5785285 1139680 160 5786139 1175295 165 5786965 1210910 170 5787762 1246525 175 6063381 1282140 180 6075401 1317755 185 6086678 1353370 190 5790619 1388985 195 5791252 1424600 200 5791856 Mesmo que o 135 de alfa sinalize anteriormente que seria o valor que maximizaria esta função temos um problema ao rodar a nlm com este steptol Last global step failed to locate a point lower than x Either x is an approximate local minimum of the function the function is too nonlinear for this algorithm or steptol is too large Por isso decidiuse seguir com o alfa a maximizado 100 como o valor com maior otimização de resultados 26 Escreva textualmente a expressão matemática para o modelo estimado utilizando as estimativas para os parâmetros a partir do chute de α que maximizou lθ Somente considerada se itens 23 24 e 25 estiverem corretos l θ y 1198 ln 2π 2 1198 ln σ 1 2 0000411933 2 i1 1198 y i 7123 1 e 3355150005337 2 27 Para esse modelo estime os errospadrão dos estimadores dos parâmetros sua matriz de variânciacovariância dos estimadores θ e seus respectivos pvalores para os testes de significância individual de cada parâmetro Apresente os códigos e resultados através do console do R Encontrar a hessianamatriz identidade library numDeriv hess hessian funcMLstarlnLmaxestimate hess 1 2 3 4 1 5804091e05 1190770e01 9486028e01 1153139e14 2 1190770e01 3571818e02 2423255e05 1795240e10 3 9486028e01 2423255e05 1831409e08 1502994e07 4 1153139e14 1795240e10 1502994e07 0000000e00 Matriz Var Cov Foi utilizado o valor de regularização para auxiliar no calculo da inversa lambda 1e 7 Valor de regularização ajuste conforme necessário hessreg hess lambda 4 Tamanho da diagonal de hess varcovarML solve hessreg varcovarML 1 2 3 4 1 1286606e05 2266775e01 9663462e02 1286832e05 2 2266775e01 3135686e02 5323132e05 2269906e01 3 9663462e02 5323132e05 1259472e07 9668775e02 4 1286832e05 2269906e01 9668775e02 2628706e06 Erro padrão varML abs diag varcovarML varML 1 1286606e05 3135686e02 1259472e07 2628706e06 epML varML051length varML 1 epML 1 3586929e02 1770787e01 3548904e04 4Teste de significancia z thetah 13epML z 1 2184402 1812542 1332205 pvalor 2 1 pnorm abs z pvalor 1 0 0 0 28 025 Interprete os pvalores calculados no item anterior Somente considerada se itens 23 24 25 e 28 estiverem corretos Ho pvalor5Parâmetros estatisticamentes iguais a zero H1 pvar5Parâmetros estatisticamentes diferentes de zero pvalor α 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada alfa ser menor que zero é de aproximadamente 0 pvalor β 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada beta ser menor que zero é de aproximadamente 0 pvalor γ 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada gama ser menor que zero é de aproximadamente 0 pvalor σ 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada sigma ser menor que zero é de aproximadamente 0 29 Utilizando das estimativas para os parâmetros a partir do chute de α que maximizou lθ obtenha as estimativas para ŷt para o período de 21042020 a 01082023 Em seguida faça o plot da curva de crescimento logístico ajustada aos dados Apresente os resultados dos códigos através do console do R Somente considerada se itens 23 24 e 25 estiverem corretos Yh ah 1 exp bh gh tstar Yh date2asDatedadosVargemdatawhichdadosVargemdata20200421dadosVargemdata20230801 plot ystardate2xaxtn ylab Casos acumulados xlab Data main Casos Acumulados COVID Vargem Grande type p linesYhdate2colpurple 210 Calcule e apresente através do console do R os resultados para o MSE RMSE MAE e MAPE das estimativas de ŷt geradas no item 27 Somente considerada se itens 23 24 25 26 27 e 210 estiverem corretos residlog ystarYh MSE sum residlog 2 Tstar MSE 1 1105966 RMSE sqrtMSE RMSE 1 3325607 MAE sumabs residlog Tstar MAE 1 2767209 MAPE sumabs residlog ystar Tstar MAPEmean absolute percentage error 1 06813478 T em uma diferença média do valor observado para o valor calculado de 413 entre observado e simulado 211 Calcule os valores previstos de yt para 02082023 até 05092023 Somente considerada se itens 23 24 25 26 27 e 210 estiverem corretos Dado que a função de previsão de yt é yt α 1 e β γ t Até dia 509 são 34 dias e teríamos que somar o máximo de dias in sample mais os 34 dias e calcular a sequência Yhx ah 1 exp bh gh max tstar seq 135 Yhx Criar tabela Defina a data de início e a data de término datainicio asDate 20230802 datafim asDate 20230905 Cri ando a sequência de datas Datasolicitada seq datainicio datafim by days Formate as datas para o formato desejado dd mm yyyy DataPrevista format Datasolicitada format dmY DataPrevista Crie um dataframe combinando os dois vetores tabela dataframe Data DataPrevista Yhx Yhx Data Yhx 1 02082023 6799356 2 03082023 6801001 3 04082023 6802638 4 05082023 6804267 5 06082023 6805889 6 07082023 6807502 7 08082023 6809108 8 09082023 6810705 9 10082023 6812295 10 11082023 6813878 11 12082023 6815452 12 13082023 6817019 13 14082023 6818579 14 15082023 6820130 15 16082023 6821674 16 17082023 6823211 17 18082023 6824740 18 19082023 6826262 19 20082023 6827776 20 21082023 6829283 21 22082023 6830782 22 23082023 6832275 23 24082023 6833759 24 25082023 6835237 25 26082023 6836707 26 27082023 6838171 27 28082023 6839627 28 29082023 6841075 29 30082023 6842517 30 31082023 6843952 31 01092023 6845380 32 02092023 6846800 33 03092023 6848214 34 04092023 6849621 212 Construa um gráfico comparando os valores observados e previstos de yt no item acima yobs Yt whichdate20230802date20230905 yobs Yt whichdate20230802date20230904 yobs 1 7123 7123 7126 7126 7126 7126 7126 7126 7126 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 21 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7153 7153 7153 7153 yobs yobs isna yobs yobs 1 7123 7123 7126 7126 7126 7126 7126 7126 7126 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 21 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7153 7153 7153 7153 Ajustar um modelo de regressão linear tendencia lm yobs Yhx tendencia plot yobsYhxylab Casos observado xlab Casos Previstos main Obsv X Prevs typep lines yobsYhx abline tendencia col red 213 Calcule o MSE RMSE MAE e o MAPE das previsões de yt de 02082023 até 05092023 Interprete o MAPE dessas previsões residlog yobsYhx MSE sum residlog 2 Thx MSE 1 2204925 RMSE sqrtMSE RMSE 1 4695663 MAE sumabs residlog Thx MAE 1 4662351 MAPE sumabs residlog Yhx Thx MAPE 1 000657649 T em uma diferença média do valor observado para o valor calculado de 000657 entre observado e simulado 2 35 pontos Considere o modelo de crescimento de Gompertz α e e βγt ut e que ytN yt σ 2 31 025 Escreva textualmente a expressão matemática para a função log da verossimilhança lθ para a estimação dos parâmetros do modelo l θ y n ln 2π 2 n ln σ 1 2 σ 2 i1 n y i α e e βγt 2 32 025 Escreva textualmente a expressão matemática para o vetor de parâmetros θ da função lθ para esse modelo Quais seriam bons chutes iniciais para a estimativa dos parâmetros β e σ θ α β γ σ Sigma é o desvio padrão desta distribuição sigmag var Yt T 1 T 05 sigmag ag max Yt ag 1 7153 bg loglog ag logmin Yt bg 1 2101954 33 Não pontua Escreva implemente a função lθ no R para o modelo logístico Apresente o resultado do console do R funcMLgomp functiontheta Defino meus coeficiente a theta 1 a alpha b theta 2 b beta g theta3 g gamma sigma var Yt length Yt 1length Yt 05 Criase um vetor de residuos u Yt aexpexpbg Yt Escrevo a fun o de log da verossimilhanca lnL com uma distribuio normal lnL 05length Yt log2pi05length Yt logabs sigmag 05sumu2sigmag2 na frente da funcao para maximizala no algoritmo de otimizacao Agora calculo o valor da funcao return lnL lnLmaxlog nlm funcMLgomp chutes steptol 1e4 gradtol 1e6 printlevel 1 iterlim 100000 ndigit 12 fscale 1 iteration 0 Step 1 0 0 0 Parameter 1 7153000000 2101954 0000000 Function Value 1 7929873 Gradient 1 1976538e04 1156810e01 6236320e04 iteration 105 Parameter 1 7153000e03 2113364e00 7297128e04 Function Value 1 609394 Gradient 1 5090611e03 2425766e01 3167652e04 Successive iterates within tolerance Current iterate is probably solution lnLmaxlogestimate 1 7153000e03 2113364e00 7297128e04 34Considerando diferentes chutes iniciais para α partindo de max yt acrescentando 5 10 15 e com o uso da função nlm com opção printlevel 1 estime os parâmetros da função lθ para o período de 21042020 a 01082023 Busque convergência das estimativas calibrando steptol gradtol e iterlim se necessário Para cada um dos chutes iniciais para α apresente os valores da função do log da verossimilhança e dos parâmetros estimados Apresente os códigos e resultados através do console do R Somente considerada se itens 33 e 34 estiverem corretos O valor da função para alfa max ystar lnLmaxlogestimate 1 7153000e03 2113364e00 7297128e04 sigmag1 var ystar Tstar 1 Tstar 05 sigmag1 1 242758 ag1max ystar ag1 1 7123 bg1loglogag1logmin ystar bg1 1 1882285 gg10 theta1cag1bg1gg1 funcMLgomp1functiontheta Defino meus coeficiente a theta 1 a alpha b theta 2 b beta g theta3 g gamma sigma var ystar length ystar 1length ystar 05 Criase um vetor de residuos u ystar aexpexpbg tstar lnL 05lengthystarlog2pi05lengthystarlogabssigmag1 05sumu2sigmag12 na frente da função para maximizala no algoritmo de otimização Agora calculo o valor da funcao return lnL lnLmaxlog1 nlm funcMLgomp1 theta1 steptol 1e4 gradtol 1e6 printlevel 1 iterlim 100000 ndigit 12 fscale 1 iteration 0 Step 1 0 0 0 Parameter 1 7123000000 1882285 0000000 Function Value 1 7582765 Gradient 1 9862676e04 4614518e01 3883552e04 iteration 12 Parameter 1 7123000e03 1883003e00 3730551e03 Function Value 1 5794977 Gradient 1 1504136e02 9668077e00 2761316e03 Successive iterates within tolerance Current iterate is probably solution theta1lnLmaxlog1estimate Valuefunc1lnLmaxlog1minimum theta1 1 7123000e03 1883003e00 3730551e03 Valuefunc1 1 5794977 Para alfa mais 5 sigmag1 var ystar Tstar 1 Tstar 05 sigmag1 Manter o mesmo sigma 1 242758 ag2max ystar 105 ag2 1 747915 bg2loglogag2logmin ystar bg2 1 1889686 gg20 theta2cag2bg2gg2 funcMLgomp2functiontheta2 Defino meus coeficientes a theta21 a alpha b theta22 b beta g theta23 g gamma sigma var ystar length ystar 1length ystar 05 Criase um vetor de residuos u ystar aexpexpbg tstar lnL 05lengthystarlog2pi05lengthystarlogabssigmag1 05sumu2sigmag12 na frente da funcao para maximizala no algoritmo de otimizacao Agora calculo o valor da funcao return lnL lnLmaxlog2 nlm funcMLgomp2 theta2 steptol 1e6 gradtol 1e8 printlevel 1 iterlim 1000000 ndigit 12 fscale 1 iteration 0 Step 1 0 0 0 Parameter 1 7479150000 1889686 0000000 Function Value 1 7582765 Gradient 1 9393025e04 4648795e01 3912486e04 iteration 23 Parameter 1 7479150e03 1724237e00 3349727e03 Function Value 1 5790175 Gradient 1 1088910e02 1002206e05 3436980e03 Successive iterates within tolerance Current iterate is probably solution tail ystar 1 7123 7123 7123 7123 7123 7123 theta2vlnLmaxlog2estimate Valuefunc2lnLmaxlog2minimum Comparação theta chute ao theta por Gompertz theta2v 1 7479150e03 1724237e00 3349727e03 Valuefunc2 1 5790175 theta2 1 7479150000 1889686 0000000 Para o alfa com 10 ag3max ystar 11 ag3 1 78353 bg3loglogag3logmin ystar bg3 1 1896691 gg30 theta3cag3bg3gg3 funcMLgomp3functiontheta3 Defino meus coeficiente a theta31 a alpha b theta32 b beta g theta33 g gamma sigma var ystar length ystar 1length ystar 05 Criase um vetor de residuos u ystar aexpexpbg tstar Escrevo a fun o de log da verossimilhanca lnL com uma distribuio normal lnL 05lengthystarlog2pi05lengthystarlogabssigmag1 05sumu2sigmag12 na frente da funcao para maximizala no algoritmo de otimizacao Agora calculo o valor da funcao return lnL lnLmaxlog3 nlm funcMLgomp3 theta3 steptol 1e4 gradtol 1e6 printlevel 1 iterlim 100000 ndigit 12 fscale 1 iteration 0 Step 1 0 0 0 Parameter 1 7835300000 1896691 0000000 Function Value 1 7582765 Gradient 1 8966070e04 4681476e01 3940075e04 iteration 8 Parameter 1 7835300e03 1896687e00 3462962e03 Function Value 1 5789136 Gradient 1 2752714e03 1367501e01 3255676e03 Successive iterates within tolerance Current iterate is probably solution tail ystar 1 7123 7123 7123 7123 7123 7123 theta3vlnLmaxlog3estimate Valuefunc3lnLmaxlog3minimum theta3v 1 7835300e03 1896687e00 3462962e03 Valuefunc3 1 5789136 theta3 1 7835300000 1896691 0000000 Para um alfa mais 15 sigmag1 var ystar Tstar 1 Tstar 05 sigmag1 1 242758 ag4max ystar 115 ag4 1 819145 bg4loglogag4logmin ystar bg4 1 190334 gg40 theta4cag4bg4gg4 funcMLgomp4functiontheta4 Defino meus coeficiente a theta41 a alpha b theta42 b beta g theta43 g gamma sigma var ystar length ystar 1length ystar 05 Criase um vetor de residuos u ystar aexpexpbg tstar Escrevo a fun o de log da verossimilhanca lnL com uma distribuio normal lnL 05lengthystarlog2pi05lengthystarlogabssigmag1 05sumu2sigmag12 na frente da funcao para maximizala no algoritmo de otimizacao Agora calculo o valor da funcao return lnL lnLmaxlog4 nlm funcMLgomp4 theta4 steptol 1e4 gradtol 1e6 printlevel 1 iterlim 100000 ndigit 12 fscale 1 iteration 0 Step 1 0 0 0 Parameter 1 819145000 190334 000000 Function Value 1 7582765 Gradient 1 8576240e04 4712704e01 3966439e04 iteration 7 Parameter 1 8191450e03 1903326e00 3321772e03 Function Value 1 5788463 Gradient 1 1678667e03 2876072e01 2546578e03 Successive iterates within tolerance Current iterate is probably solution tail ystar 1 7123 7123 7123 7123 7123 7123 theta4vlnLmaxlog4estimate Valuefunc4lnLmaxlog4minimum theta4v 1 8191450e03 1903326e00 3321772e03 Valuefunc4 1 5788463 Chutes para o Y max Valor da função maximizada Max y 5794977 Max y 105 5790175 Max y 11 5789136 Max y 115 5788463 Como podese notar entre as fun ções maximizadas agravando 5 10 e 15 de alfa o max y115 apresenta uma função menor No entanto para encontrar um valor melhor para estimar essa base foi criado um loop com o objetivo de ver se há um parâmetro maior que para uma verossimilhança melhor ex 35 35Qual dos chutes para α maximizou l θ Você considera que o algoritmo de otimização foi bemsucedido nas estimativas dos parâmetros para cada um dos chutes de α utilizados Justifique devidamente sua resposta Somente considerada se itens 33 34 e 35 estiverem corretos No r acloop max ystar seq125005 bcloop c funclog c for i in 1length acloop bcloop i log acloop i min ystar 1 gcloop 0 chutesloop c acloop i bcloop i gcloop sigma1 cat Valor para i acloop i acloop 1n funclog i nlm funcMLstar chutesloopsteptol 1e4 gradtol 1e6 printlevel 1 iterlim 100000 ndigit 12 fscale 1 minimum E criando a tabela loop2 dataframe Valoresalfa max ystar seq 125005 Valordaagravo seq 125005 Valordafunção funclog 131 loop2 Valoresalfa Valordaagravo Valordafunção 1 712300 100 5843334 2 747915 105 5851183 3 783530 110 5861925 4 819145 115 5875018 5 854760 120 5889885 6 890375 125 5906036 7 925990 130 5922992 8 961605 135 5781840 9 997220 140 5782651 10 1032835 145 5783532 11 1068450 150 57 79 411 12 1104065 155 5785285 13 1139680 160 5786139 14 1175295 165 5786965 15 1210910 170 5787762 16 1246525 175 6063381 17 1282140 180 6075401 18 1317755 185 6086678 19 1353370 190 5790619 20 1388985 195 5791252 21 1424600 200 5791856 22 1460215 205 5792429 23 1495830 210 5792977 24 1531445 215 5793498 25 1567060 220 5793995 26 1602675 225 5794469 27 1638290 230 5794922 28 1673905 235 5795353 29 1709520 240 5795767 30 1745135 245 5796162 31 1780750 250 5796540 Com o melhor chute sigmag1 var ystar Tstar 1 Tstar 05 sigmag1 1 242758 ag5max ystar 17 ag5 1 121091 bg5loglogag5logmin ystar bg5 1 1959972 gg50 theta5cag5bg5gg5 funcMLgomp5functiontheta5 Defino meus coeficiente a theta51 a alpha b theta52 b beta g theta53 g gamma sigma var ystar length ystar 1length ystar 05 Criase um vetor de residuos u ystar aexpexpbg tstar Escrevo a fun o de log da verossimilhanca lnL com uma distribuio normal lnL 05lengthystarlog2pi05lengthystarlogabssigmag1 05sumu2sigmag12 na frente da funcao para maximizala no algoritmo de otimizacao Agora calculo o valor da funcao return lnL lnLmaxlog5 nlm funcMLgomp5 theta5 steptol 1e4 gradtol 1e6 printlevel 1 iterlim 100000 ndigit 12 fscale 1 iteration 0 Step 1 0 0 0 Parameter 1 12109100000 1959972 0000000 Function Value 1 7582765 Gradient 1 5801574e04 4987300e01 4198303e04 iteration 14 Parameter 1 1210893e04 1391858e00 1846769e03 Function Value 1 5779872 Gradient 1 5048631e04 4080972e01 1125287e03 Last global step failed to locate a point lower than x Either x is an approximate local minimum of the function the function is too nonlinear for this algorithm or steptol is too large Irá ser utilizado o valor de 115 de alfa para realizar os chutes a seguir como o 170 de alfa não esta trazendo uma função ideal 36 025 Escreva textualmente a expressão matemática para o modelo estimado utilizando as estimativas para os parâmetros a partir do chute de α que maximizou l θ Somente considerada se itens 33 34 35 e 36 estiverem corretos l θ y 1198 ln 2 π 2 1198 ln 242758 1 2 242758 2 i 1 1198 y i 819145 e e 1903340003321 t 2 l θ y 1198 ln 2 π 2 1198 ln 242758 1 2 242758 2 i 1 1198 y i 819145 e e 1903340003321 t 2 37 025 Para esse modelo estime os errospadrão dos estimadores dos parâmetros sua matriz de variânciacovariância dos estimadores θ e seus respectivos pvalores para os testes de significância individual de cada parâmetro Apresente os códigos e resultados através do console do R agc lnLmaxlog1estimate1 bgc lnLmaxlog1estimate2 ggc lnLmaxlog1estimate3 sigmagc sigmag1 library numDeriv hess hessianfuncMLgomp4c agc bgc ggc sigmagc hess 1 2 3 4 1 6652651e05 1380302e01 1031417e02 4971584e21 2 1380302e01 5331352e02 3207362e05 4169822e17 3 1031417e02 3207362e05 2164100e08 1562207e14 4 4971584e21 4169822e17 1562207e14 0000000e00 lambda 1e7 Valor de regularização ajuste conforme necessário hessreg abs hess lambda 4 Tamanho da diagonal de hess varcovarML solve hessreg varcovarML 1 2 3 4 1 7374297e04 1893200e01 6320482e02 7376184e04 2 1893200e01 2216759e02 4187713e05 1895413e01 3 6320482e02 4187713e05 9680959e08 6324660e02 4 7376184e04 1895413e01 6324660e02 2573781e06 varML abs diag varcovarML varML 1 3729536e05 4569042e02 2378527e07 2873165e06 epML varML051length varML 1 epML 1 6106993e02 2137532e01 4877014e04 z thetah 13epML z 1 12830040 15015565 9694188 pvalor 2 1 pnorm abs z pvalor 1 0 0 0 38 025 Interprete os pvalores calculados no item anterior Somente considerada se itens 33 34 35 36 37 e 38 estiverem corretos Ho pvalor5Parâmetros estatisticamentes iguais a zero H1 pvar5Parâmetros estatisticamentes diferentes de zero pvalor α 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada alfa ser menor que zero é de aproximadamente 0 pvalor β 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada beta ser menor que zero é de aproximadamente 0 pvalor γ 0 5 isto é rejeitamos a hipótese nula pois a probabilidade estimada gama ser menor que zero é de aproximadamente 0 39 05 Utilizando das estimativas para os parâmetros a partir do chute de α que maximizou lθ obtenha as estimativas para ŷt com t 1 T Em seguida faça o plot da curva de crescimento de Gompertz ajustada aos dados Apresente os resultados dos códigos através do console do R Somente considerada se itens 33 34 35 36 e 37 estiverem corretos lnLmaxlog5estimate Yh5 lnLmaxlog5estimate1expexplnLmaxlog5estimate2lnLmaxlog5estimate3tstar Yh5 date2asDatedadosVargemdatawhichdadosVargemdata20200421dadosVargemdata20230801 plot ystardate2 ylab Casos acumulados xlab Data main Casos Acumulados COVID Vargem Grande type p linesYh5date2colblue 310 05 Calcule e apresente através do console do R os resultados para o MSE RMSE MAE e MAPE das estimativas de ŷt geradas no item 27 Somente considerada se itens 23 24 25 26 27 e 210 estiverem corretos residlog ystarYh5 MSE sum residlog 2 Tstar MSE 1 1105966 RMSE sqrtMSE RMSE 1 3325607 MAE sumabs residlog Tstar MAE 1 376 1023 MAPE sumabs residlog ystar Tstar MAPEmean absolute percentage error 1 02452222 311 Calcule os valores previstos de yt para 02082023 até 05092023 Somente considerada se itens 23 24 25 26 27 e 210 estiverem corretos 1 º Criado sequencia de Y previstos a partir dos chute que defini como os que maximizam esta função Yprevisão agc exp exp bgcggc max tstar seq 134 Yprevisão 1 6608148 6609994 6611834 6613668 6615495 6617316 6619131 6620940 6622742 6624538 6626328 12 6628111 6629889 6631660 6633425 6635184 6636937 6638684 6640425 6642160 6643889 6645612 23 6647329 6649039 6650744 6652444 6654137 6655824 6657505 6659181 6660851 6662515 6664173 6665826 2º Defin a data de início e a data de término da base e crio a sequencia datainicio asDate 20230802 datafim asDate 2023090 5 Crie a sequência de datas Datasolicitada seq datainicio datafim by days Formate as datas para o formato desejado dd mm yyyy DataPrevista2 format Datasolicitada format dmY DataPrevista2 1 02082023 03082023 04082023 05082023 06082023 07082023 08082023 8 09082023 10082023 11082023 12082023 13082023 14082023 15082023 15 16082023 17082023 18082023 19082023 20082023 21082023 22082023 22 23082023 24082023 25082023 26082023 27082023 28082023 29082023 29 30082023 31082023 01092023 02092023 03092023 04092023 3º Cri o um data frame combinando os dois vetores tabela2 dataframe Data DataPrevista2 Yprevisão Yprevisão Exiba a tabela printtabela2 Data Yprevisão 1 02082023 6608148 2 03082023 6609994 3 04082023 6611834 4 05082023 6613668 5 06082023 6615495 6 07082023 6617316 7 08082023 6619131 8 09082023 6620940 9 10082023 6622742 10 11082023 6624538 11 12082023 6626328 12 13082023 6628111 13 14082023 6629889 14 15082023 6631660 15 16082023 6633425 16 17082023 6635184 17 18082023 6636937 18 19082023 6638684 19 20082023 6640425 20 21082023 6642160 21 22082023 6643889 22 23082023 6645612 23 24082023 6647329 24 25082023 6649039 25 26082023 6650744 26 27082023 6652444 27 28082023 6654137 28 29082023 6655824 29 30082023 6657505 30 31082023 6659181 31 01092023 6660851 32 02092023 6662515 33 03092023 6664173 34 04092023 6665826 4º Crio um vetor com os valores observados de dados acumulados yobs2 Yt whichdate20230802date2023090 5 yobs2 1 7123 7123 7126 7126 7126 7126 7126 7126 7126 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 21 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7153 7153 7153 7153 yobs2 yobs isna yobs yobs2 1 7123 7123 7126 7126 7126 7126 7126 7126 7126 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 21 7137 7137 7137 7137 7137 7137 7137 7137 7137 7137 7153 7153 7153 7153 312 05 Construa um gráfico comparando os valores observados e previstos de yt no item acima 5º Desidiuse criar uma linha de tendencia para o observado e o previsionado Ajustar um modelo de regressão linear tendencia lm yobs2 Yprevisão tendencia Call lm formula yobs2 Yprevisão Coefficients Intercept Yprevisão 44757717 04008 plotyobs2Yprevisão abline tendencia col red 313 05 Calcule o MSE RMSE MAE e o MAPE das previsões de yt de 02082023 até 05092023 thx1 1length Yprevisão Thx1 length Yprevisão residlog yobs2Yprevisão MSE sum residlog 2Thx1 MSE 1 2484006 RMSE sqrtMSE RMSE 1 498398 MAE sumabs residlog Thx1 MAE 1 4982741 MAPE sum abs residlog Yprevisão Thx1 MAPE 1 007507382 3 1 ponto Qual dos modelos estimados apresentou o melhor ajustamento aos dados Justifique sua resposta Com base no melhor modelo você considera que a curva de casos acumulados por COVID19 está próxima do achatamento total ou ainda há indícios de continuidade deescimento O melhor modelo para o ajustamento da curva de Vargem Grande é o modelo de Gompertz Já que o modelo se ajusta mais a curva da cidade ao longo do período analisado em comparação ao modelo logístico Além disso há menor porcentagem de erro em relação ao estimado pelo modelo em comparação ao ocorrido como conta o MAPE na tabela a seguir MSE RMSE MAE MAPE Modelo logistico 159799 400 334 54 Modelo de Gompertz 1827871 428 376 245 Notase que o Modelo de gompertz apresenta um erro de 243 enquanto o logístico apresenta 413 Comparando os modelos em relação ao observado No r date2asDatedadosVargemdatawhichdadosVargemdata20200421dadosVargemdata20230801 plot ystardate2 ylab Casos acumulados xlab Data main Casos Acumulados COVID Vargem Grande type p lines Yh5date2colbluelinha para o modelo Gompertz lines Yhlogdate2col purple linha para o modelo logístico Modelo de gompertz Azul Logístico Roxo Modelo de gompertz Azul Logístico Roxo Comparando o modelo para previsão plotyobs2tprev ylim c7080 7400 linesYprevisãotprev colblue lines Yhxtprev colpurple No modelo de previsão calculado podese observar que o logístico se encaixa mais no para na estimação da curva No entanto Gompertz teve um melhor ajustamento com um valor de função de 5788463 enquanto o logístico tem 5 78 1 136