rm(list=ls())
setwd("C:/DAVE2/CIDE/Doc/Econometria II/Lab_4")
options(scipen=999) 

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.0.5
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.0.5
## Warning: package 'tibble' was built under R version 4.0.5
## Warning: package 'tidyr' was built under R version 4.0.5
## Warning: package 'readr' was built under R version 4.0.5
## Warning: package 'dplyr' was built under R version 4.0.5
## Warning: package 'forcats' was built under R version 4.0.5
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(janitor)
## Warning: package 'janitor' was built under R version 4.0.5
## 
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(modelsummary)
## Warning: package 'modelsummary' was built under R version 4.0.5
library(nnet)
library(mlogit)
## Loading required package: dfidx
## Warning: package 'dfidx' was built under R version 4.0.5
## 
## Attaching package: 'dfidx'
## The following object is masked from 'package:stats':
## 
##     filter
library(sandwich)
## Warning: package 'sandwich' was built under R version 4.0.5

Modelos multinomiales

#Preliminares
#rm(list=ls())
#Cargamos paquetes necesarios
pacman::p_load(ggplot2,nnet,mlogit, foreign,reshape2,xtable,stargazer,tibble,knitr,dplyr,stats,margins,pscl,devtools,PoEdata,MCMCpack,survival,Epi,mclogit,mnlogit,gmnl)
#Los ajustes de modelos logit multinomiales: mlogit y nnet
# Si presentas problemas para correr el archivo .rmd ve a file>save with encoding>UTF-8

MODELOS MULTINOMIALES

Un modelo con solo variables específicas individuales a veces se denomina modelo logit multinomial, uno con solo variables específicas alternativas un modelo logit condicional y otro con ambos tipos de variables un modelo logit mixto.

Logit multinomial

Ejemplo 1. Las elecciones ocupacionales de las personas pueden verse influidas por las ocupaciones de sus padres y su propio nivel educativo. Podemos estudiar la relación de la elección de ocupación de uno con el nivel de educación y la ocupación del padre. Las opciones ocupacionales serán la variable de resultado que consta de categorías de ocupaciones.

Ejemplo 2. Un biólogo puede estar interesado en las elecciones de alimentos que hacen los caimanes. Los caimanes adultos pueden tener preferencias diferentes a las de los jóvenes. La variable de resultado aquí serán los tipos de comida, y las variables predictoras podrían ser el tamaño de los caimanes y otras variables ambientales.

Ejemplo 3. Los estudiantes que ingresan a la escuela secundaria hacen elecciones de programas entre programa general, programa vocacional y programa académico. Su elección puede ser modelada usando su puntuación de escritura y su estatus económico social.

Ejemplo Clase

Es importante tener los paquetes: * nnet: más usado para logit multinomial, pero también para muchas otras cosas más * mlogit: un poco más confuso de usar, más usado para el logit condicional * sandwich: para trabajar con matríces de varianzas y construir versiones robustas

Importmaos la based dedaos que utilizamos para estimar los modelos probit y logit, son los datos de Herriges & Kling (1999) usados y provistos por Cameron & Trivedi (2005):

data_fishing<-read_csv(
  "./fishing_data.csv",
  locale = locale(encoding = "latin1")) %>% 
  clean_names() %>% 
  mutate(income=income/1000)
## Rows: 1182 Columns: 17
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## dbl (17): client, mode, price, crate, dbeach, dpier, dprivate, dcharter, pbe...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

Analizamos la base de datos:

head(data_fishing)
## # A tibble: 6 x 17
##   client  mode price  crate dbeach dpier dprivate dcharter pbeach ppier pprivate
##    <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>  <dbl> <dbl>    <dbl>
## 1      1     4 183.  0.539       0     0        0        1  158.  158.     158. 
## 2      2     4  34.5 0.467       0     0        0        1   15.1  15.1     10.5
## 3      3     3  24.3 0.241       0     0        1        0  162.  162.      24.3
## 4      4     2  15.1 0.0789      0     1        0        0   15.1  15.1     55.9
## 5      5     3  41.5 0.108       0     0        1        0  107.  107.      41.5
## 6      6     4  63.9 0.398       0     0        0        1  192.  192.      28.9
## # ... with 6 more variables: pcharter <dbl>, qbeach <dbl>, qpier <dbl>,
## #   qprivate <dbl>, qcharter <dbl>, income <dbl>

Vamos a ver la estruictura de datos:

str(data_fishing)
## spec_tbl_df [1,182 x 17] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ client  : num [1:1182] 1 2 3 4 5 6 7 8 9 10 ...
##  $ mode    : num [1:1182] 4 4 3 2 3 4 1 4 3 3 ...
##  $ price   : num [1:1182] 182.9 34.5 24.3 15.1 41.5 ...
##  $ crate   : num [1:1182] 0.5391 0.4671 0.2413 0.0789 0.1082 ...
##  $ dbeach  : num [1:1182] 0 0 0 0 0 0 1 0 0 0 ...
##  $ dpier   : num [1:1182] 0 0 0 1 0 0 0 0 0 0 ...
##  $ dprivate: num [1:1182] 0 0 1 0 1 0 0 0 1 1 ...
##  $ dcharter: num [1:1182] 1 1 0 0 0 1 0 1 0 0 ...
##  $ pbeach  : num [1:1182] 157.9 15.1 161.9 15.1 106.9 ...
##  $ ppier   : num [1:1182] 157.9 15.1 161.9 15.1 106.9 ...
##  $ pprivate: num [1:1182] 157.9 10.5 24.3 55.9 41.5 ...
##  $ pcharter: num [1:1182] 182.9 34.5 59.3 84.9 71 ...
##  $ qbeach  : num [1:1182] 0.0678 0.1049 0.5333 0.0678 0.0678 ...
##  $ qpier   : num [1:1182] 0.0503 0.0451 0.4522 0.0789 0.0503 ...
##  $ qprivate: num [1:1182] 0.26 0.157 0.241 0.164 0.108 ...
##  $ qcharter: num [1:1182] 0.539 0.467 1.027 0.539 0.324 ...
##  $ income  : num [1:1182] 7.08 1.25 3.75 2.08 4.58 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   client = col_double(),
##   ..   mode = col_double(),
##   ..   price = col_double(),
##   ..   crate = col_double(),
##   ..   dbeach = col_double(),
##   ..   dpier = col_double(),
##   ..   dprivate = col_double(),
##   ..   dcharter = col_double(),
##   ..   pbeach = col_double(),
##   ..   ppier = col_double(),
##   ..   pprivate = col_double(),
##   ..   pcharter = col_double(),
##   ..   qbeach = col_double(),
##   ..   qpier = col_double(),
##   ..   qprivate = col_double(),
##   ..   qcharter = col_double(),
##   ..   income = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

Ahora vamos a estimar el modelo:

mmultilogit_1 <- multinom(
  mode ~ income,
  data = data_fishing)
## # weights:  12 (6 variable)
## initial  value 1638.599935 
## iter  10 value 1477.505846
## final  value 1477.150570 
## converged

Una vez estimado podemos obtener un resumen de las estimaciones:

summary(mmultilogit_1)
## Call:
## multinom(formula = mode ~ income, data = data_fishing)
## 
## Coefficients:
##   (Intercept)      income
## 2   0.8141701 -0.14340453
## 3   0.7389569  0.09190030
## 4   1.3413284 -0.03164588
## 
## Std. Errors:
##   (Intercept)     income
## 2   0.2286325 0.05328822
## 3   0.1967314 0.04066362
## 4   0.1945173 0.04184622
## 
## Residual Deviance: 2954.301 
## AIC: 2966.301

La tabla 15.2 dice que todos los coeficientes, excepto \(\beta_{I4}\) son estadísticamente significativos al 5%.

Para poder calcular los p values asociados tenremos que utilizar el siguiente comando:

z <- summary(mmultilogit_1)$coefficients/summary(mmultilogit_1)$standard.errors
pv <- (1 - pnorm(abs(z)))*2 
pv
##            (Intercept)      income
## 2 0.000369384199215173 0.007121444
## 3 0.000172532839256112 0.023820471
## 4 0.000000000005360823 0.449503854

El paquete usado para estimar los logit multinomial es nnet, que es una abreviatura para neural networks

Las redes neuronales, una aplicación en constante desarrollo del aprendizaje automatizado pueden emplearse, entre muchas otras cosas, para problemas de clasificación. Pueden encontrar màs informaciòn al respecto en este link. Finalmente podemos considerar nuestro modelos multinomiales como problemas de clasificaciòn.

Ejemplo 2

library("sas7bdat")
## Warning: package 'sas7bdat' was built under R version 4.0.3

Importamos los datos:

math =read.table("http://www.utstat.utoronto.ca/~brunner/data/legal/mathcat.data.txt")
head(math)
##   hsgpa hsengl hscalc   course passed outcome
## 1  78.0     80    Yes Mainstrm     No  Failed
## 2  66.0     75    Yes Mainstrm    Yes  Passed
## 3  80.2     70    Yes Mainstrm    Yes  Passed
## 4  81.7     67    Yes Mainstrm    Yes  Passed
## 5  86.8     80    Yes Mainstrm    Yes  Passed
## 6  76.7     75    Yes Mainstrm    Yes  Passed

Modificamos los datos para que la funcion de regresion logistica multinomial pueda procesarlos. Para hacer esto, necesitamos expandir la variable de resultado (y) como lo hariamos para la codificacion ficticia de una variable categ?rica para incluirla en la regresion multiple estandar.

math0 = math[,c(1,5)]; 
head(math0)
##   hsgpa passed
## 1  78.0     No
## 2  66.0    Yes
## 3  80.2    Yes
## 4  81.7    Yes
## 5  86.8    Yes
## 6  76.7    Yes
long0 = mlogit.data(math0,shape="wide",choice="passed")
head(long0)
## ~~~~~~~
##  first 10 observations out of 788 
## ~~~~~~~
##    hsgpa passed chid alt   idx
## 1   78.0   TRUE    1  No  1:No
## 2   78.0  FALSE    1 Yes 1:Yes
## 3   66.0  FALSE    2  No  2:No
## 4   66.0   TRUE    2 Yes 2:Yes
## 5   80.2  FALSE    3  No  3:No
## 6   80.2   TRUE    3 Yes 3:Yes
## 7   81.7  FALSE    4  No  4:No
## 8   81.7   TRUE    4 Yes 4:Yes
## 9   86.8  FALSE    5  No  5:No
## 10  86.8   TRUE    5 Yes 5:Yes
## 
## ~~~ indexes ~~~~
##    chid alt
## 1     1  No
## 2     1 Yes
## 3     2  No
## 4     2 Yes
## 5     3  No
## 6     3 Yes
## 7     4  No
## 8     4 Yes
## 9     5  No
## 10    5 Yes
## indexes:  1, 2
simple0 = mlogit(passed ~ 0 | hsgpa, data=long0)
simple0 = mlogit(passed ~ 0 | hsgpa, data=long0); summary(simple0)
## 
## Call:
## mlogit(formula = passed ~ 0 | hsgpa, data = long0, method = "nr")
## 
## Frequencies of alternatives:choice
##      No     Yes 
## 0.40102 0.59898 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.000119 
## successive function values within tolerance limits 
## 
## Coefficients :
##                   Estimate Std. Error z-value            Pr(>|z|)    
## (Intercept):Yes -15.210112   1.998398 -7.6112 0.00000000000002709 ***
## hsgpa:Yes         0.197734   0.025486  7.7587 0.00000000000000866 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -221.72
## McFadden R^2:  0.16436 
## Likelihood ratio test : chisq = 87.221 (p.value = < 0.000000000000000222)

Ahora implementamos un modelo multinomial

outcome = as.character(math$outcome)
for(j in 1:length(outcome))
   {if(outcome[j]=='Disappeared') outcome[j]='Gone'}
math$outcome = factor(outcome); table(outcome)
## outcome
## Failed   Gone Passed 
##     61     97    236

Seleccionamos los datos:

math1 = math[,c(1,6)] #datos de cols 1 al 6
long1 = mlogit.data(math1,shape="wide",choice="outcome")

Que tenemos en long1?

head(long1)
## ~~~~~~~
##  first 10 observations out of 1182 
## ~~~~~~~
##    hsgpa outcome chid    alt    idx
## 1   78.0    TRUE    1 Failed 1:iled
## 2   78.0   FALSE    1   Gone 1:Gone
## 3   78.0   FALSE    1 Passed 1:ssed
## 4   66.0   FALSE    2 Failed 2:iled
## 5   66.0   FALSE    2   Gone 2:Gone
## 6   66.0    TRUE    2 Passed 2:ssed
## 7   80.2   FALSE    3 Failed 3:iled
## 8   80.2   FALSE    3   Gone 3:Gone
## 9   80.2    TRUE    3 Passed 3:ssed
## 10  81.7   FALSE    4 Failed 4:iled
## 
## ~~~ indexes ~~~~
##    chid    alt
## 1     1 Failed
## 2     1   Gone
## 3     1 Passed
## 4     2 Failed
## 5     2   Gone
## 6     2 Passed
## 7     3 Failed
## 8     3   Gone
## 9     3 Passed
## 10    4 Failed
## indexes:  1, 2

Que tenemos en math?

head(math)
##   hsgpa hsengl hscalc   course passed outcome
## 1  78.0     80    Yes Mainstrm     No  Failed
## 2  66.0     75    Yes Mainstrm    Yes  Passed
## 3  80.2     70    Yes Mainstrm    Yes  Passed
## 4  81.7     67    Yes Mainstrm    Yes  Passed
## 5  86.8     80    Yes Mainstrm    Yes  Passed
## 6  76.7     75    Yes Mainstrm    Yes  Passed

Hacemos el ajuste multinomial

simple1 = mlogit(outcome ~ 0 | hsgpa, data=long1)
summary(simple1)
## 
## Call:
## mlogit(formula = outcome ~ 0 | hsgpa, data = long1, method = "nr")
## 
## Frequencies of alternatives:choice
##  Failed    Gone  Passed 
## 0.15482 0.24619 0.59898 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 1.09E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##                      Estimate Std. Error z-value      Pr(>|z|)    
## (Intercept):Gone     1.904226   2.744979  0.6937        0.4879    
## (Intercept):Passed -13.393056   2.570453 -5.2104 0.00000018845 ***
## hsgpa:Gone          -0.018816   0.035775 -0.5260        0.5989    
## hsgpa:Passed         0.186437   0.033018  5.6465 0.00000001637 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -326.96
## McFadden R^2:  0.11801 
## Likelihood ratio test : chisq = 87.497 (p.value = < 0.000000000000000222)

Como podemos evalua la probabilidad estimada para un estudiante con HSGPA = 90

betahat1 = simple1$coefficients
betahat1
##   (Intercept):Gone (Intercept):Passed         hsgpa:Gone       hsgpa:Passed 
##         1.90422575       -13.39305637        -0.01881621         0.18643711 
## attr(,"names.sup.coef")
## character(0)
## attr(,"fixed")
##   (Intercept):Gone (Intercept):Passed         hsgpa:Gone       hsgpa:Passed 
##              FALSE              FALSE              FALSE              FALSE 
## attr(,"sup")
## character(0)

Y

gpa = 90
L1 = betahat1[1] + betahat1[3]*gpa # Gone
L2 = betahat1[2] + betahat1[4]*gpa # Passed
denom = 1+exp(L1)+exp(L2)
pihat1 = exp(L1)/denom # Gone
pihat2 = exp(L2)/denom # Passed
pihat3 = 1/denom # Failed
rbind(pihat1,pihat2,pihat3)
##        (Intercept):Gone
## pihat1       0.03883621
## pihat2       0.92970789
## pihat3       0.03145590

Podemos tener un modelo mas completo:

long = mlogit.data(math[,c(1:4,6)],shape="wide",choice="outcome")
fullmod = mlogit(outcome ~ 0 | hsgpa+hsengl+hscalc+course, data=long)
summary(fullmod)
## 
## Call:
## mlogit(formula = outcome ~ 0 | hsgpa + hsengl + hscalc + course, 
##     data = long, method = "nr")
## 
## Frequencies of alternatives:choice
##  Failed    Gone  Passed 
## 0.15482 0.24619 0.59898 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.000216 
## successive function values within tolerance limits 
## 
## Coefficients :
##                          Estimate  Std. Error z-value      Pr(>|z|)    
## (Intercept):Gone        2.5734410   2.8288386  0.9097       0.36297    
## (Intercept):Passed    -14.0411854   2.7005870 -5.1993 0.00000020003 ***
## hsgpa:Gone             -0.0079779   0.0413277 -0.1930       0.84693    
## hsgpa:Passed            0.2157706   0.0382179  5.6458 0.00000001644 ***
## hsengl:Gone            -0.0067241   0.0251049 -0.2678       0.78882    
## hsengl:Passed          -0.0399811   0.0228733 -1.7479       0.08047 .  
## hscalcYes:Gone         -0.3902775   0.6742796 -0.5788       0.56272    
## hscalcYes:Passed        1.0009683   0.8215247  1.2184       0.22306    
## courseElite:Gone       -2.0666545   0.9836801 -2.1009       0.03565 *  
## courseElite:Passed      0.6032839   0.8044316  0.7500       0.45328    
## courseMainstrm:Gone    -0.6834686   0.5560854 -1.2291       0.21905    
## courseMainstrm:Passed   0.4086564   0.6339142  0.6447       0.51915    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -312.26
## McFadden R^2:  0.15766 
## Likelihood ratio test : chisq = 116.89 (p.value = < 0.000000000000000222)

Ejemplo 3

Los estudiantes que ingresan a la escuela secundaria hacen elecciones de programas entre programa general, programa vocacional y programa academico. Su eleccion puede ser modelada usando su puntuacion de escritura y su nivel economico social.

Datos

ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
head(ml)
##    id female    ses schtyp     prog read write math science socst       honors
## 1  45 female    low public vocation   34    35   41      29    26 not enrolled
## 2 108   male middle public  general   34    33   41      36    36 not enrolled
## 3  15   male   high public vocation   39    39   44      26    42 not enrolled
## 4  67   male    low public vocation   37    37   42      33    32 not enrolled
## 5 153   male middle public vocation   39    31   40      39    51 not enrolled
## 6  51 female   high public  general   42    36   42      31    39 not enrolled
##   awards cid
## 1      0   1
## 2      0   1
## 3      0   1
## 4      0   1
## 5      0   1
## 6      0   1

Variables:

Variable categorica multinomial:

prog=program type (general program, vocational program and academic program)

Predictoras:

ses= social economic status

write=writing score(variable continua)

#Estadisticas descriptivas
with(ml, table(ses, prog))
##         prog
## ses      general academic vocation
##   low         16       19       12
##   middle      20       44       31
##   high         9       42        7
with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
##                 M       SD
## general  51.33333 9.397775
## academic 56.25714 7.943343
## vocation 46.76000 9.318754
#ajuste 1

test1<- multinom(prog ~ ses + write, data = ml)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.985215
## final  value 179.981726 
## converged
summary(test1) #salida para cada una de la ecuaciones
## Call:
## multinom(formula = prog ~ ses + write, data = ml)
## 
## Coefficients:
##          (Intercept) sesmiddle   seshigh       write
## academic   -2.851973 0.5332914 1.1628257  0.05792480
## vocation    2.366097 0.8246384 0.1802176 -0.05567514
## 
## Std. Errors:
##          (Intercept) sesmiddle   seshigh      write
## academic    1.166437 0.4437319 0.5142215 0.02141092
## vocation    1.174251 0.4901237 0.6484508 0.02333135
## 
## Residual Deviance: 359.9635 
## AIC: 375.9635
stargazer(test1,type = "text")  # ¿y esto que? Interprtemos enseguida
## 
## ==============================================
##                       Dependent variable:     
##                   ----------------------------
##                      academic      vocation   
##                        (1)            (2)     
## ----------------------------------------------
## sesmiddle             0.533         0.825*    
##                      (0.444)        (0.490)   
##                                               
## seshigh              1.163**         0.180    
##                      (0.514)        (0.648)   
##                                               
## write                0.058***      -0.056**   
##                      (0.021)        (0.023)   
##                                               
## Constant             -2.852**       2.366**   
##                      (1.166)        (1.174)   
##                                               
## ----------------------------------------------
## Akaike Inf. Crit.    375.963        375.963   
## ==============================================
## Note:              *p<0.1; **p<0.05; ***p<0.01

Interpretacion

  1. Como tenemos tres categorias vamos a estimar dos vectores de parámetros \(\beta_1\) y \(\beta_2\)

  2. En general si tenemos \(J\) categoias estimaremos \(J-1\) vectores de parametros

  3. Para estimar estos parametros se hace de la funcion de verosimilitud apropiada

  4. En este ejemplo tipo de programa tiene tres categorias y la categoria base es programa general.

  5. En el logit multinomial la probabilidad de que el individuo \(i\) elija la categoria \(j\) es para \(j=1,...J-1\)

\[P_{ij}=\dfrac{e^{x_i\beta_j}}{1+\sum_{k=1}^{J-1}e^{x_i\beta_k}}\]

Y para la categoria base J:

\[P_{iJ}=\dfrac{1}{1+\sum_{k=1}^{J-1}e^{x_i\beta_k}}\]

  1. La descision por alguna de la categorias depende del vector de caracteristicas \(x_i\).

  2. En nuestro ejemplo

\[P_{i1}=\dfrac{e^{x_i\beta_1}}{1+e^{x_i\beta_1}+e^{x_i\beta_2}}\]

\[P_{i2}=\dfrac{e^{x_i\beta_2}}{1+e^{x_i\beta_1}+e^{x_i\beta_2}}\]

\[P_{i3}=\dfrac{1}{1+e^{x_i\beta_1}+e^{x_i\beta_2}}\]

Si podemos obtener los sig cocientes:

\[\dfrac{P_{i1}}{P_{i3}}=e^{x_i\beta_1}\]

\[\dfrac{P_{i2}}{P_{i3}}=e^{x_i\beta_2}\]

y entonces

\[\dfrac{P_{i1}}{P_{i2}}=e^{x_i(\beta_1-\beta_2)}\]

y si tomas logaritmos llegas a que:

\[\ln\left(\dfrac{P_{i1}}{P_{i3}}\right)=x_i\beta_1\]

\[\ln\left(\dfrac{P_{i2}}{P_{i3}}\right)=x_i\beta_2\]

\[\ln\left(\dfrac{P_{i1}}{P_{i2}}\right)=x_i(\beta_1-\beta_2)\]

En concreto para nuestro caso

\[\ln\left(\dfrac{P_{i1}}{P_{i3}}\right)=\beta_1^1+\beta_2^1*(ses=middle)+\beta_3^1*(ses=high)+\beta_4^1*write\]

\[\ln\left(\dfrac{P_{i2}}{P_{i3}}\right)=\beta_1^2+\beta_2^2*(ses=middle)+\beta_3^2*(ses=high)+\beta_4^2*write\]

\[\ln\left(\dfrac{P_{i1}}{P_{i2}}\right)=(\beta_1^1-\beta_1^2)+(\beta_2^1-\beta_2^2)*(ses=middle)+(\beta_3^1-\beta_3^2)*(ses=high)+(\beta_4^1-\beta_4^2)*write\]

O de manera áun más explícita, que para el \(i\)-esimo individuo:

\[\ln\left(\dfrac{P(prog=academic)}{P(prog=general)}\right)=\beta_1^1+\beta_2^1*(ses=middle)+\beta_3^1*(ses=high)+\beta_4^1*write\]

\[\ln\left(\dfrac{P(prog=vocation)}{P(prog=general)}\right)=\beta_1^2+\beta_2^2*(ses=middle)+\beta_3^2*(ses=high)+\beta_4^2*write\]

\[\ln\left(\dfrac{P(prog=academic)}{P(prog=vocation)}\right)=(\beta_1^1-\beta_1^2)+(\beta_2^1-\beta_2^2)*(ses=middle)+(\beta_3^1-\beta_3^2)*(ses=high)+(\beta_4^1-\beta_4^2)*write\]

  1. Los betas miden el efecto de cambios en \(x_i\) sobre el logaritmo de la razon de probabilidades(relative risk ratio) de elegir el programa j=1,2 respecto al programa general

\(\beta_4^1\) mide el efecto de cambios en el score de escritura sobre el log de la razon de probabilidades de elegir el programa academico respecto a el programa general (si aumenta el score en write es mas probable que un alumno elija el programa academico en vez de el programa general) \(\beta_4^2\) mide el efectode cambios en el score de escritura sobre el log de la razon de probabilidades de elegir el programa vocacional respecto a el programa general(si aumenta el score en write es menos probable que un alumno elija el programa vocacional en lugar de el programa general)

Si queremos ver el efecto de un cambio en el score de escritura sobre la decision de elegir programa academico sobre el programa vocacional el efecto es \((\beta_4^1-\beta_4^2)\)

  1. El analisis solo permite compara entre pares de alternativas
  2. El cambio de base no deberia alterar las conclusiones del modelo, solo tener cuidado con la interpretacion

##Efectos marginales en logit multinomial 11. Es el efecto de un cambio en la variable sobre la probabilidad absoluta de caer en alguna de las categorias. 12. A diferencia de las betas su interpretacion no es relativa sino absoluta: el impaacto sobre la posibilidad de estar en una categoria 13. no depende de la base

  1. Efecto marginal de un cambio en el score de escritura sobre la probabilidades de que el individuo \(i\) elija el programa academico o el programa vocacional \[\dfrac{\partial P_{ij}}{\partial x_i}=P_{ij}(\beta_j-\sum_{j}P_{ij}\beta_j)\]

\[\dfrac{\partial P_{i1}}{\partial write}=P_{i1}(\beta_1^w-P_{i1}\beta_1^w-P_{i2}\beta_2^w)\]

\[\dfrac{\partial P_{i2}}{\partial write}=P_{i2}(\beta_2^w-P_{i1}\beta_1^w-P_{i2}\beta_2^w)\] 15. El signo del efecto marginal puede ser distinto del signo del parametro beta

Probabilidades relativas

  1. Exp(B) es la razon de probabilidades asociada con cada predictor. Esperamos que los predictores que aumentan el logit muestren Exp (B) mayor que 1.0, aquellos predictores que no tienen un efecto en el logit mostraran un Exp (B) de 1.0 y los predictores que disminuyen el logit tendran valores Exp (B) menos de 1.0.
exp(coef(test1))
##          (Intercept) sesmiddle  seshigh     write
## academic  0.05773033  1.704533 3.198960 1.0596353
## vocation 10.65572511  2.281056 1.197478 0.9458464

Ante un cambio unitario en score de write,es un 5% mas probable que se elija un programa academico respecto aun programa general.

Por lo tanto, la probabilidad relativa de optar por un programa de tipo academico en lugar de un programa general es un 70% mas alta para los jovenes de nivel socieconomico medio que para los de nivel socieconomico bajo con el mismo puntaje en escritura.

Ojo: Solo la probabilidad relativa de elegir programa academico sobre programa general es mayor. No estamos diciendo que los jovenes de estatus medio tengan mas probabilidad de elegir programa academico.

La probabilidad relativa de elegir progrma tipo academico en lugar de un programa general para los alumnos de estatus socieconomico alto es mas del triple de la probabilidad relativa correspondiente para los chavos de nivel socieconomico bajo con el mismo nivel de score de escritura.

El efecto marginal promedio de un joven de status alto sobre el programa academico: -0.16. Esto significa que la probabilidad de elegir un programa academico es, en promedio, 16 puntos porcentuales mas baja para los de estatus alto que para los de estatus bajo con el mismo nivel escritura

pp <- fitted(test1)
head(pp)
##     general  academic  vocation
## 1 0.3382355 0.1482852 0.5134793
## 2 0.1806255 0.1202128 0.6991617
## 3 0.2367932 0.4186802 0.3445267
## 4 0.3508282 0.1726975 0.4764743
## 5 0.1689350 0.1001332 0.7309318
## 6 0.2377813 0.3533635 0.4088552
#summary(pp)
#str(pp)

###Asi cambias a otra categoria base

#Ajuste 2:Con otra categoria base
ml$prog2<-relevel(ml$prog, ref = "academic")
test2<-multinom(prog2 ~ ses + write, data = ml)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.982880
## final  value 179.981726 
## converged
#summary(test): salida para cada una de la ecuaciones
stargazer(test2,type = "text")  #resumen
## 
## ==============================================
##                       Dependent variable:     
##                   ----------------------------
##                      general       vocation   
##                        (1)            (2)     
## ----------------------------------------------
## sesmiddle             -0.533         0.291    
##                      (0.444)        (0.476)   
##                                               
## seshigh              -1.163**       -0.983*   
##                      (0.514)        (0.596)   
##                                               
## write               -0.058***      -0.114***  
##                      (0.021)        (0.022)   
##                                               
## Constant             2.852**       5.218***   
##                      (1.166)        (1.164)   
##                                               
## ----------------------------------------------
## Akaike Inf. Crit.    375.963        375.963   
## ==============================================
## Note:              *p<0.1; **p<0.05; ***p<0.01

Interpretacion

\[\ln\left(\dfrac{P(prog=general)}{P(prog=academic)}\right)=2.852-0.533*(ses=middle) -1.163*(ses=high)-0.058*write\]

\[\ln\left(\dfrac{P(prog=vocation)}{P(prog=academic)}\right)= 5.218+0.291*(ses=middle)-0.983*(ses=high)-0.114*write\]

\[\ln\left(\dfrac{P(prog=general)}{P(prog=vocation)}\right)=(2.852-5.218)+(-0.533 -0.291)*(ses=middle)+(-1.163-(-0.983))*(ses=high)+(-0.058-(-0.114))*write\]

Un incremento unitario en el score de escritura lleva a una decremento de 0.058 en el logaritmo de la razón de probabilidades de elegir el programa general respecto al programa academico.

Ante un incremento unitario en el score de escritura se tiene un decremento de 0.114 en el logaritmo de la razón de probabilidades de elegir el programa vocacional respecto al programa academico.

El logaritmo de la razón de probabilidades de elegir un programa general respecto a elegir un programa academico decrese en 1.163 si pasamos de un nivel socieconomico bajo a uno nivel socieconomico alto.

p valor de los parametros estimados

#valor Z
z <- summary(test2)$coefficients/summary(test2)$standard.errors
z
##          (Intercept)  sesmiddle   seshigh     write
## general     2.445214 -1.2018081 -2.261334 -2.705562
## vocation    4.484769  0.6116747 -1.649967 -5.112689
# test z a dos colas: p-valor
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p
##           (Intercept) sesmiddle    seshigh           write
## general  0.0144766100 0.2294379 0.02373856 0.0068189015731
## vocation 0.0000072993 0.5407530 0.09894976 0.0000003176045
## extraer los coeficientes del model y exponenciarlos
exp(coef(test2))
##          (Intercept) sesmiddle   seshigh     write
## general     17.32582 0.5866769 0.3126026 0.9437172
## vocation   184.61262 1.3382809 0.3743123 0.8926116
#Algunos valores predichos
head(pp <- fitted(test2))
##    academic   general  vocation
## 1 0.1482764 0.3382454 0.5134781
## 2 0.1202017 0.1806283 0.6991700
## 3 0.4186747 0.2368082 0.3445171
## 4 0.1726885 0.3508384 0.4764731
## 5 0.1001231 0.1689374 0.7309395
## 6 0.3533566 0.2377976 0.4088458
dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test2, newdata = dses, "probs")
##    academic   general  vocation
## 1 0.4396845 0.3581917 0.2021238
## 2 0.4777488 0.2283353 0.2939159
## 3 0.7009007 0.1784939 0.1206054
dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
                                                                                   3))
##almacena el predicho probabilidades para cada valor de ses and write
pp.write <- cbind(dwrite, predict(test2, newdata = dwrite, type = "probs", se = TRUE))

##calcular las probabilidades medias dentro de cada nivel de ses
by(pp.write[, 3:5], pp.write$ses, colMeans)
## pp.write$ses: high
##  academic   general  vocation 
## 0.6164315 0.1808037 0.2027648 
## ------------------------------------------------------------ 
## pp.write$ses: low
##  academic   general  vocation 
## 0.3972977 0.3278174 0.2748849 
## ------------------------------------------------------------ 
## pp.write$ses: middle
##  academic   general  vocation 
## 0.4256198 0.2010864 0.3732938
## usando melt y  ggplot2 
lpp  <-  melt (pp.write,  id.vars  =  c ( "ses" ,  "write" ),  value.name  =  "probabilidad" ) 
head (lpp)   # ver las primeras combinaciones
##   ses write variable probabilidad
## 1 low    30 academic   0.09843588
## 2 low    31 academic   0.10716868
## 3 low    32 academic   0.11650390
## 4 low    33 academic   0.12645834
## 5 low    34 academic   0.13704576
## 6 low    35 academic   0.14827643

Graficamos las probabilidades predichas a travess de valores write para cada nivel de ses Hacemos un plot por tipo de programa

ggplot(lpp, aes(x = write, y = probabilidad, colour = ses)) + geom_line(size=1)+ facet_grid(variable ~., scales = "free")

El modelo logit condicional

R ofrece varias alternativas para ajustar un modelo logit condicional

  1. Puedes usar la función mlogit() del paquete mlogit

  2. Puedes usar la función gmnl() del paquete gmnl

  3. Puedes usar la función mclogit() del paquete mclogit

  4. Puedes usar la función MCMCmnl() del paquete MCMCpack

  5. Puedes usar la función clogistic() del paquete Epi

Primera forma usando la función mlogit del pkg mlogit

data("Electricity", package = "mlogit")
Electr <- mlogit.data(Electricity, id.var = "id", choice = "choice",varying = 3:26, shape = "wide", sep = "")
head(Electr,10)
## ~~~~~~~
##  first 10 observations out of 17232 
## ~~~~~~~
##    choice id alt pf cl loc wk tod seas chid idx
## 1   FALSE  1   1  7  5   0  1   0    0    1 1:1
## 2   FALSE  1   2  9  1   1  0   0    0    1 1:2
## 3   FALSE  1   3  0  0   0  0   0    1    1 1:3
## 4    TRUE  1   4  0  5   0  1   1    0    1 1:4
## 5   FALSE  1   1  7  0   0  1   0    0    2 2:1
## 6   FALSE  1   2  9  5   0  1   0    0    2 2:2
## 7    TRUE  1   3  0  1   1  0   1    0    2 2:3
## 8   FALSE  1   4  0  5   0  0   0    1    2 2:4
## 9   FALSE  1   1  9  5   0  0   0    0    3 3:1
## 10  FALSE  1   2  7  1   0  1   0    0    3 3:2
## 
## ~~~ indexes ~~~~
##    chid id alt
## 1     1  1   1
## 2     1  1   2
## 3     1  1   3
## 4     1  1   4
## 5     2  1   1
## 6     2  1   2
## 7     2  1   3
## 8     2  1   4
## 9     3  1   1
## 10    3  1   2
## indexes:  1, 1, 2
clogit1 <- mlogit(choice ~ pf + cl + loc + wk + tod + seas | 0, data = Electr)
summary(clogit1)
## 
## Call:
## mlogit(formula = choice ~ pf + cl + loc + wk + tod + seas | 0, 
##     data = Electr, method = "nr")
## 
## Frequencies of alternatives:choice
##       1       2       3       4 
## 0.22702 0.26393 0.23816 0.27089 
## 
## nr method
## 4 iterations, 0h:0m:1s 
## g'(-H)^-1g = 9.34E-07 
## gradient close to zero 
## 
## Coefficients :
##        Estimate Std. Error z-value              Pr(>|z|)    
## pf   -0.6252278  0.0232223 -26.924 < 0.00000000000000022 ***
## cl   -0.1082991  0.0082442 -13.136 < 0.00000000000000022 ***
## loc   1.4422429  0.0505571  28.527 < 0.00000000000000022 ***
## wk    0.9955040  0.0447801  22.231 < 0.00000000000000022 ***
## tod  -5.4627587  0.1837125 -29.735 < 0.00000000000000022 ***
## seas -5.8400308  0.1866779 -31.284 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -4958.6

Segunda forma usando la función gmnl del pkg gmnl

clogit2<- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0, data = Electr)
summary(clogit2)
## 
## Model estimated on: Mon Sep 13 06:01:16 2021 
## 
## Call:
## gmnl(formula = choice ~ pf + cl + loc + wk + tod + seas | 0, 
##     data = Electr, method = "nr")
## 
## Frequencies of categories:
## 
##       1       2       3       4 
## 0.22702 0.26393 0.23816 0.27089 
## 
## The estimation took: 0h:0m:2s 
## 
## Coefficients:
##        Estimate Std. Error z-value              Pr(>|z|)    
## pf   -0.6252278  0.0232223 -26.924 < 0.00000000000000022 ***
## cl   -0.1082991  0.0082442 -13.136 < 0.00000000000000022 ***
## loc   1.4422429  0.0505571  28.527 < 0.00000000000000022 ***
## wk    0.9955040  0.0447801  22.231 < 0.00000000000000022 ***
## tod  -5.4627587  0.1837125 -29.735 < 0.00000000000000022 ***
## seas -5.8400308  0.1866779 -31.284 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -4958.6
## Number of observations: 4308
## Number of iterations: 4
## Exit of MLE: gradient close to zero (gradtol)

Tercera forma usando mclogit del pkg mclogit

data(Transport)
head(Transport,10)
##    transport suburb  distance cost working        prop.true resp
## 1        bus      1  8.418094  4.0     212 0.00130138501561    2
## 2        car      1  0.000000 10.0     212 0.98296149144755  207
## 3      train      1  6.956365  3.7     212 0.01573712353685    3
## 4        bus      2  4.254465  4.0     211 0.37320718128614   89
## 5        car      2  0.000000 10.0     211 0.54666211060243  104
## 6      train      2  5.480114  3.7     211 0.08013070811143   18
## 7        bus      3  6.545856  4.0     201 0.02148267975951    4
## 8        car      3  0.000000 10.0     201 0.97851051800407  197
## 9      train      3 12.117690  3.7     201 0.00000680223642    0
## 10       bus      4 15.864880  4.0     171 0.00000001865187    0
l.condicional1<-mclogit(cbind(resp,suburb)~distance+cost,data=Transport)
## 
## Iteration 1 - deviance = 39.74973 - criterion = 0.8590917
## Iteration 2 - deviance = 10.50328 - criterion = 2.758244
## Iteration 3 - deviance = 9.231325 - criterion = 0.1363107
## Iteration 4 - deviance = 9.227742 - criterion = 0.0003840654
## Iteration 5 - deviance = 9.227742 - criterion = 0.000000003446475
## converged
summary(l.condicional1)
## 
## Call:
## mclogit(formula = cbind(resp, suburb) ~ distance + cost, data = Transport)
## 
##          Estimate Std. Error z value            Pr(>|z|)    
## distance -1.43940    0.05318  -27.07 <0.0000000000000002 ***
## cost     -0.97753    0.03987  -24.52 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Null Deviance:     2734 
## Residual Deviance: 9.228 
## Number of Fisher Scoring iterations:  5 
## Number of observations:  1994
data(electors)
head(electors,10)
##                  party      class Freq time econ.left welfare  auth
## 1                green    working   22 0.04      0.40     0.0 -0.75
## 2                labor    working  463 0.04      1.00     1.5  0.00
## 3            communist    working    6 0.04      2.50     1.0  0.50
## 4         conservative    working    2 0.04     -1.00    -1.0  1.00
## 5              liberal    working    0 0.04     -1.25     0.0 -0.75
## 6  right.wing.populist    working    7 0.04      0.00     0.1  2.00
## 7                green new.middle  157 0.04      0.40     0.0 -0.75
## 8                labor new.middle  117 0.04      1.00     1.5  0.00
## 9            communist new.middle    0 0.04      2.50     1.0  0.50
## 10        conservative new.middle   31 0.04     -1.00    -1.0  1.00
l.condicional2<-mclogit(cbind(Freq,interaction(time,class))~econ.left/class+welfare/class+auth/class, random=~1|party.time,data=within(electors,party.time<-interaction(party,time)))
## Iteration 1 - deviance = 1054.511 - criterion = 0.1598497
## Iteration 2 - deviance = 923.1626 - criterion = 0.02666473
## Iteration 3 - deviance = 890.3113 - criterion = 0.006530011
## Iteration 4 - deviance = 883.0567 - criterion = 0.0005723444
## Iteration 5 - deviance = 881.4344 - criterion = 0.00001387329
## Iteration 6 - deviance = 881.2041 - criterion = 0.0000001394381
## Iteration 7 - deviance = 881.1809 - criterion = 0.00000000104388
## converged
summary(l.condicional2)
## 
## Call:
## mclogit(formula = cbind(Freq, interaction(time, class)) ~ econ.left/class + 
##     welfare/class + auth/class, data = within(electors, party.time <- interaction(party, 
##     time)), random = ~1 | party.time)
## 
## Coefficents:
##                           Estimate Std. Error z value             Pr(>|z|)    
## econ.left                 -0.12603    0.18576  -0.678                0.497    
## welfare                    2.01955    0.29136   6.932     0.00000000000416 ***
## auth                       0.11409    0.15908   0.717                0.473    
## econ.left:classnew.middle -1.81703    0.09859 -18.430 < 0.0000000000000002 ***
## econ.left:classold.middle -3.13862    0.15780 -19.890 < 0.0000000000000002 ***
## classnew.middle:welfare   -0.89578    0.06573 -13.628 < 0.0000000000000002 ***
## classold.middle:welfare   -1.47904    0.13451 -10.996 < 0.0000000000000002 ***
## classnew.middle:auth      -1.43391    0.04855 -29.535 < 0.0000000000000002 ***
## classold.middle:auth       1.44109    0.05883  24.494 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Co-)Variances:
## Grouping level: party.time 
##             Estimate   Std.Err.
## (Intercept) 3.066      0.9423  
## 
## Null Deviance:     80580 
## Residual Deviance: 881.2 
## Number of Fisher Scoring iterations:  7
## Number of observations
##   Groups by party.time: 150
##   Individual observations:  37500

Cuarta forma usando el pkg MCMCpack

Supongamos que queremos estudiar el efecto del precio en la decisión de un individuo sobre la elección de una de las tres marcas de refrescos:

  1. pepsi

  2. sevenup

  3. coke

En el modelo logit condicional, la probabilidad que el individuo \(i\) elija la marca \(j\) está dada por la siguiente ecuación:

\[\begin{equation} p_{ij}=\frac{exp(\beta_{1j}+\beta_{2}price_{ij})}{exp(\beta_{11}+\beta_{2}price_{i1})+exp(\beta_{12}+\beta_{2}price_{i2})+exp(\beta_{13}+\beta_{2}price_{i3})} \label{eq:CondLogitPij} \end{equation}\]

A diferencia del modelo logístico multinomial, el coeficiente de la variables independiente price es el mismo para todas la categorías, pero el valor de la variable independiente es diferente para cada individuo.

install_git("https://github.com/ccolonescu/PoEdata")
## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install Rtools 4.0 from https://cran.r-project.org/bin/windows/Rtools/.
## Skipping install of 'PoEdata' from a git2r remote, the SHA1 (c41bee8f) has not changed since last install.
##   Use `force = TRUE` to force installation
library("PoEdata")
data("cola")
head(cola,15)
##    id choice price feature display
## 1   1      0  1.79       0       0
## 2   1      0  1.79       0       0
## 3   1      1  1.79       0       0
## 4   2      0  1.79       0       0
## 5   2      0  1.79       0       0
## 6   2      1  0.89       1       1
## 7   3      0  1.41       0       0
## 8   3      0  0.84       0       1
## 9   3      1  0.89       1       0
## 10  4      0  1.79       0       0
## 11  4      0  1.79       0       0
## 12  4      1  1.33       1       0
## 13  5      0  1.79       0       0
## 14  5      0  1.79       0       0
## 15  5      1  1.79       0       0
N <- nrow(cola)
N3 <- N/3
price1 <- cola$price[seq(1,N,by=3)]
price2 <- cola$price[seq(2,N,by=3)]
price3 <- cola$price[seq(3,N,by=3)]

bchoice <- rep("1", N3)
for (j in 1:N3){
    if(cola$choice[3*j-1]==1) bchoice[j] <- "2"
    if(cola$choice[3*j]==1) bchoice[j] <- "3"
 }
cola.clogit <- MCMCmnl(bchoice ~
    choicevar(price1, "b2", "1")+
    choicevar(price2, "b2", "2")+
    choicevar(price3, "b2", "3"),
    baseline="3", mcmc.method="IndMH")
## Calculating MLEs and large sample var-cov matrix.
## This may take a moment...
## Inverting Hessian to get large sample var-cov matrix.
sclogit <- summary(cola.clogit)
tabMCMC <- as.data.frame(sclogit$statistics)[,1:2]
row.names(tabMCMC)<- c("b2","b11","b12")
kable(tabMCMC, digits=4, align="c",
caption="Estimaciones para el Logit Condicional para los datos 'cola' ")
Estimaciones para el Logit Condicional para los datos ‘cola’
Mean SD
b2 -2.2991 0.1382
b11 0.2839 0.0610
b12 0.1037 0.0621 
pPepsi <- 1
pSevenup <- 1.25
pCoke <- 1.10
b13 <- 0
b2  <- tabMCMC$Mean[1]
b11 <- tabMCMC$Mean[2]
b12 <- tabMCMC$Mean[3]

# La probabilidad de que el individuo i elija Pepsi:
PiPepsi <- exp(b11+b2*pPepsi)/
          (exp(b11+b2*pPepsi)+exp(b12+b2*pSevenup)+
                              exp(b13+b2*pCoke))
PiPepsi
## [1] 0.4834922
# La probabilidad de que el individuo i elija Sevenup:
PiSevenup <- exp(b12+b2*pSevenup)/
          (exp(b11+b2*pPepsi)+exp(b12+b2*pSevenup)+
                              exp(b13+b2*pCoke))
PiSevenup
## [1] 0.2272638
# La probabilidad de que el individuo i elija Coke:
PiCoke <- 1-PiPepsi-PiSevenup
PiCoke
## [1] 0.2892439

\[p_{i, \,pepsi}=0.483\] \[p_{i, \,sevenup}=0.227\] \[p_{i, \,coke}=0.289\]

Las tres probabilidades son diferentes para diferentes individuos porque diferentes individuos enfrentan diferentes precios; en un modelo más complejo se pueden incluir otros regresores, algunos de los cuales pueden reflejar características individuales.

##Quinta forma usando clogistic del pkg Epi

res.clogistic <- clogistic(choice ~ price, strata = id, data = cola)
res.clogistic
## 
## Call: 
## clogistic(formula = choice ~ price, strata = id, data = cola)
## 
## 
## 
## 
##        coef exp(coef) se(coef)     z p
## price -2.23     0.108    0.132 -16.8 0
## 
## Likelihood ratio test=333  on 1 df, p=0, n=4521

#Referencias: 1. https://stats.idre.ucla.edu/sas/dae/multinomiallogistic-regression/

  1. https://bookdown.org/ccolonescu/RPoE4/