I know this is really late, but you can do it using the pequod package.
library(pequod)
model.peq <- lmres(outcome ~ predictor*moderator1*moderator2, centered = c("predictor","moderator1","moderator2"), data=mydata)
S_slopes<-simpleSlope(model.peq, pred="predictor", mod1="moderator1", mod2="moderator2")
S_slopes
In the example, it seems like the error is really small so all the slopes are significant. If you change the data slightly, you get only one significant slope:
set.seed(123)
predictor <- rnorm(1000, 10, 5)
moderator1 <- rnorm(1000, 100, 25)
moderator2 <- rnorm(1000, 50, 20)
outcome <- predictor*moderator1*moderator2*rnorm(20, 0)/10000 # I changed this to zero from 30
mydata <- data.frame(predictor, moderator1, moderator2, outcome)
model.peq <- lmres(outcome ~ predictor*moderator1*moderator2, centered = c("predictor","moderator1","moderator2"), data=mydata)
S_slopes<-simpleSlope(model.peq, pred="predictor", mod1="moderator1", mod2="moderator2")
S_slopes
It's not clear whether you want estimates of height for each individual man and woman (more of a classification problem) or to characterize the distribution of heights of each sex. I will assume the latter. You also do not specify what additional information you are using in your model, so I will confine myself to addressing the case where you only have height data (and sex data, in the case of non-US citizens).
I recommend simply fitting a mixture of distributions to the height data from the US only, because the distributions of height in men and women are reasonably different. This would estimate the parameters of two distributions that when summed together best describe the variation in the data. The parameters of these distributions (mean and variance, since a Gaussian distribution should work fine) give you the information you are after. The R packages mixtools and mixdist let you do this; I'm sure there are many more as well.
This solution may seem odd, because it leaves out all the information you have from outside the US, where you have know the sex and height of each individual. But I think it is justified because:
1) We have a very strong prior expectation that men are on average taller than women. Wikipedia's List of average human height worldwide shows not even one country or region where women are taller than men. So the identity of the distribution with the greater mean height is not really in doubt.
2) Integrating more specific information from the non-US data will likely involve making the assumption that the covariance between sex and height is the same outside the US as inside. But this is not entirely true - the same Wikipedia list indicates that the ratio of male to female heights varies between approximately 1.04 and 1.13.
3) Your international data may be much more complicated to analyse because people in different countries have wide variation in height distributions as well. You may therefore need to consider modelling mixtures of mixtures of distributions. This may also be true in the US, but it is likely to be less of a problem than a dataset that includes the Dutch (mean height: 184 cms) and Indonesians (mean height: 158 cms). And those are country-level averages; subpopulations differ to an even degree.
Best Answer
Are the time points relevant (like time of day) or are they just multiple samples? If they are just multiple samples, you might just take the mean over time points, then run a single regression. You should not run 4 separate regressions (what would that tell you, anyway?). If the time points are relevant (I mean, they carry information, regardless of if they're relevant to you), enter it as another predictor, such that you have y ~ x1 + x2 + ... + x_time