VC covariance structure for lme model?
On the off chance, I wanted to push this question once more that someone is familiar with it. I am currently trying to contact asreml but I would like to try and accomplish VC covariance structure with lme if possible. Below are all my attempts. Thanks for any assistance,
dat=read.table("C:/subset.csv",sep=",",header=TRUE, na.strings=".")
attach(dat)
dat34=dat[Group %in% c("3", "4"),]
attach(dat34)
liver34=within(dat34, {
Group=factor(Group)
Event_name=factor(Event_name)
Died=factor(Died)
ID=factor(ID)
})
attach(liver34)
##What is should be from SAS
#CV var
#Type 3 Tests of Fixed Effects
#Effect NumDF DenDF F Value Pr > F
#Group 1 22 0.25 0.6244
#Died 1 22 6.55 0.0179
#Group*Died 1 22 4.43 0.0470
fit.1=lme(var~Group*Died,
random=~1|ID,
data=dat34)
anova(fit.1, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 227.58700 <.0001
#Group 1 22 0.18320 0.6728
#Died 1 22 3.63388 0.0698
#Group:Died 1 22 3.04103 0.0951
fit.2=lme(var~Group*Died,
data=dat34,
random=~1|ID/Died)
anova(fit.2, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 77.99004 <.0001
#Group 1 22 1.46275 0.2393
#Died 1 22 5.84535 0.0243
#Group:Died 1 22 3.04103 0.0951
fit.3=lme(var~Group*Died,
random=list(ID=pdSymm(~Event_name)),
data=dat34)
anova(fit.3, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 273.10918 <.0001
#Group 1 22 0.69692 0.4128
#Died 1 22 1.43316 0.2440
#Group:Died 1 22 5.74399 0.0255
fit.4=lme(var~Group*Died,
random=list(ID=pdSymm(~Group)),
data=dat34)
anova(fit.4, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 235.13889 <.0001
#Group 1 22 0.15878 0.6941
#Died 1 22 3.83253 0.0631
#Group:Died 1 22 3.01222 0.0966
fit.5=lme(var~Group*Died,
random=list(ID=pdSymm(~Group)),
data=dat34,
weights=varIdent(form=~1|Event_name))
anova(fit.5, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 277.16705 <.0001
#Group 1 22 0.23901 0.6298
#Died 1 22 3.99283 0.0582
#Group:Died 1 22 3.23135 0.0860
fit.6=lme(var~Group*Died,
random=list(ID=pdSymm(~Group)),
data=dat34,
weights=varIdent(form=~1|Event_name))
anova(fit.6, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 277.16705 <.0001
#Group 1 22 0.23901 0.6298
#Died 1 22 3.99283 0.0582
#Group:Died 1 22 3.23135 0.0860
fit.7=lme(var~(Group*Died),
random=list(ID=pdCompSymm(~Died)),
data=dat34)
anova(fit.7, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 85.83799 <.0001
#Group 1 22 1.60624 0.2183
#Died 1 22 4.71795 0.0409
#Group:Died 1 22 2.65379 0.1175
fit.8=lme(var~(Group*Died),
data=dat34,
random=~1|ID,
corr=corSymm())
anova(fit.8, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 119.54403 <.0001
#Group 1 22 4.58972 0.0435
#Died 1 22 8.01715 0.0097
#Group:Died 1 22 5.27470 0.0315
fit.9=lme(var~(Group*Died),
data=dat34,
random=list(ID=pdDiag(~Group*Died)),
corr=corSymm(, ~1|ID))
# Error in lme.formula(var ~ (Group * Died), data = dat34, random =
list(ID = pdDiag(~Group * :
# nlminb problem, convergence error code = 1
# message = iteration limit reached without convergence (9)
fit.10=lme(var~(Group*Died),
data=dat34,
random=list(ID=pdDiag(~Group*Died)),
corr=NULL)
anova(fit.10, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 93.90211 <.0001
#Group 1 22 1.75311 0.1991
#Died 1 22 6.84379 0.0158
#Group:Died 1 22 3.11458 0.0915
fit.11=lme(var~Group*Died,
data=dat34,
random=list(ID=pdDiag(~Group*Died)))
anova(fit.11, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 93.90211 <.0001
#Group 1 22 1.75311 0.1991
#Died 1 22 6.84379 0.0158
#Group:Died 1 22 3.11458 0.0915
fit.12=lme(var~Group*Died,
data=dat34,
random=list(ID=pdDiag(~Event_name)))
anova(fit.12, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 87.33040 <.0001
#Group 1 22 1.09661 0.3064
#Died 1 22 5.46329 0.0289
#Group:Died 1 22 2.94589 0.1001
summary(fit.12)
fit.13=lme(var~Group*Died,
data=dat34,
random=list(ID=pdDiag(~Group)))
anova(fit.13, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 77.99004 <.0001
#Group 1 22 1.46275 0.2393
#Died 1 22 5.84535 0.0243
#Group:Died 1 22 3.04103 0.0951
fit.14=lme(var~Group*Died,
data=dat34,
random=list(ID=pdDiag(~Died)))
anova(fit.14, type="marginal", adjustSigma=F)
# numDF denDF F-value p-value
#(Intercept) 1 101 85.83800 <.0001
#Group 1 22 1.60624 0.2183
#Died 1 22 4.71795 0.0409
#Group:Died 1 22 2.65379 0.1175
fit.15=lme(var~Group*Died,
data=dat34,
random=~1|ID,
corr=corCompSymm())
anova(fit.15, type="marginal", adjustSigma=F)
#same as fit.13
fit.16=lme(var~Group*Died,
data=dat34,
random=~1|ID/Event_name)
anova(fit.16, type="marginal", adjustSigma=F)
#same as fit.13
#######################
---------- Forwarded message ---------- From: Charles Determan Jr <deter088 at umn.edu> Date: Thu, May 10, 2012 at 12:59 PM Subject: Re: [R-sig-ME] VC covariance structure for lme model? To: Kevin Wright <kw.stat at gmail.com> Cc: r-sig-mixed-models at r-project.org Thanks for your feedback again Kevin, I did post the data but I will attach it again here. ?I will also look into this asreml package too. Regards, Charles
On Thu, May 10, 2012 at 11:15 AM, Kevin Wright <kw.stat at gmail.com> wrote:
Some comments: 1. Posting the data would help. 2. As I mentioned months ago, looking at the fixed effects tests is generally not satisfactory. ?It's better, IMHO, to make comparisons using the variance structures. 3. You appear to be at UMN. ?If so, you might want to try asreml-r, which is sooooooooo much cleaner and clearer about variance structures. ?You can grab the beta version here: http://www.mmontap.org/tracker Download this one: asreml_3.0-1.zip You'll need to contact VSN for a free license: http://www.vsni.co.uk/software Kevin On Thu, May 10, 2012 at 11:00 AM, Charles Determan Jr <deter088 at umn.edu> wrote:
Greetings again R users,
Some of you will likley recognize me but I hope you can help me once
more. ?I have previously attempted to replicate the UN, CS, and AR(1)
covariance structures used in SAS PROC MIXED. ?However, my efforts
have fallen short on replicating the Variance Components (VC)
structure. ?I have read that it is also known as a diagonal structure.
?Below I have copied over all the models I have tried and their output
with no success. ?Perhaps someone here will see my error or something
I have overlooked. ?I have attached the data for this particular
model. ?Thanks to all, I certainly cannot thank this help list enough.
?I you need any further information/clarification, please ask.
Cheers, Charles
dat=read.table("C:/subset.csv",sep=",",header=TRUE, na.strings=".")
attach(dat)
dat34=dat[Group %in% c("3", "4"),]
attach(dat34)
liver34=within(dat34, {
? ? ? ?Group=factor(Group)
? ? ? ?Event_name=factor(Event_name)
? ? ? ?Died=factor(Died)
? ? ? ?ID=factor(ID)
})
attach(liver34)
##What is should be from SAS
#CV var
#Type 3 Tests of Fixed Effects
#Effect ? ? ? ? ?NumDF DenDF F Value Pr > F
#Group ? ? ? ? ? ? ? ? ? 1 ? ? ?22 ? ? ?0.25 ? ?0.6244
#Died ? ? ? ? ? ? ? ? ? ?1 ? ? ?22 ? ? ?6.55 ? ?0.0179
#Group*Died ? ? ? ? ? ? ?1 ? ? ?22 ? ? ?4.43 ? ?0.0470
fit.1=lme(var~Group*Died,
? ? ? ?random=~1|ID,
? ? ? ?data=dat34)
anova(fit.1, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ? F-value p-value
#(Intercept) ? ?1 ? 101 227.58700 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ? 0.18320 ?0.6728
#Died ? ? ? ? ? 1 ? ?22 ? 3.63388 ?0.0698
#Group:Died ? ? 1 ? ?22 ? 3.04103 ?0.0951
fit.2=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=~1|ID/Died)
anova(fit.2, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 77.99004 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.46275 ?0.2393
#Died ? ? ? ? ? 1 ? ?22 ?5.84535 ?0.0243
#Group:Died ? ? 1 ? ?22 ?3.04103 ?0.0951
fit.3=lme(var~Group*Died,
? ? ? ?random=list(ID=pdSymm(~Event_name)),
? ? ? ?data=dat34)
anova(fit.3, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ? F-value p-value
#(Intercept) ? ?1 ? 101 273.10918 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ? 0.69692 ?0.4128
#Died ? ? ? ? ? 1 ? ?22 ? 1.43316 ?0.2440
#Group:Died ? ? 1 ? ?22 ? 5.74399 ?0.0255
fit.4=lme(var~Group*Died,
? ? ? ?random=list(ID=pdSymm(~Group)),
? ? ? ?data=dat34)
anova(fit.4, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ? F-value p-value
#(Intercept) ? ?1 ? 101 235.13889 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ? 0.15878 ?0.6941
#Died ? ? ? ? ? 1 ? ?22 ? 3.83253 ?0.0631
#Group:Died ? ? 1 ? ?22 ? 3.01222 ?0.0966
fit.5=lme(var~Group*Died,
? ? ? ?random=list(ID=pdSymm(~Group)),
? ? ? ?data=dat34,
? ? ? ?weights=varIdent(form=~1|Event_name))
anova(fit.5, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ? F-value p-value
#(Intercept) ? ?1 ? 101 277.16705 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ? 0.23901 ?0.6298
#Died ? ? ? ? ? 1 ? ?22 ? 3.99283 ?0.0582
#Group:Died ? ? 1 ? ?22 ? 3.23135 ?0.0860
fit.6=lme(var~Group*Died,
? ? ? ?random=list(ID=pdSymm(~Group)),
? ? ? ?data=dat34,
? ? ? ?weights=varIdent(form=~1|Event_name))
anova(fit.6, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ? F-value p-value
#(Intercept) ? ?1 ? 101 277.16705 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ? 0.23901 ?0.6298
#Died ? ? ? ? ? 1 ? ?22 ? 3.99283 ?0.0582
#Group:Died ? ? 1 ? ?22 ? 3.23135 ?0.0860
fit.7=lme(var~(Group*Died),
? ? ? ?random=list(ID=pdCompSymm(~Died)),
? ? ? ?data=dat34)
anova(fit.7, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 85.83799 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.60624 ?0.2183
#Died ? ? ? ? ? 1 ? ?22 ?4.71795 ?0.0409
#Group:Died ? ? 1 ? ?22 ?2.65379 ?0.1175
fit.8=lme(var~(Group*Died),
? ? ? ?data=dat34,
? ? ? ?random=~1|ID,
? ? ? ?corr=corSymm())
anova(fit.8, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ? F-value p-value
#(Intercept) ? ?1 ? 101 119.54403 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ? 4.58972 ?0.0435
#Died ? ? ? ? ? 1 ? ?22 ? 8.01715 ?0.0097
#Group:Died ? ? 1 ? ?22 ? 5.27470 ?0.0315
fit.9=lme(var~(Group*Died),
? ? ? ?data=dat34,
? ? ? ?random=list(ID=pdDiag(~Group*Died)),
? ? ? ?corr=corSymm(, ~1|ID))
# ?Error in lme.formula(var ~ (Group * Died), data = dat34, random =
list(ID = pdDiag(~Group * ?:
# ?nlminb problem, convergence error code = 1
# ?message = iteration limit reached without convergence (9)
fit.10=lme(var~(Group*Died),
? ? ? ?data=dat34,
? ? ? ?random=list(ID=pdDiag(~Group*Died)),
? ? ? ?corr=NULL)
anova(fit.10, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 93.90211 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.75311 ?0.1991
#Died ? ? ? ? ? 1 ? ?22 ?6.84379 ?0.0158
#Group:Died ? ? 1 ? ?22 ?3.11458 ?0.0915
fit.11=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=list(ID=pdDiag(~Group*Died)))
anova(fit.11, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 93.90211 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.75311 ?0.1991
#Died ? ? ? ? ? 1 ? ?22 ?6.84379 ?0.0158
#Group:Died ? ? 1 ? ?22 ?3.11458 ?0.0915
fit.12=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=list(ID=pdDiag(~Event_name)))
anova(fit.12, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 87.33040 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.09661 ?0.3064
#Died ? ? ? ? ? 1 ? ?22 ?5.46329 ?0.0289
#Group:Died ? ? 1 ? ?22 ?2.94589 ?0.1001
summary(fit.12)
fit.13=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=list(ID=pdDiag(~Group)))
anova(fit.13, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 77.99004 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.46275 ?0.2393
#Died ? ? ? ? ? 1 ? ?22 ?5.84535 ?0.0243
#Group:Died ? ? 1 ? ?22 ?3.04103 ?0.0951
fit.14=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=list(ID=pdDiag(~Died)))
anova(fit.14, type="marginal", adjustSigma=F)
# ? ? ? ? ? ? ? numDF denDF ?F-value p-value
#(Intercept) ? ?1 ? 101 85.83800 ?<.0001
#Group ? ? ? ? ?1 ? ?22 ?1.60624 ?0.2183
#Died ? ? ? ? ? 1 ? ?22 ?4.71795 ?0.0409
#Group:Died ? ? 1 ? ?22 ?2.65379 ?0.1175
fit.15=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=~1|ID,
? ? ? ?corr=corCompSymm())
anova(fit.15, type="marginal", adjustSigma=F)
#same as fit.13
fit.16=lme(var~Group*Died,
? ? ? ?data=dat34,
? ? ? ?random=~1|ID/Event_name)
anova(fit.16, type="marginal", adjustSigma=F)
#same as fit.13
#######################
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- Kevin Wright