[R] random effects model
arun
smartpink111 at yahoo.com
Sun Jan 6 20:23:05 CET 2013
Hi,
I am not very familiar with the geese/geeglm(). Is it from library(geepack)?
Regarding your question:
"
Can you tell me if I can use the geese or geeglm function with this data
eg: : HIBP~ time* Age
Here age is a factor with 3 levels, time: 2 levels, HIBP = yes/no.
From your original data:
BP_2b<-read.csv("BP_2b.csv",sep="\t")
head(BP_2b,2)
# CODEA Sex MaternalAge Education Birthplace AggScore IntScore Obese14
#1 1 NA 3 4 1 NA NA NA
#2 3 2 3 3 1 0 0 0
# Overweight14 Overweight21 Obese21 hibp14 hibp21
#1 NA NA NA NA NA
#2 0 1 0 0 0
If I understand your new classification:
BP.stacknormal<- subset(BP_2b,Obese14==0 & Overweight14==0 & Obese21==0 & Overweight21==0)
BP.stackObese <- subset(BP_2b,(Obese14==1& Overweight14==0 & Obese14==1&Overweight14==1)|(Obese14==1&Overweight14==1 & Obese21==1 & Overweight21==0)|(Obese14==1&Overweight14==0 & Obese21==0 & Overweight21==0)|(Obese14==0 & Overweight14==0 & Obese21==1 & Overweight21==0)|(Obese14==0 & Overweight14==0 & Obese21==1 & Overweight21==1)|(Obese14==0 & Overweight14==1 & Obese21==1 &Overweight21==1)|(Obese14==1& Overweight14==1 & Obese21==1& Overweight21==1)) #check whether there are more classification that fits to #Obese
BP.stackOverweight <- subset(BP_2b,(Obese14==0 & Overweight14==1 & Obese21==0 & Overweight21==1)|(Obese14==0 &Overweight14==1 & Obese21==0 & Overweight21==0)|(Obese14==0 & Overweight14==0 & Obese21==0 & Overweight21==1))
BP.stacknormal$Categ<-"Normal"
BP.stackObese$Categ<-"Obese"
BP.stackOverweight$Categ <- "Overweight"
BP.newObeseOverweightNormal<-na.omit(rbind(BP.stacknormal,BP.stackObese,BP.stackOverweight))
nrow(BP.newObeseOverweightNormal)
#[1] 1581
BP.stack3 <- reshape(BP.newObeseOverweightNormal,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21"),c("hibp14","hibp21")),v.names=c("Obese","Overweight","hibp"),direction="long")
library(car)
BP.stack3$time<-recode(BP.stack3$time,"1=14;2=21")
head(BP.stack3,2)
# CODEA Sex MaternalAge Education Birthplace AggScore IntScore Categ time
#8.1 8 2 4 4 1 0 0 Normal 14
#9.1 9 1 3 6 2 0 0 Normal 14
# Obese Overweight hibp
#8.1 0 0 0
Now, your formula: (HIBP~time*Age), is it MaternalAge?
If it is, it has three values
unique(BP.stack3$MaternalAge)
#[1] 4 3 5
and for time (14,21) # If it says that geese/geeglm, contrasts could be applied with factors>=2 levels, what is the problem?
If you take "Categ" variable, it also has 3 levels (Normal, Obese, Overweight).
BP.stack3$MaternalAge<-factor(BP.stack3$MaternalAge)
BP.stack3$time<-factor(BP.stack3$time)
library(geepack)
For your last question about how to get the p-values:
# Using one of the example datasets:
data(seizure)
seiz.l <- reshape(seizure,
varying=list(c("base","y1", "y2", "y3", "y4")),
v.names="y", times=0:4, direction="long")
seiz.l <- seiz.l[order(seiz.l$id, seiz.l$time),]
seiz.l$t <- ifelse(seiz.l$time == 0, 8, 2)
seiz.l$x <- ifelse(seiz.l$time == 0, 0, 1)
m1 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
data=seiz.l, corstr="exch", family=poisson)
summary(m1)
summary(m1)$mean["p"]
# p
#(Intercept) 0.0000000
#x 0.3347040
#trt 0.9011982
#x:trt 0.6236769
#If you need the p-values of the scale
summary(m1)$scale["p"]
# p
#(Intercept) 0.0254634
Hope it helps.
A.K.
----- Original Message -----
From: rex2013 <usha.nathan at gmail.com>
To: r-help at r-project.org
Cc:
Sent: Sunday, January 6, 2013 4:55 AM
Subject: Re: [R] random effects model
Hi A.K
Regarding my question on comparing normal/ obese/overweight with blood
pressure change, I did finally as per the first suggestion of stacking the
data and creating a normal category . This only gives me a obese not obese
14, but when I did with the wide format hoping to get a
obese14,normal14,overweight 14 Vs hibp 21, i could not complete any of the
models.
This time I classified obese=1 & overweight=1 as obese itself.
Can you tell me if I can use the geese or geeglm function with this data
eg: : HIBP~ time* Age
Here age is a factor with 3 levels, time: 2 levels, HIBP = yes/no.
It says geese/geeglm: contrast can be applied only with factor with 2 or
more levels. What is the way to overcome this. Can I manipulate the data to
make it work.
I need to know if the demogrphic variables affect change in blood pressure
status over time?
How to get the p values with gee model?
Thanks
On Thu, Jan 3, 2013 at 5:06 AM, arun kirshna [via R] <
ml-node+s789695n4654438h5 at n4.nabble.com> wrote:
> HI Rex,
> If I take a small subset from your whole dataset, and go through your
> codes:
> BP_2b<-read.csv("BP_2b.csv",sep="\t")
> BP.sub<-BP_2b[410:418,c(1,8:11,13)] #deleted the columns that are not
> needed
> BP.stacknormal<- subset(BP.subnew,Obese14==0 & Overweight14==0)
> BP.stackObese <- subset(BP.subnew,Obese14==1)
> BP.stackOverweight <- subset(BP.subnew,Overweight14==1)
> BP.stacknormal$Categ<-"Normal14"
> BP.stackObese$Categ<-"Obese14"
> BP.stackOverweight$Categ <- "Overweight14"
> BP.newObeseOverweightNormal<-rbind(BP.stacknormal,BP.stackObese,BP.stackOverweight)
>
> BP.newObeseOverweightNormal
> # CODEA Obese14 Overweight14 Overweight21 Obese21 hibp21 Categ
> #411 541 0 0 0 0 0 Normal14
> #415 545 0 0 1 1 1 Normal14
> #418 549 0 0 1 0 0 Normal14
> #413 543 1 0 1 1 0 Obese14
> #417 548 0 1 1 0 0 Overweight14
> BP.newObeseOverweightNormal$Categ<-
> factor(BP.newObeseOverweightNormal$Categ)
> BP.stack3 <-
> reshape(BP.newObeseOverweightNormal,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),direction="long")
>
> library(car)
> BP.stack3$time<-recode(BP.stack3$time,"1=14;2=21")
> BP.stack3 #Here Normal14 gets repeated even at time==21. Given that you
> are using the "Categ" and "time" #columns in the analysis, it will give
> incorrect results.
> # CODEA hibp21 Categ time Obese Overweight
> #541.1 541 0 Normal14 14 0 0
> #545.1 545 1 Normal14 14 0 0
> #549.1 549 0 Normal14 14 0 0
> #543.1 543 0 Obese14 14 1 0
> #548.1 548 0 Overweight14 14 0 1
> #541.2 541 0 Normal14 21 0 0
> #545.2 545 1 Normal14 21 1 1
> #549.2 549 0 Normal14 21 0 1
> #543.2 543 0 Obese14 21 1 1
> #548.2 548 0 Overweight14 21 0 1
> #Even if I correct the above codes, this will give incorrect
> results/(error as you shown) because the response variable (hibp21) gets
> #repeated when you reshape it from wide to long.
>
> The correct classification might be:
> BP_2b<-read.csv("BP_2b.csv",sep="\t")
> BP.sub<-BP_2b[410:418,c(1,8:11,13)]
> BP.subnew<-reshape(BP.sub,idvar="CODEA",timevar="time",sep="",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),direction="long")
>
> BP.subnew$time<-recode(BP.subnew$time,"1=14;2=21")
> BP.subnew<-na.omit(BP.subnew)
>
> BP.subnew$Categ[BP.subnew$Overweight==1 & BP.subnew$time==14 &
> BP.subnew$Obese==0]<-"Overweight14"
> BP.subnew$Categ[BP.subnew$Overweight==1 & BP.subnew$time==21 &
> BP.subnew$Obese==0]<-"Overweight21"
> BP.subnew$Categ[BP.subnew$Obese==1 & BP.subnew$time==14 &
> BP.subnew$Overweight==0]<-"Obese14"
> BP.subnew$Categ[BP.subnew$Obese==1 & BP.subnew$time==21 &
> BP.subnew$Overweight==0]<-"Obese21"
> BP.subnew$Categ[BP.subnew$Overweight==1 & BP.subnew$time==21&
> BP.subnew$Obese==1]<-"ObeseOverweight21"
> BP.subnew$Categ[BP.subnew$Overweight==1 & BP.subnew$time==14&
> BP.subnew$Obese==1]<-"ObeseOverweight14"
> BP.subnew$Categ[BP.subnew$Overweight==0 & BP.subnew$Obese==0
> &BP.subnew$time==14]<-"Normal14"
> BP.subnew$Categ[BP.subnew$Overweight==0 & BP.subnew$Obese==0
> &BP.subnew$time==21]<-"Normal21"
>
> BP.subnew$Categ<-factor(BP.subnew$Categ)
> BP.subnew$time<-factor(BP.subnew$time)
> BP.subnew
> # CODEA hibp21 time Obese Overweight Categ
> #541.1 541 0 14 0 0 Normal14
> #543.1 543 0 14 1 0 Obese14
> #545.1 545 1 14 0 0 Normal14
> #548.1 548 0 14 0 1 Overweight14
> #549.1 549 0 14 0 0 Normal14
> #541.2 541 0 21 0 0 Normal21
> #543.2 543 0 21 1 1 ObeseOverweight21
> #545.2 545 1 21 1 1 ObeseOverweight21
> #548.2 548 0 21 0 1 Overweight21
> #549.2 549 0 21 0 1 Overweight21
>
> #NOw with the whole dataset:
> BP.sub<-BP_2b[,c(1,8:11,13)] #change here and paste the above lines:
> head(BP.subnew)
> # CODEA hibp21 time Obese Overweight Categ
> #3.1 3 0 14 0 0 Normal14
> #7.1 7 0 14 0 0 Normal14
> #8.1 8 0 14 0 0 Normal14
> #9.1 9 0 14 0 0 Normal14
> #14.1 14 1 14 0 0 Normal14
> #21.1 21 0 14 0 0 Normal14
>
> tail(BP.subnew)
> # CODEA hibp21 time Obese Overweight Categ
> #8485.2 8485 0 21 1 1 ObeseOverweight21
> #8506.2 8506 0 21 0 1 Overweight21
> #8520.2 8520 0 21 0 0 Normal21
> #8529.2 8529 1 21 1 1 ObeseOverweight21
> #8550.2 8550 0 21 1 1 ObeseOverweight21
> #8554.2 8554 0 21 0 0 Normal21
>
> summary(lme.1 <- lme(hibp21~time+Categ+ time*Categ,
> data=BP.subnew,random=~1|CODEA, na.action=na.omit))
> #Error in MEEM(object, conLin, control$niterEM) :
> #Singularity in backsolve at level 0, block 1
> #May be because of the reasons I mentioned above.
>
> #YOu didn't mention the library(gee)
> BP.gee8 <- gee(hibp21~time+Categ+time*Categ,
> data=BP.subnew,id=CODEA,family=binomial,
> corstr="exchangeable",na.action=na.omit)
> #Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
> #Error in gee(hibp21 ~ time + Categ + time * Categ, data = BP.subnew, id =
> CODEA, :
> #rank-deficient model matrix
> With your codes, it might have worked, but the results may be inaccurate
> # After running your whole codes:
> BP.gee8 <- gee(hibp21~time+Categ+time*Categ,
> data=BP.stack3,id=CODEA,family=binomial,
> corstr="exchangeable",na.action=na.omit)
> #Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
> #running glm to get initial regression estimate
> # (Intercept) time CategObese14
> # -2.456607e+01 9.940875e-15 2.087584e-13
> # CategOverweight14 time:CategObese14 time:CategOverweight14
> # 2.087584e-13 -9.940875e-15 -9.940875e-15
> #Error in gee(hibp21 ~ time + Categ + time * Categ, data = BP.stack3, id =
> CODEA, :
> # Cgee: error: logistic model for probability has fitted value very close
> to 1.
> #estimates diverging; iteration terminated.
>
> In short, I think it would be better to go with the suggestion in my
> previous email with adequate changes in "Categ" variable (adding
> ObeseOverweight14, ObeseOverweight21 etc) as I showed here.
>
> A.K.
>
>
>
>
>
>
>
>
> ------------------------------
> If you reply to this email, your message will be added to the discussion
> below:
> http://r.789695.n4.nabble.com/random-effects-model-tp4654346p4654438.html
> To unsubscribe from random effects model, click here<http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=4654346&code=dXNoYS5uYXRoYW5AZ21haWwuY29tfDQ2NTQzNDZ8MjAyMjE1NDI0>
> .
> NAML<http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>
--
View this message in context: http://r.789695.n4.nabble.com/random-effects-model-tp4654346p4654776.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list