[R-sig-ME] VC covariance structure for lme model?
Charles Determan Jr
deter088 at umn.edu
Thu May 10 19:59:21 CEST 2012
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
More information about the R-sig-mixed-models
mailing list