[R] Covariance structure for lme
Charles Determan Jr
deter088 at umn.edu
Thu May 17 18:03:44 CEST 2012
Greetings again R users,
Some of you will likely recognize me but I hope you can help me once
more. I have tried the mixed model mailing list for this question but have
yet to find a solution. As such I hope someone will have another idea.
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)
contrasts(liver34$Old_Event_name)=contr.sum(n=6)
contrasts(liver34$Pig_group)=contr.sum(n=2)
contrasts(liver34$Died)=contr.sum(n=2)
##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
#######################
More information about the R-help
mailing list