Patsy: Contrast Coding Systems for categorical variables¶
Note
This document is based heavily on this excellent resource from UCLA.
A categorical variable of K categories, or levels, usually enters a regression as a sequence of K-1 dummy variables. This amounts to a linear hypothesis on the level means. That is, each test statistic for these variables amounts to testing whether the mean for that level is statistically significantly different from the mean of the base category. This dummy coding is called Treatment coding in R parlance, and we will follow this convention. There are, however, different coding methods that amount to different sets of linear hypotheses.
In fact, the dummy coding is not technically a contrast coding. This is because the dummy variables add to one and are not functionally independent of the model’s intercept. On the other hand, a set of contrasts for a categorical variable with k levels is a set of k-1 functionally independent linear combinations of the factor level means that are also independent of the sum of the dummy variables. The dummy coding isn’t wrong per se. It captures all of the coefficients, but it complicates matters when the model assumes independence of the coefficients such as in ANOVA. Linear regression models do not assume independence of the coefficients and thus dummy coding is often the only coding that is taught in this context.
To have a look at the contrast matrices in Patsy, we will use data from UCLA ATS. First let’s load the data.
Example Data¶
In [1]: import pandas
In [2]: url = 'http://www.ats.ucla.edu/stat/data/hsb2.csv'
In [3]: hsb2 = pandas.read_table(url, delimiter=",")
---------------------------------------------------------------------------
ConnectionRefusedError Traceback (most recent call last)
/usr/lib/python3.5/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
1253 try:
-> 1254 h.request(req.get_method(), req.selector, req.data, headers)
1255 except OSError as err: # timeout error
/usr/lib/python3.5/http/client.py in request(self, method, url, body, headers)
1106 """Send a complete request to the server."""
-> 1107 self._send_request(method, url, body, headers)
1108
/usr/lib/python3.5/http/client.py in _send_request(self, method, url, body, headers)
1151 body = _encode(body, 'body')
-> 1152 self.endheaders(body)
1153
/usr/lib/python3.5/http/client.py in endheaders(self, message_body)
1102 raise CannotSendHeader()
-> 1103 self._send_output(message_body)
1104
/usr/lib/python3.5/http/client.py in _send_output(self, message_body)
933
--> 934 self.send(msg)
935 if message_body is not None:
/usr/lib/python3.5/http/client.py in send(self, data)
876 if self.auto_open:
--> 877 self.connect()
878 else:
/usr/lib/python3.5/http/client.py in connect(self)
848 self.sock = self._create_connection(
--> 849 (self.host,self.port), self.timeout, self.source_address)
850 self.sock.setsockopt(socket.IPPROTO_TCP, socket.TCP_NODELAY, 1)
/usr/lib/python3.5/socket.py in create_connection(address, timeout, source_address)
711 if err is not None:
--> 712 raise err
713 else:
/usr/lib/python3.5/socket.py in create_connection(address, timeout, source_address)
702 sock.bind(source_address)
--> 703 sock.connect(sa)
704 return sock
ConnectionRefusedError: [Errno 111] Connection refused
During handling of the above exception, another exception occurred:
URLError Traceback (most recent call last)
<ipython-input-3-626fbea21bbc> in <module>()
----> 1 hsb2 = pandas.read_table(url, delimiter=",")
/usr/lib/python3/dist-packages/pandas/io/parsers.py in parser_f(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, escapechar, comment, encoding, dialect, tupleize_cols, error_bad_lines, warn_bad_lines, skipfooter, skip_footer, doublequote, delim_whitespace, as_recarray, compact_ints, use_unsigned, low_memory, buffer_lines, memory_map, float_precision)
653 skip_blank_lines=skip_blank_lines)
654
--> 655 return _read(filepath_or_buffer, kwds)
656
657 parser_f.__name__ = name
/usr/lib/python3/dist-packages/pandas/io/parsers.py in _read(filepath_or_buffer, kwds)
390 compression = _infer_compression(filepath_or_buffer, compression)
391 filepath_or_buffer, _, compression = get_filepath_or_buffer(
--> 392 filepath_or_buffer, encoding, compression)
393 kwds['compression'] = compression
394
/usr/lib/python3/dist-packages/pandas/io/common.py in get_filepath_or_buffer(filepath_or_buffer, encoding, compression)
184 if _is_url(filepath_or_buffer):
185 url = str(filepath_or_buffer)
--> 186 req = _urlopen(url)
187 content_encoding = req.headers.get('Content-Encoding', None)
188 if content_encoding == 'gzip':
/usr/lib/python3.5/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
161 else:
162 opener = _opener
--> 163 return opener.open(url, data, timeout)
164
165 def install_opener(opener):
/usr/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
464 req = meth(req)
465
--> 466 response = self._open(req, data)
467
468 # post-process response
/usr/lib/python3.5/urllib/request.py in _open(self, req, data)
482 protocol = req.type
483 result = self._call_chain(self.handle_open, protocol, protocol +
--> 484 '_open', req)
485 if result:
486 return result
/usr/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
442 for handler in handlers:
443 func = getattr(handler, meth_name)
--> 444 result = func(*args)
445 if result is not None:
446 return result
/usr/lib/python3.5/urllib/request.py in http_open(self, req)
1280
1281 def http_open(self, req):
-> 1282 return self.do_open(http.client.HTTPConnection, req)
1283
1284 http_request = AbstractHTTPHandler.do_request_
/usr/lib/python3.5/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
1254 h.request(req.get_method(), req.selector, req.data, headers)
1255 except OSError as err: # timeout error
-> 1256 raise URLError(err)
1257 r = h.getresponse()
1258 except:
URLError: <urlopen error [Errno 111] Connection refused>
It will be instructive to look at the mean of the dependent variable, write, for each level of race ((1 = Hispanic, 2 = Asian, 3 = African American and 4 = Caucasian)).
Treatment (Dummy) Coding¶
Dummy coding is likely the most well known coding scheme. It compares each level of the categorical variable to a base reference level. The base reference level is the value of the intercept. It is the default contrast in Patsy for unordered categorical factors. The Treatment contrast matrix for race would be
In [4]: from patsy.contrasts import Treatment
In [5]: levels = [1,2,3,4]
In [6]: contrast = Treatment(reference=0).code_without_intercept(levels)
In [7]: print(contrast.matrix)
[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
Here we used reference=0, which implies that the first level, Hispanic, is the reference category against which the other level effects are measured. As mentioned above, the columns do not sum to zero and are thus not independent of the intercept. To be explicit, let’s look at how this would encode the race variable.
In [8]: contrast.matrix[hsb2.race-1, :][:20]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-8-7f2c2108d810> in <module>()
----> 1 contrast.matrix[hsb2.race-1, :][:20]
NameError: name 'hsb2' is not defined
This is a bit of a trick, as the race category conveniently maps to zero-based indices. If it does not, this conversion happens under the hood, so this won’t work in general but nonetheless is a useful exercise to fix ideas. The below illustrates the output using the three contrasts above
In [9]: from statsmodels.formula.api import ols
In [10]: mod = ols("write ~ C(race, Treatment)", data=hsb2)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-10-726467a33d67> in <module>()
----> 1 mod = ols("write ~ C(race, Treatment)", data=hsb2)
NameError: name 'hsb2' is not defined
In [11]: res = mod.fit()