[R-sig-ME] Linear Mixed Model with Twin Data; Heteroscedasticity?‏

[Frank Savage]

Very sorry!! -- first time poster and it looks like you have to follow an odd link to get to my question. I've relayed it here below for convenience.

(Let me preface this by stating that I am somewhat of an R newb and also a LMM newb.)

I have a very, very simple research question. I have a continuous outcome variable ("HAPPY"). I have two factors that divide my sample into 4 cells (Factor1, Factor2), and I want to assess whether there are any main effects or an interactions. A simple 2 x 2 ANOVA should do the trick.

There are two issues:

1) The individuals in my sample are family members. Every single individual has a fraternal twin or an identical twin in the dataset. I know that this throws off all the independence assumptions, so I was instructed to use linear mixed models and include the family identification variable (which identifies an individual as being part of a twin pair) ("Family") in my dataset as a random effect.

Right now, this is my R syntax:
model <- lme(HAPPY ~ Factor1*Factor2, random= ~1|Family, data=mydata, method="REML", na.action=na.omit)
summary(model)
anova(model,type="marginal")

So right now, there is only one random effect. But -- theoretically, we'd expect that identical twins and fraternal twins would differ in their degree of relatedness (the variable is called "twintype"). Should this be reflected in my model, and if so...how? Conceptually, I'm struggling with how it would be specified. And practically, I'm also struggling -- putting the above syntax together was actually quite the accomplishment for me!

2) Heteroscedasticity. What is the best way to assess if heteroscedasticity is affecting my model? Right now I'm just eyeballing --
plot(fitted(model), resid(model)) and looking for cone shapes.

And if heteroscedasticity is an issue, what can I do to fix it? From what I've gathered from the google machine, I can enter a weighting command that will account for the inconstant variance. Could I modify my syntax above to read --
model <- lme(HAPPY ~ Factor1*Factor2, random= ~1|Family, data=mydata, method="REML", na.action=na.omit, weights=varPower())