Programação Linear no R.

A programação linear possui muitas aplicações nas Ciências Sociais e Exatas.

Nesse post mostrarei como solucionar problemas de programação linear no R por meio da biblioteca lpSolveAPI.

Considere o seguinte problema de programação linear:

O primeiro passo no R para a resolução desse problema é invocar a biblioteca lpSolveAPI, para isso utilizamos o seguinte comando:

#Chama a biblioteca
library(lpSolveAPI)

Como o problema possui somente duas variáveis vamos criá-lo com o nome de modelo.lp1 inicialmente sem nenhuma restrição:

#Define a criação de um modelo com 0 restrições e 2 variáveis
modelo.lp1 <- make.lp(0, 2)

#Dá o nome ao problema de programação linear
name.lp(modelo.lp1, "Exemplo 1 - Aula de MMQD 1")

Note que o comando make.lp recebe dois argumentos, o primeiro informa quantas restrições existem no problema e o segundo quantas variáveis. Inicialmente, colocamos zero restrições pois iremos adicioná-las dinamicamente nos próximos passos. Já o comando name.lp recebe dois argumentos, o primeiro argumento é o objeto modelo.lp1 que nós criamos e o segundo é o nome (ou descrição) do problema. Uma vez definida as principais características do problema de programação linear precisamos definir o domínio das variáveis e o tipo de problema (maximização ou minimização). Como nosso problema é um problema de maximação, escrevemos:

#Define as características do modelo
lp.control(modelo.lp1, sense="max")

Nesse comando, atribuímos ao modelo.lp1 o sentido de Maximização (isso é sense=''max'') para a função objetivo. Caso desejássemos Minimizar ao invés de maximizar escreveríamos sense=''min''. Como a função objetivo é escrita na forma $Minimize: Z=5x_{1} +6x_{2}$ devemos informar ao R que os coeficientes da função objetivo são, $5$ e $6$ respectivamente. Isso é feito por meio da seguinte sintaxe:

#Define a função objetivo
set.objfn(modelo.lp1, c(5,6))

Em seguida, como o problema pode assumir valores $x_{1}\geq 0,x_{2}\geq 0$ atribuímos os seguintes limites para as variáveis de decisão:

#Define os limites da região factível
set.bounds(modelo.lp1, lower = c(0,0), upper = c(Inf, Inf))

#Tipo das variáveis de decisão
set.type(modelo.lp1, c(1,2), type = c("real"))

No comando set.bounds informamos que a primeira variável possui limite inferior (lower) igual a zero e limite superior (upper) igual a infinito (Inf). O mesmo fazemos para a segunda variável, nesse caso dizemos que o limite inferior (lower) é igual a zero e o limite superior (upper) é igual a infinito (Inf). A ordem aqui importa, ou seja, caso a primeira variável estivesse contida entre zero e infinito e a segunda variável estivesse contida entre os valores 1 e 4 deveríamos escrever set.bounds(modelo.lp1, lower = c(0,1), upper = c(Inf, 4)). O comando set.type define o tipo de variável para cada uma das variáveis de decisão. É possível escolher as seguintes opções:

1. integer - A variável só pode assumir valores inteiros.
2. binary - A variável só pode assumir valores binários, isso é, 1 ou zero.
3. real - A variável pode assumir valores nos reais.

O comando set.type informa ao objeto modelo.lp1 que as variáveis $x_{1}$ e $x_{2}$ (c(1,2)) são números reais. Mas o comando anterior, set.bounds diz ao R que as variáveis devem estar entre 0 e $\infty$, portanto, $x_{1}\geq 0$ e $x_{2}\geq 0$. O passo 3 é responsável pela adição das restrições no problema de Programação Linear. Vamos então adicionar uma restrição de cada vez. A primeira restrição é dada por: $x_{1} \le 6$ a qual pode ser escrita na forma $ 1x_{1}+0x_{2} \le 6$. Assim, a adição dessa restrição no objeto modelo.lp1 é feita da seguinte forma:

#Define os coeficientes da primeira restrição
coef1 <- c(1,0)

#Adiciona a restrição
add.constraint(modelo.lp1, coef1, "<=", 6)

Para a segunda restrição $2x_{2} \le 12\Rightarrow 0x_{1}+2x_{2} \le 12$, a adição é feita da seguinte forma:

#Define os coeficientes da segunda restrição
coef2 <- c(0,2)

#Adiciona a restrição
add.constraint(modelo.lp1, coef2, "<=", 12)

Finalmente, a última restrição $3x_{1}+2x_{2}\le 18$ é adicionada fazendo-se:

#Define os coeficientes da terceira restrição
coef3 <- c(3,2)

#Adiciona a restrição
add.constraint(modelo.lp1, coef3, "<=", 18)

Uma vez que o problema tenha sido configurado, todas as informações associadas a ele estão armazenadas no objeto modelo.lp1. Podemos visualizar essas informações fazendo:

#Mostra as informações do modelo
print(modelo.lp1)

Podemos também visualizar esse problema, uma vez que é um problema bidimensional, isso é, contém somente duas variáveis de decisão:

#Plota a região factível
plot(modelo.lp1)

Finalmente, para resolvê-lo, fazemos:

#Resolve o problema primal
solve(modelo.lp1)
resultado<-get.primal.solution(modelo.lp1)
print(resultado)

Entretanto, como os resultados não estão em um formato muito explicativo, vamos transformá-los em uma tabela (isso é um data.frame) e colocar títulos nas linhas e nas colunas:

#Transforma o resultado em uma tabela
solucao<-as.data.frame(resultado)

#Coloca os nomes nas colunas e nas linhas da tabela
names(solucao)<-c("Valores")
rownames(solucao)<-c("Função Objetivo", "Variável de folga 1","Variável de folga 2",
"Variável de folga 3","Solução X1", "Solução X2")

#Mostra os resultados com os nomes
print(solucao)

Causalidade e Correlação

Causalidade é a relação entre um evento (a causa) e um segundo evento (o efeito), em que o segundo acontecimento é entendida como uma consequência do primeiro.

No uso comum, a causalidade é também a relação entre um conjunto de fatores (causas) e um fenômeno (o efeito). Qualquer coisa que afete um efeito, é denominada fator desse efeito. Um fator direto é um fator que afeta diretamente o efeito, isto é, sem quaisquer fatores intervenientes.

Compreender a relação causa-efeitoentre as variáveis ​​é de interesse primordial em muitos campos da ciência. Normalmente, a intervenção experimental é utilizada para avaliar estas relações.

Entre os métodos comuns para a avaliação da cause e efeito destacam-se:

1. Variáveis Instrumentais.
2. Equações estruturais.
3. Redes Bayesianas.
4. Modelos Gráficos.

Pacote pcalg

Suponha que temos um sistema descrito por algumas variáveis ​​e muitas observações obtidas deste sistema. Além disso, suponha que seja plausível a não existência de variáveis faltantes (omissão de variáveis) e também que no sistema não apresenta loops de feedback do sistema causal subjacente.

A estrutura causal de um sistema deste tipo pode ser convenientemente representado por um gráfico acíclico dirigido (DAG - Directed Acyclic Graph), em que cada nó representa uma variável e cada aresta representa uma causa direta.

Por exemplo, suponha o seguinte grafo:

Nesse exemplo, a variável 1 causa a variável 2 diretamente, ou seja é um fator direto para a variável 2 mas é um fator interveniente da variável 5.

Usualmente, não se conhece a relação entre as variáveis a não ser que algum modelo teórico seja desenvolvido para explicar a relação entre as variáveis, e mesmo nesse caso o modelo necessita ser validado.

O caso mais comum ocorre quando não há qualquer informação quanto ao comportamento causal dessas variáveis, nesse caso o objetivo é estimar essas relações e apresentar o grafo causal estimado.

Os modelos implementados no pacote pcalg são capazes de estimar o grafo causal das variáveis de uma base de dados mesmo que nenhum modelo seja conhecido a priori.

Os métodos implementados são o Algoritmo PC (Spirtes et al., 2000), Algoritmo FCI (Spirtes et al. ,1999), Algoritmo RFCI (Colombo et al., 2012) e o Método IDA (Maathuis et al., 2009).

Alguns pressupostos para cada um dos algoritmos implementados no pacote pcalg são apresentados abaixo:

• Algoritmo PC: Nenhuma variável oculta existe ou deixou de ser selecionada.
• Algoritmo FCI: Permite a existência de variável oculta ou omissão de variáveis.
• Algoritmo RFCI: Permite a existência de variável oculta ou omissão de variáveis.
• Método IDA: Nenhuma variável oculta existe ou deixou de ser selecionada.

Prática no R.

Antes de instalar o pacote pcalg no R é necessário instalar dois pacotes:


Para a representação de gráficos, o pacote pcalg utiliza dois outros pacotes chamados RBGL e o pacote graph.

Estes pacotes não estão disponíveis no CRAN, mas estão no outro repositório do software R, BioConductor.

Para instalá-los, siga as instruções abaixo:

#Instala os pacotes necessários.
source("http://bioconductor.org/biocLite.R") 
biocLite("RBGL")
biocLite("Rgraphviz")

Após a instalação desses pacotes, instale normalmente o pacote pcalg. Como exemplo utilizaremos dados simulados pelos autores do pacote pcalg.

Esses dados foram construídos de modo a possuírem a estrutura gráfica apresentada anteriormente. Aqui, realizaremos estimativas com base nos dados para avaliar se os métodos gráficos implementados no pcalg conseguem identificar corretamente a estrutura gráfica original.

Devemos inicialmente habilitar o pacote pcalg e em seguida ler o conjunto de dados que desejamos explorar:

#Habilita o pacote pcal
library("pcalg")

#Lê os dados simulados
dados.df<-read.csv("http://dl.dropbox.com/u/36068691/DadosPCALG.csv",sep=" ")

O primeiro método que utilizaremos é o Algoritmo PC:

#Gera as Estatísticas Suficientes
estatisticas  <-  list(C  =  cor(dados.df),  n  =  nrow(dados.df))

#Executa o Algoritmo PC
pc.fit  <-  pc(estatisticas,  indepTest  =  gaussCItest, p  =  ncol(dados.df),  alpha  =  0.01)

#Plota os resultados
plot(pc.fit,  main  =  "Algoritmo PC.")

Na primeira parte do comando a lista estatisticas armazena a Matriz de Correlações e o tamanho da amostra (número de observações na base). Essas informações são suficientes para a execução do Algoritmo PC. Em seguida a função pc(.) executa o Algoritmo PC. O comando "gaussCItest" testa a independência condicional assumindo que as variáveis possuem distribuição normal a um nível de significância de 0.01. O resultado obtido fornece o seguinte grafo:

Note que o grafo gerado está muito próximo da estrutura verdadeira, os desvios ocorrem devido ao tamanho amostral e flutuações estocásticas. Por exemplo, a relação entre as variáveis 1 e 6 foi apresentada corretamente, mas a relação entre as variáveis 1 e 2 apresenta uma relação bidirecional, o que não é verdade para os dados simulados. O próximo algoritmo a ser utilizado é o Algoritmo FCI. O Algoritmo FCI é uma generalização do Algoritmo PC, no sentido de permitir a existência de muitas possíveis variáveis latentes (omitidas) no modelo. A execução do Algoritmo FCI é similar ao do Algoritmo PC:

#Executa o Algoritmo FCI
fci.fit  <-  fci(estatisticas,  indepTest  =  gaussCItest, p  =  ncol(dados.df),  alpha  =  0.01)

#Plota os resultados
plot(fci.fit,  main  =  "Algoritmo FCI.")

O grafo obtido usando esse algoritmo foi:

Note que nesse caso, algumas relações foram identificadas mas a direção não pode ser estimada. Por exemplo, a variável 1 causa a variável 6, mas a relação entre a variável 1 e a variável 2 não pode ser identificada... A única informação que temos é que possivelmente as variáveis 1 e 2 se relacionam mas não sabemos como. Para o Algoritmo RFCI utilizamos a seguinte sintaxe:

#Executa o Algoritmo RFCI
rfci.fit  <-  rfci(estatisticas,  indepTest  =  gaussCItest, p  =  ncol(dados.df),  alpha  =  0.01)

#Plota os resultados
plot(rfci.fit,  main  =  "Algoritmo RFCI.")

O Algoritmo RFCI é similar ao Algoritmo FCI. No entanto, esse é mais rápido em sua execução. O grafo produzido, nesse caso, é idêntico ao Algoritmo FCI:

Por fim, o último algoritmo é o Método IDA. Esse método, diferente dos demais algoritmos mensura o grau de relação entre as variáveis. Por exemplo, se desejamos saber qual o grau de efeito entre as variáveis 1 e 6 (variável 1 causando a variável 6) utilizamos o seguinte código:

#Executa o Método IDA para o Algoritmo PC
ida(1, 6, cov(dados.df), pc.fit@graph)

O comando retorna os seguintes valores [0.7536376 0.5487757]. Isso significa que a magnitude do efeito da variável 1 sobre a variável 6 é algum valor nesse intervalo. Uma vez que ambos os valores são maiores do que zero, podemos concluir que a variável 1 apresenta um efeito positivo causal sobre a variável 6. Esses valores refletem o efeito do aumento em uma unidade na variável 1 em relação a variável 6. Em alguns casos quando se conhece determinadas relações entre as variáveis, é possível fixar relações e permitir que os algoritmos busquem apenas as relações de causalidade desconhecidas, nesses casos utiliza-se os comandos fixedGaps e fixedEdges.

Modelos de mediação.

Em muitos campos da ciência o objetivo dos pesquisadores não é apenas estimar o efeito causal de um tratamento, mas também a compreensão do processo no qual o tratamento afeta o resultado de maneira causal.

A análise de mediação causal é frequentemente usada para avaliar os potenciais mecanismos causais.

Estudiosos de diversas disciplinas estão cada vez mais interessados ​​em identificar mecanismos causais, indo além da estimativa dos efeitos causais e explorando por completo o modelo causal.

Uma vez que certas variáveis foram identificadas como responsáveis pelo efeito causal associado a um determinado resultado, o próximo passo é entender como essas variáveis ​​exercem influência.

Nesses casos, o procedimento padrão para analisar os mecanismos causais na pesquisa aplicada é chamada análise de mediação. Nessas análises um conjunto de modelos de regressão, são ajustados e, em seguida, as estimativas dos "efeitos de mediação" são calculados a partir dos modelos ajustados (por exemplo, Haavelmo 1943; Baron e Kenny 1986; Shadish, Cook e Campbell 2001; MacKinnon 2008).

Pacote mediation.

O pacote de mediation do R permite aos usuários:

1. Investigar o papel dos mecanismos causais utilizando diferentes tipos de dados e modelos estatísticos.
2. Explorar como os resultados mudam quando os pressupostos são relaxados (análise de sensibilidade).
3. Calcular medidas de interesse em diversos projetos de pesquisa.

A prática corrente na análise de mediação atualmente são inferências baseadas em modelos. Em um delineamento experimental, a variável tratamento é randomizada e as variáveis de mediação e os resultados são observados.

Imai et ai. (2010) mostram que uma gama de modelos paramétricos e semi-paramétrico podem ser utilizados para estimar o efeito médio causal média e mediação causal e outras quantias de interesse.

Aplicação usando o pacote mediation.

Para ilustrar o procedimento utilizaremos como exemplo os dados obtidos por Brader, Valentino e Suhat (2008).

Brader et al. (2008) conduziram um experimento casualizado onde os indivíduos foram expostos a diferentes histórias sobre a imigração. O objetivo era investigar como o seu enquadramento os influencia as atitudes e comportamentos políticos em relação à política de imigração.

#Habilita o pacote mediation
library("mediation")

#Lê os dados
dados.df<-read.csv("http://dl.dropbox.com/u/36068691/dadosMediacao.csv")

Os autores postularam a ansiedade (variável emo) como variável mediadora do efeito causal para o enquadramento na opinião pública. O primeiro passo é ajustar o modelo de mediação para a medida de ansiedade (variável emo), essa é modelada como uma função da variável tratamento (treat) e co-variáveis ​​de pré-tratamento como (idade - age, educ, gênero - gender e renda - income).

# Modelo de mediação
med.fit <- lm(emo ~ treat + age + educ + gender + income, data = dados.df)

Em seguida, modelamos a variável resultado, a qual é uma variável binária que indica se o participante concordou em enviar uma carta sobre a política de imigração ao seu membro do Congresso (variável cong_mesg). As variáveis ​​explicativas do modelo incluem a variável de mediação, a variável de tratamento, e o mesmo conjunto de variáveis de pré-tratamento ​​como as utilizadas no modelo de mediação.

# Modelo para o resultado
out.fit <- glm(cong_mesg ~ emo + treat + age + educ + gender + income, data = dados.df, family = binomial("probit"))

Neste exemplo, espera-se que o tratamento aumente a resposta emocional dos entrevistados, que por sua vez é postulado fazer com que os respondentes sejam mais propensos a enviar uma carta ao seu membro do Congresso. Inicialmente usamos um modelo de regressão linear e regressão de probit para os modelos de mediação e o modelo para o resultado, respectivamente. Vamos agora usar o modelo de mediação para estimar os efeitos médios causais diretos (ACME - Average Causal Mediation effects). Para isso, basta executar o seguinte código:

#Habilita o pacote sandwich
library(sandwich)

#Estima os efeitos médios
med.out <- mediate(med.fit, out.fit, treat = "treat", mediator = "emo", robustSE = TRUE)

#Apresenta as estimativas
summary(med.out)

Como argumentos para esta função, é necessário especificar os modelos (neste caso, med.fit e out.fit), bem como o nome da variável de tratamento e da variável de mediação, que são representados como treat = "treat" e mediator = "emo", respectivamente. Além disso, usamos a matriz de variâncias e covariâncias robustas para a heterocedasticidade oriunda do pacote de sandwich. Ao executar o comando, o seguinte resultado surge:

Nesse exemplo, utilizou-se 1000 simulações para que o erro-padrão das estimativas fossem obtidos. Neste exemplo, os efeitos médios causais diretos (ACME - Average Causal Mediation effects) estimados são estatisticamente significantes e portanto, diferentes de zero, mas as estimativas dos efeitos diretos para médio e total não são. Em outras palavras, os resultados sugerem que o tratamento no ensaio pode ter aumentado a resposta emocional, que por sua vez tornou os indivíduos mais propensos a enviar uma mensagem ao seu congressista. Aqui, uma vez que a variável resultado é binária todos os efeitos estimados são expressos como uma alteração na probabilidade do respondente enviar uma mensagem ao Congresso. O pacote mediation apresenta outras opções como a possibilidade de análises gráficas, análises por segmentos e múltiplas entradas. Para maiores detalhes veja Tingley et. al (2012).

Otimização de portfólio.

O pacote do R denominado Tawny fornece uma maneira simples para otimizar carteiras de investimento de forma a minimizar o risco em uma carteira.

Este otimizador pode executar a otimização de carteiras para diversos períodos temporais simultaneamente.

A ideia é: "não colocar todos os seus ovos em uma mesma cesta".

Para fins de ilustração, considere um subconjunto de ativos do S&P 500. Abaixo estão os códigos em R para otimizar a carteira de investimentos:

#Habilita a biblioteca Tawny 
library(tawny)

#Lê o subconjunto de dados do S&P 500
data(sp500.subset)
dados <- create(TawnyPortfolio, sp500.subset, window=190)

# Otimizando por meio de uma janela temporal de tamanho 190
pesos <- optimizePortfolio(dados, create(SampleFilter) )