Warning message

The matter of status Red when analyzing threshold traits has been discussed endlessly on these forums. Try searching the forums for more information.

• Did you mean there's something wrong with the script?

Yes, there are several things wrong with the script. The fact that both the trait of interest and the putative moderator are binary variables requires considerations that aren't relevant to the continuous-only case.

The script correctly casts variables "smk_cur1" and "smk_cur2" to ordered factors (through `mxFactor()`). In R, variables of type "factor" are a compound data-type, consisting of a vector of integers (one unique integer value per level), as well as a vector of character strings (one unique string per level). The problem is that R begins enumerating levels of the factor at 1, not 0. So, if you use "smk_cur1" and "smk_cur2" as definition variables in MxAlgebras, they will be involved in calculations as 1s and 2s, not 0s and 1s as you might expect. Since the 0-and-1 numeric coding greatly simplifies interpretation, I suggest adding some lines like the following to the script, prior to the calls to `mxFactor()`:

mzDataF$smk_cur1_num <- as.numeric(mzDataF$smk_cur1)
mzDataF$smk_cur2_num <- as.numeric(mzDataF$smk_cur2)

Do likewise for the DZ dataset. Then, use those new variables with "_num" appended to their names as definition variables, e.g.

mod_tw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.smk_cur1_num', name="Mod1")

. Once you've created the "_num" variables, you can allow for incomplete twin pairs if you don't filter them out earlier in the script (as is currently done in lines 27 & 28). You can do something like the following,

mzDataF$coronary1[is.na(mzDataF$smk_cur1_num)] <- NA
mzDataF$smk_cur1_num[is.na(mzDataF$smk_cur1_num)] <- -999.0

, and likewise for MZ twin #2 and for the DZ twins.

To treat current smoking status as a moderator is equivalent to regressing "coronary" onto interactions between smoking status and the common and unique *A*, *C*, and *E* variables. As a rule, you also need to model the main effect of all variables appearing in interaction terms. In this case, that means you need to adjust for current smoking in the same way you currently adjust for age and sex.

The biggest problem with the script is that the model is unidentified, because the scale of the liabilities underlying smoking status and coronary aren't identified. As explained in Medland et al. (2009), fixing the liability-scale variance to a constant at some arbitrary value of the moderator is sufficient for model identification in the binary case. Replace lines 150 to 167 in the script with this:

myVarA <- mxAlgebra(name = "A0", expression = a %*% t(a))
myVarC <- mxAlgebra(name = "C0", expression = c %*% t(c))
myVarE <- mxAlgebra(name = "E0", expression = e %*% t(e))
myVar <- mxAlgebra( A0+C0+E0, name="V0" )
# Constraint on variance of Binary variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var_constraint <- mxConstraint( expression=diag2vec(V0)==Unv1, name="Var1" )

The MxConstraint will fix the liability-scales variances of both binary variables to 1.0 when the moderator is zero. Make sure `myVarA`, `myVarC`, `myVarE`, `myVar`, `matUnv`, and `var_constraint` all go into either your MZ or DZ model.

One minor thing: the script currently fixes the regression intercept at 1.5. There's nothing wrong with that, but most people fix it to zero to make interpretation easier.