Carregando bibliotecas

# "limpar" ambiente
rm(list = ls())

# Carregar bibliotecas
library(readxl)
library(readr) 
library(mfx)
#library(caret)
library(pROC)
#library(ResourceSelection)
#library(modEvA)
library(foreign)
#library(stargazer) #Hlavac, Marek (2022) 
library(ggplot2)
library(glm2)
library(PerformanceAnalytics)
library(nnet) #logit multinomial
library(car)
library(reshape2)
library(lmtest)
library(MASS)
library(pROC)
library(dplyr)
#library(margins)
#library(broom)
library(curl)
#library(faraway)
library(corrplot)
library(psych)
library(factoextra)
library(FactoMineR)
library(stats)

\(m1\) são os dados com todos os tipos de contrato: Concluded, Active, Cancelled e Distressed. Fonte dos dados: https://ppi.worldbank.org/en/visualization. O objetivo é prever a probabilidade de cada status dos contratos.

# Carregar dados m1 contém os contratos ativos 
m1 <- read_excel("C:/Users/kassy/OneDrive - caen.ufc.br/1. CAEN/3. Doutorado/Thesis/Test banco mundial/corrigidaBDprincipal.xlsx")

m1 <- subset(m1, Country != "United Arab Emirates")

# Verificar missing values
colSums(is.na(m1))

##                   Region                  Country              IncomeGroup 
##                        0                        0                        0 
##               IDA Status   Financial closure year  Financial closure Month 
##                        0                        0                        0 
##             Project name             RelatedNames              Type of PPI 
##                        0                     5886                        0 
##                   catppp           Subtype of PPI           Project status 
##                        0                        0                        0 
##             status_dummy           Primary sector                  sectorR 
##                        0                        0                        0 
##                Subsector                  Segment                 Location 
##                        0                      411                      167 
##           ContractPeriod     GovtGrantingContract        DirectGovtSupport 
##                        0                        0                        0 
##   DirectGovtSupportValue      InDirectGovtSupport InDirectGovtSupportValue 
##                     7277                        0                     7509 
##             Total Equity           InvestmentYear                    covid 
##                     5867                        0                        0 
##           PercentPrivate         FeesToGovernment           PhysicalAssets 
##                        0                        0                        0 
##          TotalInvestment         TotalInvestment1                 inv/time 
##                        0                        0                        0 
##             CapacityType                 Capacity               Technology 
##                        0                        0                        0 
##          RelatedProjects              BidCriteria                   catbid 
##                     8084                        0                        0 
##              AwardMethod             NumberOfBids                 Sponsors 
##                        0                     4644                        7 
##         Sponsors Country      Main Revenue Source     Other Revenue Source 
##                        7                        8                        0 
##      MultiLateralSupport         BiLateralSupport         TotalDebtFunding 
##                      723                      784                        0 
##     DebtEquityGrantRatio             ProjectBanks      UnsolicitedProposal 
##                     6070                        0                        0 
##         PublicDisclosure 
##                        0

#colSums(is.na(m2))

Definição das variáveis m1.

Variável dependente \(statusm1\) (Categórica multinomial)

O status dos projetos no banco de dados pode ser uma das seguintes opções:

Ativo (Active): para projetos em construção (ou prestes a iniciar a construção) ou operacionais

Concluidos (Concluded): para os quais o período do contrato expirou e não foi renovado nem prorrogado nem pelo governo nem pelo operador

Projetos cancelados (Cancelled): dos quais o setor privado tenha saído de uma das seguintes formas: - vender ou transferir sua participação econômica de volta ao governo antes de cumprir os termos do contrato

retirar toda a direção e pessoal da empresa
cessar operação, prestação de serviços ou construção por 15% ou mais do período de licença ou concessão, após a revogação da licença ou repúdio do contrato

Projetos em dificuldades (Distressed): em que o governo ou a operadora solicitaram a rescisão do contrato ou estão em arbitragem internacional.

O status do contrato dito como Concluido é definido como baseline para a variável categórica dependente.

Categorização dos tipos de contratos para escrever as equações como:

\(ln(\frac{P(status=Active)}{P(status=Concluded)})\), \(ln(\frac{P(status=Cancelled)}{P(status=Concluded)})\), \(ln(\frac{P(status=Distressed)}{P(status=Concluded)})\)

m1$`Project status` <- as.factor(m1$`Project status`)
m1$`Project status` <- relevel(m1$`Project status`, ref = "Concluded")

statusm1=m1$`Project status`
is.factor(m1$`Project status`)

## [1] TRUE

#with(m1, table(m1$catstatus, m1$`Project status`))

# Crie uma tabela cruzada com as contagens
tabela_contagem <- table(m1$`Project status`)

# Crie um data frame com a contagem e os percentuais
tabela_completa <- data.frame(Contagem = tabela_contagem, Percentual = prop.table(tabela_contagem) * 100)

# Exiba a tabela
print(tabela_completa)

##   Contagem.Var1 Contagem.Freq Percentual.Var1 Percentual.Freq
## 1     Concluded           107       Concluded       1.2721436
## 2        Active          8067          Active      95.9101177
## 3     Cancelled           176       Cancelled       2.0924979
## 4    Distressed            61      Distressed       0.7252408

#write.csv(tabela_completa, "tabela_frequencia.csv", row.names = FALSE)

#with(m1, table(m1$catstatus, m1$`Project status`))

# Crie uma tabela cruzada com as contagens
tabela_contagem <- table(m1$`Type of PPI`)

# Crie um data frame com a contagem e os percentuais
tabela_completa <- data.frame(Contagem = tabela_contagem, Percentual = prop.table(tabela_contagem) * 100)

# Exiba a tabela
print(tabela_completa)

##                   Contagem.Var1 Contagem.Freq               Percentual.Var1
## 1                    Brownfield          1948                    Brownfield
## 2                   Divestiture           658                   Divestiture
## 3            Greenfield project          5403            Greenfield project
## 4 Management and lease contract           402 Management and lease contract
##   Percentual.Freq
## 1       23.160147
## 2        7.823089
## 3       64.237308
## 4        4.779455

#tirar tipos de PPI
#m1 <- subset(m1, `Type of PPI` != "Brownfield")
#m1 <- subset(m1, `Type of PPI` != "Divestiture")

#with(m1, table(m1$catstatus, m1$`Project status`))

# Crie uma tabela cruzada com as contagens
tabela_contagem <- table(m1$`Primary sector`)

# Crie um data frame com a contagem e os percentuais
tabela_completa <- data.frame(Contagem = tabela_contagem, Percentual = prop.table(tabela_contagem) * 100)

# Exiba a tabela
print(tabela_completa)

##                                         Contagem.Var1 Contagem.Freq
## 1                                              Energy          4635
## 2      Information and communication technology (ICT)           237
## 3                               Municipal Solid Waste           415
## 4                                           Transport          1985
## 5 water and sewerage - treatment plants and utilities          1139
##                                       Percentual.Var1 Percentual.Freq
## 1                                              Energy       55.106408
## 2      Information and communication technology (ICT)        2.817739
## 3                               Municipal Solid Waste        4.934015
## 4                                           Transport       23.600048
## 5 water and sewerage - treatment plants and utilities       13.541791

Definição das variáveis controle m1

(\(ITm1\), Investimento Total);

(\(percentprivm1\), Percentual da Participação Privada);

(\(tm1\) Período de Contrato - O período de tempo medido em anos em que os termos de um contrato estão em vigor);

(\(km1\), Investimentos em Ativos Físicos);

(\(debtm1\), Financiamento total da dívida).

ITm1=m1$TotalInvestment1
percentprivm1=m1$PercentPrivate
tm1=m1$ContractPeriod
t3=tm1^3
t2=tm1^2
km1=m1$PhysicalAssets
debtm1=m1$TotalDebtFunding
inv_t=m1$`inv/time`
#catstatusm1=m1$catstatus

sd(tm1)

## [1] 7.834893

sd(inv_t)

## [1] 27.06472

# Alto desvio padrão
hist(tm1)

# TESTE - Inclusão do Número de Bids. Hipótese: Relação positiva entre a quantidade de bids e o insucesso dos projetos.

nbids=m1$NumberOfBids

Variáveis categóricas m1

(\(catpppm1\), Tipo de PPI) Os tipos de contratos para a definição das variáveis categóricas não ordenadas:

Brownfield (Brownfield são semelhantes aos Greenfields, exceto que, em vez de construir um novo ativo, a entidade privada assume um ativo existente e geralmente faz uma melhoria nele (reabilitação) ou o expande. Eles geralmente assumem as operações do ativo existente primeiro e, em seguida, realizam o investimento de capital. Tal como a Greenfields, a entidade privada normalmente tem responsabilidade operacional por um determinado período de tempo, durante o qual recupera o seu investimento da operação do projeto, após o que o projeto pode reverter para o setor público.)
Divestiture (Uma entidade privada compra uma participação acionária em uma empresa estatal por meio de uma venda de ativos, oferta pública ou programa de privatização em massa. Uma alienação, como uma concessão, dá ao operador privado total responsabilidade pela operação, manuntenção e investimento.)
Greenfield project (Uma entidade privada ou uma joint venture público-privada constrói e explora uma nova instalação durante o período especificado no contrato de projeto. A entidade privada assume grande parte do risco financeiro e operacional e recupera seus investimentos ao longo da vida do projeto.)
Management and Lease Contracts (Uma entidade privada assume a gestão de um bem público por um período fixo, enquanto as decisões de propriedade e investimento permanecem com o Estado.)

O tipo Brownfield é usado como baseline para a categoria.

(\(catbidm1\)) Tipos de Leilão para os contratos e a definição das variáveis categóricas não ordenadas:

(vazio) Not Applicable Not Available Other
Highest new investment
Highest percentage of revenue share with government
Highest price paid to government
Lowest cost of construction or operation
Lowest government payments
Lowest subsidy required
Lowest tariff

A informação “1.(vazio) Not Applicable Not Available Other” é usada como baseline da categoria.

catpppm1=m1$catppp
catbidm1=m1$catbid
catpppm1 = factor(catpppm1)
catbidm1 = factor(catbidm1)
inv.year = m1$InvestmentYear
inv.year = factor(inv.year)
#catsubtypem1=m1$catsubtype
#catsetorm1=m1$catsetor
#catsetor2m1=m1$catsetor2
#catsetorm1=m1$catsetor3

Definição setorial categórica

A Base de Dados de Participação Privada em Projetos de Infraestrutura são divididas em setores com seus correspodentes subsetores, para a definição das variáveis categóricas não ordenadas:

energy - electricity and natural gas
transport - airports, ports, railways, and roads
water and sewerage - treatment plants and utilities
ICT – ICT backbone
municipal solid waste -collection and transport, treatment/disposal, integrated municipal solid waste

O setor energy é definido como baseline para a variável categórica setorial.

m1$sectorR <- as.factor(m1$sectorR)
m1$sectorR <- relevel(m1$sectorR, ref = "Enrg")
sector = m1$sectorR

is.factor(m1$sectorR)

## [1] TRUE

Definição de categórica para o metódo de adjudicar usado pelo governo no contrato

Competitive bidding
Competitive negotiation
Direct negotiation
License Scheme

A método usado como baseline é a categoria Not Available.

# FACTOR
award = m1$AwardMethod
award = factor(award)
award <- relevel(award, ref = "Not Available")

Categoria de Renda

income = m1$IncomeGroup
income = factor(income)
income <- relevel(income, ref = "Lower middle income")

#catpppm1 = factor(catpppm1)
#catbidm1 = factor(catbidm1)
#catsetorm1 = factor(catsetorm1)
#catsetor2m1 = factor(catsetor2m1)
#sector = factor(sector)
#catsetorm1 = factor(catsetorm1)
#catstatusm1 = factor(catstatusm1)
#catsubtypem1 = factor(catsubtypem1) 
covid = m1$covid
covid = factor(covid)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ## Criação de fatores%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#Data frame com as variáveis númericas
data1 <- data.frame(inv_t, tm1, debtm1,t2 )
datacat <- data.frame(catbidm1, catpppm1, sector)

matcor <- cor(data1)

print(matcor, digits = 2)

##         inv_t    tm1 debtm1     t2
## inv_t   1.000 -0.091  0.557 -0.068
## tm1    -0.091  1.000  0.048  0.927
## debtm1  0.557  0.048  1.000  0.043
## t2     -0.068  0.927  0.043  1.000

corrplot(matcor, method = "circle")

cortest.bartlett(data1)

## $chisq
## [1] 19817.87
## 
## $p.value
## [1] 0
## 
## $df
## [1] 6

#Ho: A matriz de correlação da população é uma matriz identidade, ou seja as variáveis não são correlacionadas na população.

#H1: A matriz de correlação da população não é uma matriz identidade, ou seja as variáveis são correlacionadas na população.

KMO(data1)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data1)
## Overall MSA =  0.5
## MSA for each item = 
##  inv_t    tm1 debtm1     t2 
##   0.49   0.50   0.49   0.50

A estatística KMO maior que 0,5 também concorda quanto ao fato de que a análise fatorial pode ser considerada uma técnica apropriada para analisar a matriz de correlação.

Análise Fatorial

#cor = TRUE: as componentes principais serão geradas a partir da matriz de correlação.
#cor = FALSE: as componentes principais serão geradas a partir da matriz de covariância

fit <- princomp(data1, cor=TRUE)
fit

## Call:
## princomp(x = data1, cor = TRUE)
## 
## Standard deviations:
##    Comp.1    Comp.2    Comp.3    Comp.4 
## 1.3929774 1.2465199 0.6582690 0.2692289 
## 
##  4  variables and  8411 observations.

summary(fit)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4
## Standard deviation     1.3929774 1.2465199 0.6582690 0.26922891
## Proportion of Variance 0.4850965 0.3884529 0.1083295 0.01812105
## Cumulative Proportion  0.4850965 0.8735494 0.9818789 1.00000000

A função \(summary(fit)\) mostra a aplicação da análise de componentes principais. O fator 1 responde por 81,05% da variância total. Da mesma forma, o segundo fator responde por 12,62% da variância total, sendo que os dois primeiros fatores respondem por 93,67% da variância total. Várias considerações devem integrar a análise do número de fatores que devem ser usados na análise.

Determinação do número de fatores

Autovalores Como o autovalor representa a quantidade de variância associada ao fator, incluem-se apenas os fatores com variância maior que 1.

Percentagem da variância Determina que o núumero de fatores extraídos seja de no mínimo 60% da variância

screeplot(fit)

plot(fit,type="lines")

ACP

PCAdente <- principal(data1, nfactors=2, n.obs=344,rotate="none", scores=TRUE)

PCAdente

## Principal Components Analysis
## Call: principal(r = data1, nfactors = 2, rotate = "none", n.obs = 344, 
##     scores = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##          PC1  PC2   h2    u2 com
## inv_t  -0.17 0.87 0.78 0.216 1.1
## tm1     0.98 0.07 0.96 0.037 1.0
## debtm1 -0.01 0.89 0.79 0.214 1.0
## t2      0.98 0.08 0.96 0.039 1.0
## 
##                        PC1  PC2
## SS loadings           1.94 1.55
## Proportion Var        0.49 0.39
## Cumulative Var        0.49 0.87
## Proportion Explained  0.56 0.44
## Cumulative Proportion 0.56 1.00
## 
## Mean item complexity =  1
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.09 
##  with the empirical chi square  822.47  with prob <  NA 
## 
## Fit based upon off diagonal values = 0.96

PCAdentevarimax <- principal(data1, nfactors=3,
            n.obs=344, rotate="varimax",scores=TRUE)

PCAdentevarimax

## Principal Components Analysis
## Call: principal(r = data1, nfactors = 3, rotate = "varimax", n.obs = 344, 
##     scores = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##          RC1  RC2   RC3   h2      u2 com
## inv_t  -0.05 0.29  0.95 1.00 7.6e-05 1.2
## tm1     0.98 0.03 -0.05 0.96 3.6e-02 1.0
## debtm1  0.03 0.96  0.29 1.00 4.1e-05 1.2
## t2      0.98 0.01 -0.02 0.96 3.6e-02 1.0
## 
##                        RC1  RC2  RC3
## SS loadings           1.93 1.00 1.00
## Proportion Var        0.48 0.25 0.25
## Cumulative Var        0.48 0.73 0.98
## Proportion Explained  0.49 0.25 0.25
## Cumulative Proportion 0.49 0.75 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 3 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.01 
##  with the empirical chi square  22.17  with prob <  NA 
## 
## Fit based upon off diagonal values = 1

PCAdentevarimax$values

## [1] 1.94038594 1.55381177 0.43331808 0.07248421

PCAdentevarimax$loadings

## 
## Loadings:
##        RC1    RC2    RC3   
## inv_t          0.293  0.955
## tm1     0.980              
## debtm1         0.956  0.292
## t2      0.982              
## 
##                  RC1   RC2   RC3
## SS loadings    1.927 1.001 0.999
## Proportion Var 0.482 0.250 0.250
## Cumulative Var 0.482 0.732 0.982

biplot(PCAdentevarimax)

pca=PCA(data1, graph = TRUE)

fviz_eig(pca, addlabels=TRUE, ylim = c(0,100))

variaveis=get_pca_var(pca)
head(variaveis$coord)

##               Dim.1      Dim.2       Dim.3        Dim.4
## inv_t  -0.169065045 0.86923648  0.46450890  0.008742880
## tm1     0.978742931 0.07258061  0.02021152  0.190750690
## debtm1 -0.005832947 0.88661022 -0.46243598 -0.006422999
## t2      0.976642824 0.08303066  0.05739352 -0.189685762

pca <- princomp(data1)

summary(m1$ContractPeriod)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   25.00   25.33   25.33   30.00   99.00

summary(statusm1)

##  Concluded     Active  Cancelled Distressed 
##        107       8067        176         61

summary(percentprivm1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00  100.00  100.00   92.04  100.00  100.00

summary(inv_t)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    0.001    1.197    3.855   10.247    9.589 1423.460

summary(catpppm1)

##    1    2    3    4 
## 1948  658 5403  402

summary(sector)

##  Enrg   ICT   MSW   Trp Water 
##  4635   237   415  1985  1139

sd(m1$ContractPeriod)

## [1] 7.834893

sd(percentprivm1)

## [1] 17.92253

sd(inv_t)

## [1] 27.06472

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ESTIMAÇÃO DO LOGIT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Logit Multinomial

\(ln(\frac{P(status=Active)}{P(status=Concluded)}) = b_0 + ITm1 +km1 + tm1 + debtm1 + percentprivm1 + catpppm1 + sector + catbidm1 + award\);

\(ln(\frac{P(status=Cancelled)}{P(status=Concluded)}) = b_0 + ITm1 +km1 + tm1 + debtm1 + percentprivm1 + catpppm1 + sector + catbidm1 + award\);

\(ln(\frac{P(status=Distressed)}{P(status=Concluded)}) = b_0 + ITm1 +km1 + tm1 + debtm1 + percentprivm1 + catpppm1 + sector + catbidm1 + award\).

Estimação do modelo

logm1 = multinom(statusm1 ~ tm1 + t2 + inv_t + sector + catpppm1, data = m1)

## # weights:  48 (33 variable)
## initial  value 11660.121871 
## iter  10 value 2146.770160
## iter  20 value 1604.590937
## iter  30 value 1501.413700
## iter  40 value 1393.941749
## iter  50 value 1348.406350
## iter  60 value 1345.867315
## final  value 1345.861111 
## converged

#+ catpppm1 tirei no dia 10/06/2025

#t3 que deu certo: tm1 + t3 + ITm1 + sector + percentprivm1

#logm1 = multinom(statusm1 ~ tm1 + pca$scores[,1] + catpppm1 + tm1 * catpppm1 + catpppm1 +  tm1 * sector +  sector + tm1 * catbidm1 + catbidm1 + award, data=m1)
#ITm1 + km1 + tm1 + debtm1 + percentprivm1 

#+ ITm1 + km1 + debtm1 + income + percentprivm1 + catbidm1 + income + catpppm1
summary(logm1)

## Call:
## multinom(formula = statusm1 ~ tm1 + t2 + inv_t + sector + catpppm1, 
##     data = m1)
## 
## Coefficients:
##            (Intercept)       tm1           t2      inv_t  sectorICT sectorMSW
## Active       -2.964156 0.4708203 -0.004163465 0.09248866 -0.6805172  43.52247
## Cancelled    -5.475753 0.4096986 -0.003323614 0.09193002  0.6283779  41.77611
## Distressed   -6.146477 0.4559160 -0.003892040 0.09403359  0.5698552 -23.68227
##             sectorTrp sectorWater catpppm12  catpppm13  catpppm14
## Active     -0.8768102  -0.1029293  27.14981  2.1981761 -0.9136416
## Cancelled  -0.6553114   0.5424767  27.48827  0.9483161 -1.1047812
## Distressed -1.7616301  -0.4341678  25.87199 -0.2142774 -2.8826977
## 
## Std. Errors:
##            (Intercept)        tm1           t2      inv_t sectorICT
## Active      0.23845731 0.02599960 0.0006422415 0.01231957 0.1703574
## Cancelled   0.19962872 0.02588962 0.0006420580 0.01254584 0.1555656
## Distressed  0.07759853 0.02322122 0.0007204774 0.01277797 0.1355723
##               sectorMSW sectorTrp sectorWater catpppm12 catpppm13  catpppm14
## Active     1.637249e-02 0.1644795   0.1627522 0.1588936 0.1334487 0.21690741
## Cancelled  1.637249e-02 0.1962046   0.1848374 0.1760423 0.1891366 0.25637849
## Distressed 5.965721e-28 0.1860192   0.2953300 0.2258742 0.1753464 0.01991714
## 
## Residual Deviance: 2691.722 
## AIC: 2757.722

#install.packages("mlogit")
#library(mlogit)
#install.packages("dplyr")
#library(dplyr)
#install.packages("tidyverse")
#library(tidyverse)

##
#m1 <- m1 %>% mutate(id = row_number())

# 4. Garanta que a variável de status é fator
#m1$`Project status` <- as.factor(m1$`Project status`)
#m1$`Project status` <- relevel(m1$`Project status`, ref = "Concluded")  # ajuste conforme seu interesse

# 5. Pegue os níveis (alternativas possíveis)

#long_df <- long_df %>%
 # mutate(
  #  t2 = tm1^2,
    # outras transformações aqui, se necessário
  #)


#alternativas <- levels(m1$`Project status`)

# 6. Expanda a base: uma linha por projeto por alternativa
#long_df <- m1 %>%
 # expand_grid(statusm1 = alternativas) %>%
  #mutate(chosen = ifelse(statusm1 == `Project status`, 1, 0))

# 7. Prepare o data.frame no formato dfidx (para o mlogit)
#mlogit_df <- dfidx(long_df, idx = c("id", "statusm1"), choice = "chosen")

# 8. Estime o modelo logit multinomial
#modelo <- mlogit(chosen ~ tm1 + t2 + inv_t + sector + percentprivm1 + catpppm1, 
 #                data = mlogit_df)

# 9. Veja o resumo do modelo
#summary(modelo)

Testar heterocedasticidade

A proposta seria testar heterocedasticidade no modelo afim de verificar se a variância do período de contrato é a mesma ao longo de toda a amostra. Esperamos que o resultado seja heterocedasticidade para \(tm1\), que poderiamos interpretrar como o ambiente em que são firmados os contratos são não-ergodicos, ou seja, possuem alta variabilidade no tempo, que pode ser concluido como um ambiente incerto sobre os contratos, assim demonstrando a incompletude contratual.

# Obtenha as probabilidades previstas para cada categoria
#probs <- fitted(logm1)

# Construa a matriz de respostas (one-hot encoding)
# Suponha que a variável resposta seja 'Tipo'
#resposta <- model.response(model.frame(logm1))

# Converta para matriz dummy (um 1 para a categoria observada, 0 para as outras)
#Y <- model.matrix(~ resposta - 1)

# Calcule resíduos de Pearson manualmente
#residuos_pearson <- (Y - probs) / sqrt(probs * (1 - probs))

# Para visualizar, pegue uma categoria específica (ex: a primeira)
#residuos_cat1 <- residuos_pearson[,1]

# Agora o plot funciona — ex: contra a variável x1
#plot(tm1, residuos_cat1, main = "Resíduos (categoria 1) vs t", xlab = "t", ylab = "Resíduos")
#abline(h = 0, col = "red")

#residuos <- residuals(logm1, type = "pearson")
#modelo_residual = lm(residuos^2 ~ tm1, data=m1)
#white_test <- bptest(modelo_residual, studentize = FALSE)
#print(white_test)

Propósito do Modelo Residual O objetivo principal do modelo residual é verificar a presença de heterocedasticidade. Especificamente, ele testa se a variância dos resíduos do modelo original está relacionada a uma ou mais das variáveis explicativas. Se essa relação for significativa, isso sugere que a variância dos resíduos não é constante, indicando heterocedasticidade.

Interpretação Se os coeficientes do modelo residual são significativos: Isso indica que a variância dos resíduos do modelo original está relacionada às variáveis explicativas, sugerindo a presença de heterocedasticidade. Se os coeficientes do modelo residual não são significativos: Isso sugere que a variância dos resíduos é constante, indicando homocedasticidade.

O valor-p muito baixo (< 2.2e-16) indica forte evidência contra a hipótese nula de homocedasticidade, sugerindo que há heterocedasticidade nos resíduos do modelo.

Método Stepwise

O método auxilia em selecionar as variáveis importantes para o modelo.

#step1 <- step(logm1, direction = 'both')

O modelo se ajusta melhor sem \(k_{m1}\)

Modelo estimado novamente sem a variável \(k_{m1}\).

#logm1 = multinom(statusm1 ~ ITm1 + tm1 + catpppm1 + sector + catbidm1 + award, data=m1)

#summary(logm1)

#teste_data = data.frame(ITm1, tm1, debtm1, percentprivm1)
#, catpppm1, sector, catbidm1, award
# Instale e carregue o pacote 'usdm'
#install.packages("usdm")
#library(usdm)

# Crie um conjunto de dados apenas com variáveis numéricas
#numeric_data <- m1[, c("ITm1", "tm1", "debtm1", "percentprivm1")]


# Ajuste o modelo logm1 novamente usando apenas as variáveis numéricas
#logm1_numeric <- multinom(statusm1 ~ ITm1 + tm1 + debtm1 + percentprivm1, data = teste_data)

# Calcule os VIFs
#vif_results <- vif(logm1_numeric)
#print(vif_results)

#vif(logm1)

#coef_logm1 <- tidy(logm1)
#print(coef_logm1, n = 69)

Análise da razão estatística de máxima verossimilhança (LR statistic)

lr_test <- lrtest(logm1)

## # weights:  8 (3 variable)
## initial  value 11660.121871 
## iter  10 value 2141.624399
## iter  20 value 1784.936020
## iter  20 value 1784.936020
## iter  20 value 1784.936020
## final  value 1784.936020 
## converged

print(lr_test)

## Likelihood ratio test
## 
## Model 1: statusm1 ~ tm1 + t2 + inv_t + sector + catpppm1
## Model 2: statusm1 ~ 1
##   #Df  LogLik  Df  Chisq Pr(>Chisq)    
## 1  33 -1345.9                          
## 2   3 -1784.9 -30 878.15  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#coef_test(logm1)

O teste de razão de verossimilhança compara a verossimilhança dos dois modelos para determinar se o Model 1, que inclui variáveis independentes, é significativamente melhor em ajustar os dados do que o Model 2, que é um modelo muito mais simples sem nenhuma variável independente.

O valor de Chisq é a diferença nas log-verossimilhanças entre os dois modelos. O valor p (Pr(>Chisq)) é usado para testar a hipótese nula de que o Model 2 é tão bom quanto o Model 1 em ajustar os dados. No caso:

O valor de Chisq é 1040.7. O valor p associado a este teste é menor do que 2.2e-16, o que indica que há uma diferença significativa entre os dois modelos. Portanto, com base no resultado do teste de razão de verossimilhança, podemos concluir que o Model 1, que inclui todas as variáveis independentes, é significativamente melhor em ajustar os dados do que o Model 2, que não inclui nenhuma variável independente além do intercepto. Isso sugere que pelo menos uma das variáveis independentes em Model 1 é significativa para explicar a variável dependente statusm1.

*Pseudo \(R^2\)

Semelhante ao coeficiente de determinação \(R^2\) da regressão múltipla, a medida de pseudo \(R^2\) representam o ajuste geral do modelo proposto. Sua interpretação, portanto, é semelhante à regressão múltipla. Abaixo segue o cálculo do pseudo \(R^2\).

\(R{^2}_{LOGIT} = \frac{-2LL_{nulo}-(-2LL_{modelo})}{-2LL_{nulo}}\)

Lembrando que o valor -2LL representa -2 vezes o logaritmo do valor de verossimilhança, onde a verossimilhança do modelo nulo é comparado com o modelo completo. Abaixo mostra-se o código para calcular os indicadores, sendo que constam as medidas de pseudo \(R^2\) estipuladas por Cox e Snell, Nagelkerke e McFadden.

# Calcular a log-verossimilhança do modelo completo (logm1)
log_likelihood_model <- logLik(logm1)

# Calcular a log-verossimilhança do modelo nulo (intercept apenas)
# Você pode criar um modelo nulo usando a função multinom com apenas ~ 1
null_model <- multinom(statusm1 ~ 1, data = m1)

## # weights:  8 (3 variable)
## initial  value 11660.121871 
## iter  10 value 2141.624399
## iter  20 value 1784.936020
## iter  20 value 1784.936020
## iter  20 value 1784.936020
## final  value 1784.936020 
## converged

log_likelihood_null <- logLik(null_model)

# Calcular o R² de McFadden
R2_McFadden <- 1 - (log_likelihood_model / log_likelihood_null)

# Exibir o resultado
print(R2_McFadden)

## 'log Lik.' 0.2459892 (df=33)

Estatítica \(Z\)

z = summary(logm1)$coefficients/summary(logm1)$standard.errors

z

##            (Intercept)      tm1        t2    inv_t sectorICT     sectorMSW
## Active       -12.43055 18.10875 -6.482710 7.507460 -3.994645  2.658268e+03
## Cancelled    -27.42969 15.82482 -5.176501 7.327529  4.039311  2.551604e+03
## Distressed   -79.20867 19.63359 -5.402030 7.359038  4.203331 -3.969725e+28
##            sectorTrp sectorWater catpppm12 catpppm13   catpppm14
## Active     -5.330817  -0.6324298  170.8679 16.472068   -4.212127
## Cancelled  -3.339938   2.9348857  156.1459  5.013920   -4.309181
## Distressed -9.470150  -1.4701108  114.5416 -1.222024 -144.734514

\(p\)-valor

Nível de significância.

# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2 
 
p

##            (Intercept) tm1           t2        inv_t    sectorICT sectorMSW
## Active               0   0 9.008971e-11 6.039613e-14 6.479136e-05         0
## Cancelled            0   0 2.260860e-07 2.344791e-13 5.360842e-05         0
## Distressed           0   0 6.589109e-08 1.851852e-13 2.630161e-05         0
##               sectorTrp sectorWater catpppm12    catpppm13    catpppm14
## Active     9.777197e-08 0.527106047         0 0.000000e+00 2.529770e-05
## Cancelled  8.379696e-04 0.003336706         0 5.333208e-07 1.638605e-05
## Distressed 0.000000e+00 0.141531754         0 2.216987e-01 0.000000e+00

#logitor(formula = statusm1 ~ ITm1 + tm1 + debtm1 + percentprivm1 + catpppm1 + sector + catbidm1 + award, data = m1)

Razão de chances

O modelo de regressão logística, porém, traz os resultados dos estimadores na forma logarítma, ou seja, o log das chances.

##Razão de chances?
## extract the coefficients from the model and exponentiate
exp(coef(logm1))

##            (Intercept)      tm1        t2    inv_t sectorICT    sectorMSW
## Active     0.051603979 1.601307 0.9958452 1.096901  0.506355 7.972057e+18
## Cancelled  0.004187073 1.506364 0.9966819 1.096288  1.874567 1.390379e+18
## Distressed 0.002141012 1.577618 0.9961155 1.098597  1.768011 5.187037e-11
##            sectorTrp sectorWater    catpppm12 catpppm13  catpppm14
## Active     0.4161081   0.9021907 618035901206 9.0085676 0.40106106
## Cancelled  0.5192803   1.7202621 866971199066 2.5813591 0.33128334
## Distressed 0.1717646   0.6478036 172212237366 0.8071245 0.05598353

Algumas conclusões específicas:

Chances Relativas (Odds Ratios) para a categoria “Cancelled”:

\(ITm1\): Um aumento de uma unidade em “\(ITm1\)” está associado a um aumento de aproximadamente 1.5% nas chances de um contrato ser cancelado.

\(tm1\): Aumentar \(tm1\), ou seja, o tempo de contrato em uma unidade está relacionado a um aumento de cerca de 31% nas chances de um contrato ser cancelado.

\(debtm1\): Um aumento de uma unidade em “\(debtm1\)” está associado a uma diminuição de aproximadamente 1.6% nas chances de um contrato ser cancelado.

\(percentprivm1\): Um aumento de uma unidade em “\(percentprivm1\)” tem um efeito mínimo nas chances de um evento ser cancelado (aproximadamente 0.2% de aumento), além disso o resultado não possui significância estatística.

Categorias e Setores:

\(catpppm12\), \(catpppm13\): Cada uma dessas categorias em comparação com a referência tem uma influência significativa nas chances de um evento ser cancelado. Contudo a categoria \(catpppm14\) (Management and Lease Contracts: Uma entidade privada assume a gestão de um bem público por um período fixo, enquanto as decisões de propriedade e investimento permanecem com o Estado.) possui mais influência para o evento ser concluido do que cancelado.

Setores (sector): A categoria “\(sectorMSW\)” tem um impacto particularmente forte, aumentando drasticamente as chances de um evento ser cancelado em comparação com a categoria de referência. Os outros setores \(ICT\) e \(Water\) possuem também possuem alta chance de serem cancelados, 4.61 vezes as chances de ocorrer no setor ICT e 3.06 vezes as chances de serem cancelados no setor \(Water\). O único setor que há uma chance maior de ser concluido do que cancelado é o setor de Transportes “\(sectorTRP\)”, com 0.644 de chances de cancelamento, contudo esse resultado só possui significância estatística ao nível de 10%.

Categorias de Bid e Award:

\(catbidm12\), \(catbidm13\), \(catbidm14\), \(catbidm15\), \(catbidm16\), \(catbidm17\), \(catbidm18\): Cada uma dessas categorias em comparação com a referência tem uma influência significativa nas chances de um evento ser cancelado, porém \(catbidm15\) não possui significância estatística.

Award Categories (Competitive bidding, Competitive negotiation, Direct negotiation, License scheme): Cada uma dessas categorias em comparação com a referência tem um impacto nas chances de um evento ser cancelado. Especificamente, “Competitive negotiation” e “License scheme” têm efeitos notáveis. A categoria Competitive negotiation não possui significância estatítica.

# Dados para a categoria Cancelled
cancelled_results <- data.frame(
  Variable = colnames(exp(coef(logm1))),
  Odds_Ratio = exp(coef(logm1)["Cancelled", ]),
  P_Value = p["Cancelled", ]
)

# Exibir a tabela
cancelled_results #n=23

Efeito no status “Distressed”:

Comparado com a categoria “Cancelled”, a categoria “Distressed” tem uma chance 0.43 vezes menor de ser cancelada, indicando uma menor propensão ao cancelamento.

Para o status “Distressed” boa parte dos resultados são semelhantes aos do status “Cancelled”, contudo há algumas diferenças:

\(catppm13\) (Greenfield project) possui mais chances de ser concluido do que chances de um status distressed, a uma ordem de 0.64;

\(sectorMSW\) que possuia uma das maiores chances entre os setores de ser cancelado, quando considerado (\(distressed/concluded\)) possui maiores chances de ser concluido do que distressed, a ordem de 1.645388e-05.

As categorias \(catbidm12\), \(catbidm13\), \(catbidm15\), \(catbidm16\) e \(catbidm17\) também possui alta chance de serem concluidos do que distressed.

\(awardCompetitive negotiation\) foi o único da categoria que modificou a razão de chances e passou a ter mais chances de ser concluido, quando comparado a dinamica do status cancelled.

# Dados para a categoria Distressed
distressed_results <- data.frame(
  Variable = colnames(exp(coef(logm1))),
  Odds_Ratio = exp(coef(logm1)["Distressed", ]),
  P_Value = p["Distressed", ]
)

# Exibir a tabela
distressed_results

# Dados para a categoria Distressed
#concluded_results <- data.frame(
 # Variable = colnames(exp(coef(logm1))),
  #Odds_Ratio = exp(coef(logm1)["Concluded", ]),
  #P_Value = p["Concluded", ]
#)

# Exibir a tabela
#concluded_results

Matriz de Confusão

# Instale e carregue a biblioteca caret se ainda não o fez
# install.packages("caret")
# library(caret)

# Obtenha as previsões do modelo
predictions <- predict(logm1, newdata = m1, type = "probs")

# Obtenha as classes previstas
predicted_classes <- max.col(predictions, "first")

# Crie a matriz de confusão
conf_matrix <- table(predicted_classes, statusm1)

# Exiba a matriz de confusão
conf_matrix

##                  statusm1
## predicted_classes Concluded Active Cancelled Distressed
##                 1        56     15         0          0
##                 2        51   8052       176         61

Concluded: 50 observações foram previstas corretamente como “Concluded” (linhas), enquanto 57 observações foram erroneamente previstas como “Concluded” quando na verdade eram “Active” (coluna “Active”).

Active: 8048 observações foram previstas corretamente como “Active”, 20 foram erroneamente previstas como “Active” quando na verdade eram “Concluded”, 176 foram erroneamente previstas como “Active” quando na verdade eram “Cancelled”, e 61 foram erroneamente previstas como “Active” quando na verdade eram “Distressed”.

Cancelled e Distressed: O padrão se repete para as outras classes.

A diagonal principal da matriz representa as previsões corretas (verdadeiros positivos), enquanto as entradas fora da diagonal principal representam os erros de previsão. Isso fornece uma visão geral do desempenho do modelo para cada classe.

Linha 1 (predicted_classes = 1): Essa linha representa as previsões feitas pelo modelo em que ele classificou as observações como classe 1. Em termos do seu problema, a classe 1 parece ser “Concluded”. Portanto, a entrada “50” na coluna “Concluded” significa que o modelo previu corretamente 50 observações como “Concluded”. As outras entradas (20 em “Active”, 0 em “Cancelled” e 0 em “Distressed”) representam previsões incorretas para essas classes.

Linha 2 (predicted_classes = 2): Esta linha representa as previsões feitas pelo modelo em que ele classificou as observações como classe 2. Parece que, neste caso, a classe 2 é “Active”. Portanto, a entrada “8048” na coluna “Active” significa que o modelo previu corretamente 8048 observações como “Active”. As outras entradas representam previsões incorretas para as outras classes.

Cada linha na matriz de confusão representa as previsões feitas pelo modelo para uma classe específica, e as colunas representam as classes verdadeiras. As entradas na diagonal principal representam previsões corretas, e as entradas fora da diagonal principal representam erros de previsão.

A coluna “Taxa de Precisão (%)” mostra a porcentagem de observações corretamente classificadas em relação ao total de observações reais para cada categoria e a taxa global. Portanto, a tabela fornece uma visão detalhada de como o modelo de logit multinomial está se saindo na classificação das diferentes categorias em comparação com as classificações reais. As taxas de precisão indicam a eficácia do modelo em cada categoria e globalmente.

A taxa de precisão, ou taxa de acerto, é uma métrica que indica a proporção de observações corretamente classificadas pelo modelo em relação ao total de observações. A fórmula básica para calcular a taxa de precisão é:

\[ \textrm{Taxa de precisão}(\%) = \frac{\textrm{Número de Observações Corretamente Classificadas}}{\textrm{Total de Observações}} \]

# Obter as previsões do modelo
predictions <- predict(logm1, newdata = m1, type = "class")

# Criar uma matriz de confusão
conf_matrix <- table(predictions, statusm1)

# Calcular as taxas de precisão para cada classe
taxa_precisao <- diag(conf_matrix) / rowSums(conf_matrix) * 100

# Calcular a taxa de precisão global
taxa_precisao_global <- sum(diag(conf_matrix)) / sum(conf_matrix) * 100
#taxa_precisao_global_arredondada <- round(taxa_precisao_global, digits = 15)
# Criar uma matriz para armazenar as taxas de precisão
taxas_matrix <- matrix(NA, nrow = 2, ncol = ncol(conf_matrix))
taxas_matrix[1, 1:ncol(conf_matrix)] <- taxa_precisao
taxas_matrix[2, 1:ncol(conf_matrix)] <- c(NA, taxa_precisao_global)

# Remover entradas "NaN"
taxas_matrix[is.na(taxas_matrix)] <- 0

# Adicionar os nomes das linhas e colunas
rownames(taxas_matrix) <- c("Taxa de Precisao", "Taxa de Precisao Global")
colnames(taxas_matrix) <- colnames(conf_matrix)

# Imprimir a matriz de confusão com as taxas de precisão
print(taxas_matrix)

##                         Concluded   Active Cancelled Distressed
## Taxa de Precisao         78.87324 96.54676         0    0.00000
## Taxa de Precisao Global   0.00000 96.39757         0   96.39757

Para a classe “Active”:

\(\textrm{Taxa de Precisão} = \frac{8048}{57+8048+176+61}=96.48\%\)

Para a classe “Concluded”:

\(\textrm{Taxa de precisão}=\frac{50}{50+20+0+0)}=71.43\%\)

Para as outras classes (“Cancelled” e “Distressed”), não há previsões corretas, então a taxa de precisão é zero.

A taxa de precisão global é a média dessas taxas de precisão.

# Obter as previsões do modelo
predictions <- predict(logm1, newdata = m1, type = "class")

# Criar uma matriz de confusão
conf_matrix <- table(predictions, statusm1)

# Calcular as taxas de precisão para cada classe
taxa_precisao <- diag(conf_matrix) / rowSums(conf_matrix) * 100

# Calcular a taxa de precisão global
taxa_precisao_global <- sum(diag(conf_matrix)) / sum(conf_matrix) * 100

# Adicionar as taxas de precisão como uma linha na matriz de confusão
conf_matrix_com_taxas <- rbind(conf_matrix, "Taxa de Precisao" = taxa_precisao, "Taxa de Precisao Global" = taxa_precisao_global)

# Adicionar os nomes das linhas
rownames(conf_matrix_com_taxas) <- c(levels(factor(statusm1)), "Taxa de Precisao", "Taxa de Precisao Global")

# Imprimir a matriz de confusão com as taxas de precisão
print(conf_matrix_com_taxas)

##                         Concluded     Active Cancelled Distressed
## Concluded                56.00000   15.00000   0.00000    0.00000
## Active                   51.00000 8052.00000 176.00000   61.00000
## Cancelled                 0.00000    0.00000   0.00000    0.00000
## Distressed                0.00000    0.00000   0.00000    0.00000
## Taxa de Precisao         78.87324   96.54676       NaN        NaN
## Taxa de Precisao Global  96.39757   96.39757  96.39757   96.39757

classes_previstas <- predict(logm1, m1, type = "class")

# Criar a matriz de confusão
conf_matrix <- table(classes_previstas, statusm1)

# Calcular o percentual de acerto
conf_matrix <- prop.table(conf_matrix, margin = 1)

# Adicionar uma linha total à matriz de confusão
conf_matrix <- rbind(conf_matrix, colSums(conf_matrix))

# Imprimir a matriz de confusão
conf_matrix

##              Concluded    Active  Cancelled  Distressed
## Concluded  0.788732394 0.2112676 0.00000000 0.000000000
## Active     0.006115108 0.9654676 0.02110312 0.007314149
## Cancelled          NaN       NaN        NaN         NaN
## Distressed         NaN       NaN        NaN         NaN
##                    NaN       NaN        NaN         NaN

Curva ROC

A Curva ROC (Receiver Operating Characteristic Curve) associada ao modelo logístico mensura a capacidade de predição do modelo proposto, através das predições da sensibilidade e da especificidade. Segundo Fawcett (2006) esta técnica serve para visualizar, organizar e classificar o modelo com base na performance preditiva.

A curva ROC é produzida bi-dimensionalmente pela obtenção da relação entre a taxa dos verdadeiros positivos do modelo e da taxa dos falsos positivos preditos, como mostra a Figura a) abaixo:

Curva ROC

Desta forma, o ponto inferior esquerdo (0,0) significa que não é predita uma classificação positiva; no canto oposto do gráfico (1,1) classifica os resultados incondicionalmente positivos e; o ponto (0,1) representa uma excelente classificação. Quanto mais ao noroeste do gráfico o ponto estiver melhor, assim sendo o ponto B da Figura a) classifica melhor os resultados que o ponto C.

# Transformar o problema em vários problemas de classificação binária
# "logm1" é o objeto que representa o modelo logit multinomial
# e "m1" é o conjunto de dados correspondente

#pacote para trabalhar com a lógica de classificação

# função para ajustar modelos one-vs-rest
fit_one_vs_rest <- function(class_label) {
  binary_data <- m1
  binary_data$target <- ifelse(statusm1 == class_label, 1, 0)
  binary_model <- glm(target ~ ITm1 + km1 + tm1 + percentprivm1 + catpppm1 + sector + debtm1 + catbidm1, data = binary_data, family = binomial)
  return(binary_model)
}

# Ajuste um modelo para cada classe
class_labels <- unique(statusm1)
binary_models <- lapply(class_labels, fit_one_vs_rest)

# Calcular as métricas para cada modelo binário
roc_list <- lapply(binary_models, function(model) {
  roc_obj <- roc(statusm1, predict(model, type = "response"))
  return(roc_obj)
})

# Plotar a curva ROC
#plot(roc_list[[1]]$specificities, roc_list[[1]]$sensitivities, type = "l", col = 1, 
 #    main = "Curvas ROC para o modelo logit multinomial", xlim = c(0, 1), ylim = c(0, 1), reverse = #TRUE)

plot(roc_list[[1]], col = 1, main = "Curvas ROC para o modelo logit multinomial")
for (i in 2:length(roc_list)) {
  lines(roc_list[[i]], col = i)
}
legend("bottomright", legend = class_labels, col = 1:length(roc_list), lty = 1)

No gráfico, cada curva representa um dos quatro estados possíveis da variável dependente: ativo, cancelado, concluído ou em dificuldade. A curva ROC para o estado Active é a curva mais alta, o que significa que esse estado é o mais facilmente classificado pelo modelo. A curva ROC para o estado em Distressed é a curva mais baixa, o que significa que esse estado é o mais difícil de classificar pelo modelo.

A área sob a curva ROC (AUC) é uma medida da capacidade geral do modelo de classificação. Um AUC de 1,0 indica um modelo perfeito, enquanto um AUC de 0,5 indica um modelo aleatório. No gráfico, o AUC para o estado Active é 0,92, o que indica um bom desempenho do modelo para esse estado. O AUC para o estado em Distressed é 0,70, o que indica um desempenho menos bom do modelo para esse estado.

Em geral, o gráfico indica que o modelo logit multinomial é capaz de classificar os quatro estados da variável dependente com algum grau de sucesso. No entanto, o modelo é mais eficaz na classificação do estado Active do que do estado em Distressed.

#exp(cbind(OR=coef(logm1), confint(logm1)))

Probabilidades estimadas

“fitted(logm1)”: Este trecho extrai as probabilidades estimadas com base no modelo multinomial ajustado, armazenando essas probabilidades na variável ‘pp’.

“head(pp)”: Esta parte do código mostra as primeiras linhas dos resultados das probabilidades estimadas, exibindo os valores de probabilidade para as primeiras observações. A função ‘head()’ é usada para exibir apenas as primeiras linhas, tornando mais fácil inspecionar os resultados iniciais.

A saída do código mostra as probabilidades estimadas para as quatro categorias (‘Concluded’, ‘Active’, ‘Cancelled’, ‘Distressed’) para as primeiras observações nos dados. Cada linha corresponde a um conjunto de probabilidades estimadas para uma observação específica, indicando a probabilidade de pertencer a cada uma das quatro categorias. Isso pode ser útil para entender como o modelo está atribuindo as probabilidades para cada categoria com base nas variáveis preditoras fornecidas.

head(pp <- fitted(logm1))

##      Concluded    Active   Cancelled  Distressed
## 1 0.0005148894 0.9869965 0.009423937 0.003064715
## 2 0.0043401862 0.9815383 0.011016929 0.003104594
## 3 0.0045558431 0.9813250 0.011017776 0.003101395
## 4 0.0001990151 0.9884275 0.008374322 0.002999152
## 5 0.0002019425 0.9884246 0.008375036 0.002998412
## 6 0.0002152592 0.9883586 0.008425603 0.003000559

Referências

Fawcett, Tom. 2006. «An introduction to ROC analysis». Pattern Recognition Letters 27: 861–74. https://doi.org/10.1016/j.irbm.2014.09.001.

#GPT
#dcatstatus <- data.frame(statusm1 = rep(c("Active", "Cancelled", "Concluded", "Distressed"), each = 10), ITm1 = rep(c(0:39), times = 1))

# Criar um data.frame com as variáveis necessárias para as previsões

#dcatstatus <- data.frame(statusm1 = rep(c("Active", "Cancelled", "Concluded", "Distressed"), each = 813), ITm1 = rep(c(240:3492), times = 1), km1 = rep(c(180:3022),  times = 1), percentprivm1 = rep(c(80:100), times = 162), tm1 = rep(c(20:99), times = 41), debtm1 = rep(c(202:1321), times = 2.9), 
 #                        catbidm1 = rep(c("1","2","3","4","5","6","7","8"), each = 406), 
  #                       catpppm1 = rep(c("1","2","3","4"), each = 812.5), 
   #                      catsetorm1 = rep(c("1","2","3","4","5"), each = 650) )

#GPT length
#dcatstatus <- data.frame(
 # statusm1 = rep(c("Active", "Cancelled", "Concluded", "Distressed"), each = 813),
#  ITm1 = rep(c(0:34920), length.out = 813),
 # km1 = rep(c(180:3022), length.out = 813),
  #percentprivm1 = rep(c(20:100), length.out = 813),
  #tm1 = rep(c(1:99), length.out = 813),
  #debtm1 = rep(c(1:13210), length.out = 813),
  #catbidm1 = rep(c("1","2","3","4","5","6","7","8"), length.out = 813),
  #catpppm1 = rep(c("1","2","3","4"), length.out = 813),
  #catsetorm1 = rep(c("1","2","3","4","5"), length.out = 813)
#)



#pp.inv <- cbind(dcatstatus, predict(logm1, newdata = dcatstatus, type = "probs", se = TRUE))

#by(pp.inv[,10:13], pp.inv$sector, colMeans)
#```

#```{r}
#linvp <- melt(pp.inv, id.vars = c("statusm1", "ITm1"), value.name = "probability")

#head(linvp)

#ggplot(pp.inv, aes(x=ITm1, y = probability, colour = catstatusm1)) + geom_line() + facet_grid(variable ~., scales = "free")

1. Logit Multinomial m1 - Banco Mundial

Kassyo Costa (CAEN/UFC)

22-06-2025

Carregando bibliotecas

Definição das variáveis m1.

Variável dependente \(statusm1\) (Categórica multinomial)

Categorização dos tipos de contratos para escrever as equações como:

Definição das variáveis controle m1

Variáveis categóricas m1

(\(catpppm1\), Tipo de PPI) Os tipos de contratos para a definição das variáveis categóricas não ordenadas:

(\(catbidm1\)) Tipos de Leilão para os contratos e a definição das variáveis categóricas não ordenadas:

Definição setorial categórica

Definição de categórica para o metódo de adjudicar usado pelo governo no contrato

Categoria de Renda

Análise Fatorial

Logit Multinomial

Estimação do modelo

Testar heterocedasticidade

Método Stepwise

Modelo estimado novamente sem a variável \(k_{m1}\).

Análise da razão estatística de máxima verossimilhança (LR statistic)

*Pseudo \(R^2\)

Estatítica \(Z\)

\(p\)-valor

Razão de chances

Matriz de Confusão

Curva ROC

Probabilidades estimadas

Referências