Interaction effects in patsy with patsy.dmatrices giving duplicate columns for ":" as with "+" , or "*"
Asked Answered
G

2

5

I have a dataframe with columns, both of which I intend to treat as categorical variables.

the first column is country , which has values such as SGP, AUS, MYS etc. The second column is time of day, which has values in 24 hour format such as 00, 11, 14, 15 etc. event is a binary variable that has 1/0 flags. I understand that to categorize them , I need to use patsy before running the Logistic regression. This, I build using dmatrices.

Usecase : Consider only interaction effects of country & time_day (along with other attributes say "operating system")

f= 'event_int ~ time_day:country'
y,X = patsy.dmatrices(f, df, return_type='dataframe')
X.columns
Index([u'Intercept', u'country[T.HKG]', u'country[T.IDN]', u'country[T.IND]', u'country[T.MYS]', u'country[T.NZL]', u'country[T.PHL]', u'country[T.SGP]', u'time_day[T.02]:country[AUS]', u'time_day[T.03]:country[AUS]', u'time_day[T.04]:country[AUS]', u'time_day[T.05]:country[AUS]', u'time_day[T.06]:country[AUS]', u'time_day[T.07]:country[AUS]', u'time_day[T.08]:country[AUS]', u'time_day[T.09]:country[AUS]', u'time_day[T.10]:country[AUS]', u'time_day[T.11]:country[AUS]', u'time_day[T.12]:country[AUS]', u'time_day[T.NA]:country[AUS]', u'time_day[T.02]:country[HKG]', u'time_day[T.03]:country[HKG]', u'time_day[T.04]:country[HKG]', u'time_day[T.05]:country[HKG]', u'time_day[T.06]:country[HKG]', u'time_day[T.07]:country[HKG]', u'time_day[T.08]:country[HKG]', u'time_day[T.09]:country[HKG]', u'time_day[T.10]:country[HKG]', u'time_day[T.11]:country[HKG]', u'time_day[T.12]:country[HKG]', u'time_day[T.NA]:country[HKG]', u'time_day[T.02]:country[IDN]', u'time_day[T.03]:country[IDN]', u'time_day[T.04]:country[IDN]', u'time_day[T.05]:country[IDN]', u'time_day[T.06]:country[IDN]', u'time_day[T.07]:country[IDN]', u'time_day[T.08]:country[IDN]', u'time_day[T.09]:country[IDN]', u'time_day[T.10]:country[IDN]', u'time_day[T.11]:country[IDN]', u'time_day[T.12]:country[IDN]', u'time_day[T.NA]:country[IDN]', u'time_day[T.02]:country[IND]', u'time_day[T.03]:country[IND]', u'time_day[T.04]:country[IND]', u'time_day[T.05]:country[IND]', u'time_day[T.06]:country[IND]', u'time_day[T.07]:country[IND]', u'time_day[T.08]:country[IND]', u'time_day[T.09]:country[IND]', u'time_day[T.10]:country[IND]', u'time_day[T.11]:country[IND]', u'time_day[T.12]:country[IND]', u'time_day[T.NA]:country[IND]', u'time_day[T.02]:country[MYS]', u'time_day[T.03]:country[MYS]', u'time_day[T.04]:country[MYS]', u'time_day[T.05]:country[MYS]', u'time_day[T.06]:country[MYS]', u'time_day[T.07]:country[MYS]', u'time_day[T.08]:country[MYS]', u'time_day[T.09]:country[MYS]', u'time_day[T.10]:country[MYS]', u'time_day[T.11]:country[MYS]', u'time_day[T.12]:country[MYS]', u'time_day[T.NA]:country[MYS]', u'time_day[T.02]:country[NZL]', u'time_day[T.03]:country[NZL]', u'time_day[T.04]:country[NZL]', u'time_day[T.05]:country[NZL]', u'time_day[T.06]:country[NZL]', u'time_day[T.07]:country[NZL]', u'time_day[T.08]:country[NZL]', u'time_day[T.09]:country[NZL]', u'time_day[T.10]:country[NZL]', u'time_day[T.11]:country[NZL]', u'time_day[T.12]:country[NZL]', u'time_day[T.NA]:country[NZL]', u'time_day[T.02]:country[PHL]', u'time_day[T.03]:country[PHL]', u'time_day[T.04]:country[PHL]', u'time_day[T.05]:country[PHL]', u'time_day[T.06]:country[PHL]', u'time_day[T.07]:country[PHL]', u'time_day[T.08]:country[PHL]', u'time_day[T.09]:country[PHL]', u'time_day[T.10]:country[PHL]', u'time_day[T.11]:country[PHL]', u'time_day[T.12]:country[PHL]', u'time_day[T.NA]:country[PHL]', u'time_day[T.02]:country[SGP]', u'time_day[T.03]:country[SGP]', u'time_day[T.04]:country[SGP]', u'time_day[T.05]:country[SGP]', u'time_day[T.06]:country[SGP]', u'time_day[T.07]:country[SGP]', u'time_day[T.08]:country[SGP]', u'time_day[T.09]:country[SGP]', ...], dtype='object')

I hoped to see only the column names with BOTH country & time_day, but this is not the case. I could manually take a subset by specifying X = X.ix[:,range(7,len(X.columns))] , but this would mean HARDCODING this for each dataset.

My understanding was A*B differs from A:B in the sense that it does not list out A+B Interesting thing though is that I do not see A, ie Categorical values of time_day alone, in the output above.

Also, when I do the following, to explicitly exclude "country" alone from the "X" dataframe, it doesn’t work, and I get the same output as above.

f='event_int ~ time_day:country-country'
y,X = patsy.dmatrices(f, df, return_type='dataframe')
X.columns
Index([u'Intercept', u'country[T.HKG]', u'country[T.IDN]', u'country[T.IND]', u'country[T.MYS]', u'country[T.NZL]', u'country[T.PHL]', u'country[T.SGP]', u'time_day[T.02]:country[AUS]', u'time_day[T.03]:country[AUS]', u'time_day[T.04]:country[AUS]', u'time_day[T.05]:country[AUS]', u'time_day[T.06]:country[AUS]', u'time_day[T.07]:country[AUS]', u'time_day[T.08]:country[AUS]', u'time_day[T.09]:country[AUS]', u'time_day[T.10]:country[AUS]', u'time_day[T.11]:country[AUS]', u'time_day[T.12]:country[AUS]', u'time_day[T.NA]:country[AUS]', u'time_day[T.02]:country[HKG]', u'time_day[T.03]:country[HKG]', u'time_day[T.04]:country[HKG]', u'time_day[T.05]:country[HKG]', u'time_day[T.06]:country[HKG]', u'time_day[T.07]:country[HKG]', u'time_day[T.08]:country[HKG]', u'time_day[T.09]:country[HKG]', u'time_day[T.10]:country[HKG]', u'time_day[T.11]:country[HKG]', u'time_day[T.12]:country[HKG]', u'time_day[T.NA]:country[HKG]', u'time_day[T.02]:country[IDN]', u'time_day[T.03]:country[IDN]', u'time_day[T.04]:country[IDN]', u'time_day[T.05]:country[IDN]', u'time_day[T.06]:country[IDN]', u'time_day[T.07]:country[IDN]', u'time_day[T.08]:country[IDN]', u'time_day[T.09]:country[IDN]', u'time_day[T.10]:country[IDN]', u'time_day[T.11]:country[IDN]', u'time_day[T.12]:country[IDN]', u'time_day[T.NA]:country[IDN]', u'time_day[T.02]:country[IND]', u'time_day[T.03]:country[IND]', u'time_day[T.04]:country[IND]', u'time_day[T.05]:country[IND]', u'time_day[T.06]:country[IND]', u'time_day[T.07]:country[IND]', u'time_day[T.08]:country[IND]', u'time_day[T.09]:country[IND]', u'time_day[T.10]:country[IND]', u'time_day[T.11]:country[IND]', u'time_day[T.12]:country[IND]', u'time_day[T.NA]:country[IND]', u'time_day[T.02]:country[MYS]', u'time_day[T.03]:country[MYS]', u'time_day[T.04]:country[MYS]', u'time_day[T.05]:country[MYS]', u'time_day[T.06]:country[MYS]', u'time_day[T.07]:country[MYS]', u'time_day[T.08]:country[MYS]', u'time_day[T.09]:country[MYS]', u'time_day[T.10]:country[MYS]', u'time_day[T.11]:country[MYS]', u'time_day[T.12]:country[MYS]', u'time_day[T.NA]:country[MYS]', u'time_day[T.02]:country[NZL]', u'time_day[T.03]:country[NZL]', u'time_day[T.04]:country[NZL]', u'time_day[T.05]:country[NZL]', u'time_day[T.06]:country[NZL]', u'time_day[T.07]:country[NZL]', u'time_day[T.08]:country[NZL]', u'time_day[T.09]:country[NZL]', u'time_day[T.10]:country[NZL]', u'time_day[T.11]:country[NZL]', u'time_day[T.12]:country[NZL]', u'time_day[T.NA]:country[NZL]', u'time_day[T.02]:country[PHL]', u'time_day[T.03]:country[PHL]', u'time_day[T.04]:country[PHL]', u'time_day[T.05]:country[PHL]', u'time_day[T.06]:country[PHL]', u'time_day[T.07]:country[PHL]', u'time_day[T.08]:country[PHL]', u'time_day[T.09]:country[PHL]', u'time_day[T.10]:country[PHL]', u'time_day[T.11]:country[PHL]', u'time_day[T.12]:country[PHL]', u'time_day[T.NA]:country[PHL]', u'time_day[T.02]:country[SGP]', u'time_day[T.03]:country[SGP]', u'time_day[T.04]:country[SGP]', u'time_day[T.05]:country[SGP]', u'time_day[T.06]:country[SGP]', u'time_day[T.07]:country[SGP]', u'time_day[T.08]:country[SGP]', u'time_day[T.09]:country[SGP]', ...], dtype='object')

This makes me feel that ":" is a reduced form of "*" in that it misses just ONE categorical var. I think it is not able to understand that BOTH are categorical vars ?

f='event_int ~ time_day*country'
y,X = patsy.dmatrices(f, df, return_type='dataframe')
X.columns
Index([u'Intercept', u'time_day[T.02]', u'time_day[T.03]', u'time_day[T.04]', u'time_day[T.05]', u'time_day[T.06]', u'time_day[T.07]', u'time_day[T.08]', u'time_day[T.09]', u'time_day[T.10]', u'time_day[T.11]', u'time_day[T.12]', u'time_day[T.NA]', u'country[T.HKG]', u'country[T.IDN]', u'country[T.IND]', u'country[T.MYS]', u'country[T.NZL]', u'country[T.PHL]', u'country[T.SGP]', u'time_day[T.02]:country[T.HKG]', u'time_day[T.03]:country[T.HKG]', u'time_day[T.04]:country[T.HKG]', u'time_day[T.05]:country[T.HKG]', u'time_day[T.06]:country[T.HKG]', u'time_day[T.07]:country[T.HKG]', u'time_day[T.08]:country[T.HKG]', u'time_day[T.09]:country[T.HKG]', u'time_day[T.10]:country[T.HKG]', u'time_day[T.11]:country[T.HKG]', u'time_day[T.12]:country[T.HKG]', u'time_day[T.NA]:country[T.HKG]', u'time_day[T.02]:country[T.IDN]', u'time_day[T.03]:country[T.IDN]', u'time_day[T.04]:country[T.IDN]', u'time_day[T.05]:country[T.IDN]', u'time_day[T.06]:country[T.IDN]', u'time_day[T.07]:country[T.IDN]', u'time_day[T.08]:country[T.IDN]', u'time_day[T.09]:country[T.IDN]', u'time_day[T.10]:country[T.IDN]', u'time_day[T.11]:country[T.IDN]', u'time_day[T.12]:country[T.IDN]', u'time_day[T.NA]:country[T.IDN]', u'time_day[T.02]:country[T.IND]', u'time_day[T.03]:country[T.IND]', u'time_day[T.04]:country[T.IND]', u'time_day[T.05]:country[T.IND]', u'time_day[T.06]:country[T.IND]', u'time_day[T.07]:country[T.IND]', u'time_day[T.08]:country[T.IND]', u'time_day[T.09]:country[T.IND]', u'time_day[T.10]:country[T.IND]', u'time_day[T.11]:country[T.IND]', u'time_day[T.12]:country[T.IND]', u'time_day[T.NA]:country[T.IND]', u'time_day[T.02]:country[T.MYS]', u'time_day[T.03]:country[T.MYS]', u'time_day[T.04]:country[T.MYS]', u'time_day[T.05]:country[T.MYS]', u'time_day[T.06]:country[T.MYS]', u'time_day[T.07]:country[T.MYS]', u'time_day[T.08]:country[T.MYS]', u'time_day[T.09]:country[T.MYS]', u'time_day[T.10]:country[T.MYS]', u'time_day[T.11]:country[T.MYS]', u'time_day[T.12]:country[T.MYS]', u'time_day[T.NA]:country[T.MYS]', u'time_day[T.02]:country[T.NZL]', u'time_day[T.03]:country[T.NZL]', u'time_day[T.04]:country[T.NZL]', u'time_day[T.05]:country[T.NZL]', u'time_day[T.06]:country[T.NZL]', u'time_day[T.07]:country[T.NZL]', u'time_day[T.08]:country[T.NZL]', u'time_day[T.09]:country[T.NZL]', u'time_day[T.10]:country[T.NZL]', u'time_day[T.11]:country[T.NZL]', u'time_day[T.12]:country[T.NZL]', u'time_day[T.NA]:country[T.NZL]', u'time_day[T.02]:country[T.PHL]', u'time_day[T.03]:country[T.PHL]', u'time_day[T.04]:country[T.PHL]', u'time_day[T.05]:country[T.PHL]', u'time_day[T.06]:country[T.PHL]', u'time_day[T.07]:country[T.PHL]', u'time_day[T.08]:country[T.PHL]', u'time_day[T.09]:country[T.PHL]', u'time_day[T.10]:country[T.PHL]', u'time_day[T.11]:country[T.PHL]', u'time_day[T.12]:country[T.PHL]', u'time_day[T.NA]:country[T.PHL]', u'time_day[T.02]:country[T.SGP]', u'time_day[T.03]:country[T.SGP]', u'time_day[T.04]:country[T.SGP]', u'time_day[T.05]:country[T.SGP]', u'time_day[T.06]:country[T.SGP]', u'time_day[T.07]:country[T.SGP]', u'time_day[T.08]:country[T.SGP]', u'time_day[T.09]:country[T.SGP]', ...], dtype='object')

And if I were to explicitly declare them as "categorical" vars, I get this -:

f='event_int ~ C(time_day):C(country)'
y,X = patsy.dmatrices(f, df, return_type='dataframe')
X.columns
Index([u'Intercept', u'C(country)[T.HKG]', u'C(country)[T.IDN]', u'C(country)[T.IND]', u'C(country)[T.MYS]', u'C(country)[T.NZL]', u'C(country)[T.PHL]', u'C(country)[T.SGP]', u'C(time_day)[T.02]:C(country)[AUS]', u'C(time_day)[T.03]:C(country)[AUS]', u'C(time_day)[T.04]:C(country)[AUS]', u'C(time_day)[T.05]:C(country)[AUS]', u'C(time_day)[T.06]:C(country)[AUS]', u'C(time_day)[T.07]:C(country)[AUS]', u'C(time_day)[T.08]:C(country)[AUS]', u'C(time_day)[T.09]:C(country)[AUS]', u'C(time_day)[T.10]:C(country)[AUS]', u'C(time_day)[T.11]:C(country)[AUS]', u'C(time_day)[T.12]:C(country)[AUS]', u'C(time_day)[T.NA]:C(country)[AUS]', u'C(time_day)[T.02]:C(country)[HKG]', u'C(time_day)[T.03]:C(country)[HKG]', u'C(time_day)[T.04]:C(country)[HKG]', u'C(time_day)[T.05]:C(country)[HKG]', u'C(time_day)[T.06]:C(country)[HKG]', u'C(time_day)[T.07]:C(country)[HKG]', u'C(time_day)[T.08]:C(country)[HKG]', u'C(time_day)[T.09]:C(country)[HKG]', u'C(time_day)[T.10]:C(country)[HKG]', u'C(time_day)[T.11]:C(country)[HKG]', u'C(time_day)[T.12]:C(country)[HKG]', u'C(time_day)[T.NA]:C(country)[HKG]', u'C(time_day)[T.02]:C(country)[IDN]', u'C(time_day)[T.03]:C(country)[IDN]', u'C(time_day)[T.04]:C(country)[IDN]', u'C(time_day)[T.05]:C(country)[IDN]', u'C(time_day)[T.06]:C(country)[IDN]', u'C(time_day)[T.07]:C(country)[IDN]', u'C(time_day)[T.08]:C(country)[IDN]', u'C(time_day)[T.09]:C(country)[IDN]', u'C(time_day)[T.10]:C(country)[IDN]', u'C(time_day)[T.11]:C(country)[IDN]', u'C(time_day)[T.12]:C(country)[IDN]', u'C(time_day)[T.NA]:C(country)[IDN]', u'C(time_day)[T.02]:C(country)[IND]', u'C(time_day)[T.03]:C(country)[IND]', u'C(time_day)[T.04]:C(country)[IND]', u'C(time_day)[T.05]:C(country)[IND]', u'C(time_day)[T.06]:C(country)[IND]', u'C(time_day)[T.07]:C(country)[IND]', u'C(time_day)[T.08]:C(country)[IND]', u'C(time_day)[T.09]:C(country)[IND]', u'C(time_day)[T.10]:C(country)[IND]', u'C(time_day)[T.11]:C(country)[IND]', u'C(time_day)[T.12]:C(country)[IND]', u'C(time_day)[T.NA]:C(country)[IND]', u'C(time_day)[T.02]:C(country)[MYS]', u'C(time_day)[T.03]:C(country)[MYS]', u'C(time_day)[T.04]:C(country)[MYS]', u'C(time_day)[T.05]:C(country)[MYS]', u'C(time_day)[T.06]:C(country)[MYS]', u'C(time_day)[T.07]:C(country)[MYS]', u'C(time_day)[T.08]:C(country)[MYS]', u'C(time_day)[T.09]:C(country)[MYS]', u'C(time_day)[T.10]:C(country)[MYS]', u'C(time_day)[T.11]:C(country)[MYS]', u'C(time_day)[T.12]:C(country)[MYS]', u'C(time_day)[T.NA]:C(country)[MYS]', u'C(time_day)[T.02]:C(country)[NZL]', u'C(time_day)[T.03]:C(country)[NZL]', u'C(time_day)[T.04]:C(country)[NZL]', u'C(time_day)[T.05]:C(country)[NZL]', u'C(time_day)[T.06]:C(country)[NZL]', u'C(time_day)[T.07]:C(country)[NZL]', u'C(time_day)[T.08]:C(country)[NZL]', u'C(time_day)[T.09]:C(country)[NZL]', u'C(time_day)[T.10]:C(country)[NZL]', u'C(time_day)[T.11]:C(country)[NZL]', u'C(time_day)[T.12]:C(country)[NZL]', u'C(time_day)[T.NA]:C(country)[NZL]', u'C(time_day)[T.02]:C(country)[PHL]', u'C(time_day)[T.03]:C(country)[PHL]', u'C(time_day)[T.04]:C(country)[PHL]', u'C(time_day)[T.05]:C(country)[PHL]', u'C(time_day)[T.06]:C(country)[PHL]', u'C(time_day)[T.07]:C(country)[PHL]', u'C(time_day)[T.08]:C(country)[PHL]', u'C(time_day)[T.09]:C(country)[PHL]', u'C(time_day)[T.10]:C(country)[PHL]', u'C(time_day)[T.11]:C(country)[PHL]', u'C(time_day)[T.12]:C(country)[PHL]', u'C(time_day)[T.NA]:C(country)[PHL]', u'C(time_day)[T.02]:C(country)[SGP]', u'C(time_day)[T.03]:C(country)[SGP]', u'C(time_day)[T.04]:C(country)[SGP]', u'C(time_day)[T.05]:C(country)[SGP]', u'C(time_day)[T.06]:C(country)[SGP]', u'C(time_day)[T.07]:C(country)[SGP]', u'C(time_day)[T.08]:C(country)[SGP]', u'C(time_day)[T.09]:C(country)[SGP]', ...], dtype='object')

Questions :

1. How do I include ONLY interaction effects & nothing else for such variables ?

2. Why is the exclusion of country with -country not working in the second case?

Related :Statsmodels formula API (patsy): How to exclude a subset of interaction components?

Edited to sort-of troubleshoot yourself based on @Nathaniel J. Smith's answer below -:

f2='event_int ~ country:time_day'
y2,X2 = patsy.dmatrices(f2, df, return_type='dataframe')
X2.design_info.term_names
['Intercept', 'country:time_day']

f1='event_int ~ country:time_day-1'
y1,X1 = patsy.dmatrices(f1, df, return_type='dataframe')
X1.design_info.term_names
['country:time_day']
Galvanism answered 15/5, 2014 at 8:2 Comment(0)
C
10

Short answer: try event_int ~ -1 + time_day:country

Long answer:

The first thing to understand is that there are two different phases to how patsy decides to build a design matrix. First, it determines which terms to include. Terms are things like a, or a:b. (The a and b in a:b are called factors; the term a contains a single factor which is also spelled a.) Figuring out which terms exist involves expanding and simplifying the formula you give it, until you have an expression that only uses + and :. a*b expands into a + b + a:b, etc. Subtraction is an operation that happens at this stage: a + b - a simplifies to just plain b. So a*b - a expands to a + b + a:b - a which simplifies to b + a:b, but a:b - a is the same as a:b, because there's no a to subtract, so the - a is just ignored. This is why writing time_day:country - country is the same as writing time_day:country.

Then in the second phase, once patsy has decided which terms to include, it has to decide how to code those terms. It's this phase where you're running into trouble.

The general rule is that patsy goes through each term that has categorical factors in it, and figures out a set of columns it can use that will make the model flexible enough to include the specified interaction, but without being redundant with any terms that have already been added.

In this case, your trouble is being caused by the intercept term that patsy adds by default: event_int ~ time_day:country is interpreted like event_int ~ 1 + time_day:country. This tells patsy that you want to have one column represents the intercept term alone, and then a second group of columns which cover the interaction -- but that don't overlap with the intercept. The obvious approach of dummy-coding both time_day and country would be redundant (collinear) with the intercept, so patsy instead finds a somewhat complicated scheme that doesn't have this property. If you remove the intercept, then you tell patsy that it can go ahead and use the simple scheme, so it does.

The details of how patsy chooses a coding scheme are explained here: http://patsy.readthedocs.org/en/latest/formulas.html#redundancy-and-categorical-factors

The first part of that manual section has a bit too much math perhaps, but if you scroll down there are some hopefully-nice diagrams that may make it clearer what's going on (and provide some context for the math). If you search for y ~ 1 + a:b you'll see the diagram that specifically shows the situation you're ending up in when you type event_int ~ time_day:country. And if you search for y ~ 1 + a + b + a:b you'll see a picture of what's happening in the event_int ~ time_day*country case.

In addition to looking at X.columns, it's useful to look at X.design_info.term_names and X.design_info.term_slices, which show you which "terms" patsy thinks exist and which columns they correspond to. (a and a:b are terms; each one generates multiple columns.) The thick outline in the y ~ 1 + a:b diagram is intended to indicate that in this case, the single term a:b generates two sets of columns: one set of columns coding b with treatment coding, and a second set of columns coding the pairwise products of dummy-coded b and treatment-coded a.

Finally, two tips for interpreting the output you're getting: (1) you can be sure that patsy is in fact treating the factors as categorical, because the column names look like varname[something involving the var's value]. Numerical factors either look like varname or (in the rare case where you pass a 2d matrix as a predictor) varname[column index]. (2) Pay attention to the difference between country[T.HKG] and country[HKG] -- the former indicates that patsy is using reduced-rank "treatment" coding to avoid redundancy, while the latter indicates simple dummy-coding. Of course, it turns out that as far as individual columns go these are identical, but conceptually the difference is very important -- the T. mode means that it's dropped one of the columns (notice the absence of country[T.AUS]), so subsetting the columns like you considered doing would not have worked well!

Hope this helps!

Cannula answered 15/5, 2014 at 14:1 Comment(2)
Additional thanks for explaining why time_day:country - country is the same as writing time_day:country.If I understand you correct, I should AIM for getting columns with "T" , like country[T.HKG], so that they can directly feed into my model post the "y","X" dataframes are returned.Galvanism
Well, whether you want dummy coding or treatment coding (or some other kind of coding) depends on what you're doing. In general the idea is that you can directly take the y and X matrices returned by dmatrices and use them directly for linear regression (or related techniques). And 1 + a:b, -1 + a:b, 1 + a*b will all give you the same linear model in terms of predictions, $R^2$, etc. But the $\beta$ values will have different interpretations. I suggest consulting some textbooks or asking a question on Cross Validated if you need more details on this part.Cannula
W
3

It looks like we need to remove the constant to avoid the reference coding

    >>> f0 = 'event_int ~ C(time_day):C(country) - 1'
    >>> y,X0 = patsy.dmatrices(f0, df, return_type='dataframe')
    >>> print len(X0.columns)
    12
    >>> print X0.columns
    Index([u'C(time_day)[1]:C(country)[a]', u'C(time_day)[2]:C(country)[a]', 
           u'C(time_day)[3]:C(country)[a]', u'C(time_day)[4]:C(country)[a]', 
           u'C(time_day)[1]:C(country)[b]', u'C(time_day)[2]:C(country)[b]', 
           u'C(time_day)[3]:C(country)[b]', u'C(time_day)[4]:C(country)[b]', 
           u'C(time_day)[1]:C(country)[c]', u'C(time_day)[2]:C(country)[c]', 
           u'C(time_day)[3]:C(country)[c]', u'C(time_day)[4]:C(country)[c]'],
           dtype=object)

this is based on a simple example with 3 countries and 4 time periods:

import numpy as np
import pandas as pd
import patsy

time = np.tile([1, 2, 3, 4], 3)
country = np.repeat(['a', 'b', 'c'], 4)
event_int = np.random.randint(0, 2, size=len(time))

df = pd.DataFrame({'event_int':event_int, 'time_day':time, 'country':country})
Walling answered 15/5, 2014 at 13:36 Comment(3)
Both the Answers are correct, but I could accept only one. The answer above also hints at the other things I was originally missing, so accepted that one. Thank you.Galvanism
Nathaniel is the author of patsy and knows best. His answer doesn't show up as accepted yet.Walling
I did accept it - I think it was because of the Stackoverflow maintenance on Friday , since I first accepted yours & then re-read his answer - and accepted his.Galvanism

© 2022 - 2024 — McMap. All rights reserved.