Linearly Regressing Pooled Data

Question

I'm trying to regress some simple pooled data. My data has 60 observations and three columns: Weight, Height, and Sex (female=1, male=0).

If I regress thus, Weight ~ Height + Sex, my model is fairly satisfactory, but the residuals are not homoscedastic (green errors are male, blue female):

I tried regressing on the log of Weight and/or Height, but that didn't do much. What should I do to make the residuals homescedastic and/or make my model more accurate? Any help would be appreciated.

Edit

Doing a generalized regression model gives the following.

Generalized least squares fit by REML
  Model: Weight ~ h + s 
  Data: P149 
       AIC      BIC    logLik
  514.2221 524.4374 -252.1111

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | Sex 
 Parameter estimates:
        0         1 
1.0000000 0.6685307 

Coefficients:
                 Value Std.Error   t-value p-value
(Intercept)  27.197499  51.88129  0.524226  0.6022
h             1.852382   0.75634  2.449128  0.0174
s           -25.284478   5.53300 -4.569755  0.0000

 Correlation: 
  (Intr) h     
h -0.997       
s -0.524  0.466

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.6655243 -0.6879858 -0.1839396  0.5628971  3.9857544 

Residual standard error: 22.13369 
Degrees of freedom: 60 total; 57 residual

With this s. residual plot:

Could someone please explain how precisely this model is different from a standard multiple regression model? Thanks.

Answer 1

You should use a model framework that allows you to model the different variation in males and females. Package nlme provides a generalized least squares implementation (gls) that can be used in conjunction with the varIdent function to specify the appropriate variance structure. E.g., you could try something like

library(nlme)
gls(Weight ~ Height + Sex, weights=varIdent(form = ~ 1 | Sex), data=dat)