how do I have a term for the co-occurrence of two variables in logistic regression? lets say I code variables in data:
sick: 0 or 1 for someone sick or not.
cough: 0 or 1 for someone having a cough or not.
fever: 0 or 1 for someone wearing a fever or or not.
cough_fever: 0 or 1 for someone having both a cough and fever or not
lets say I wanted to determine if having a cough, fever, or having both a cough and fever is most significantly associated with being sick. Would I set up my regression in R like this:
glm(sick ~ cough + fever + cough_fever, data = data, family='binomial')
or like this:
glm(sick ~ cough + fever + cough*fever, data = data, family='binomial')

An "interaction" term in a regression model is just another name for "product" regardless of whether you are doing ordinary least-squares, logistic, or other types of regression. If 0/1 means FALSE/TRUE (in that order) in your coding for each of cough and fever, then their product is 1 only when both are TRUE. In that case, either way will work: in principle.
There are three practical considerations in applying that principle, however, that argue against trying to construct the interaction term as a separate predictor yourself.
The first is terminological. The "*" symbol generally used for multiplication in an expression has a special meaning in an R formula. In an R formula, as Glen_b says in comments, cough*fever expands to each of the individual terms plus their product. With the ":" symbol representing just the product (as R interprets it within a formula), cough*fever in a formula expands to cough + fever + cough:fever. If R processes the last formula you showed, it expands cough*fever that way and removes the extra individual cough and fever terms before proceeding.
Second is the coding of categorical predictors. If you coded your binary predictors as 1/2 instead of 0/1 and tried to construct your own product/interaction term, you wouldn't get the same regression coefficients as if R had done the expansion of cough*fever for you. That's the case in this related question. With an unordered categorical predictor, R chooses one level to be the "reference" category and sets its internal coding to 0. A related source of confusion: the choice of 0-coded reference category that R makes can be different from what other software makes (e.g., choosing the first versus the last among category names in alphabetical order). Be careful in comparing results obtained with different statistical software systems.
Third, as Noah emphasizes in another answer, letting the software do the expansion can help with downstream processing of the regression results. For example, the rms package in R has some very useful functions for evaluating the overall significance of interactions and nonlinear terms in regressions. Those functions, however, need to know that there were interactions in the original model. If you code an interaction as a separate predictor yourself, you lose that downstream functionality.
