random effects model
HI,
Regarding question:2) Have you checked summary(m1)?
data(seizure)
???? ## Diggle, Liang, and Zeger (1994) pp166-168, compare Table 8.10
???? 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)
#Call:
#geese(formula = y ~ offset(log(t)) + x + trt + x:trt, id = id,
?# ? data = seiz.l, family = poisson, corstr = "exch")
#
#Mean Model:
# Mean Link:???????????????? log
# Variance to Mean Relation: poisson
?#Coefficients:
? # ??????????? estimate??? san.se?????? wald???????? p??????????????????
#(Intercept)? 1.34760922 0.1573571 73.3423807 0.0000000
#x??????????? 0.11183602 0.1159304? 0.9306116 0.3347040
#trt????????? 0.02753449 0.2217878? 0.0154127 0.9011982
#x:trt?????? -0.10472579 0.2134448? 0.2407334 0.6236769
--------------------------------------------------------------------
#p is the colname with the p values for the Coefficients
summary(m1)$mean["p"]
#??????????????????? p
#(Intercept) 0.0000000
#x?????????? 0.3347040
#trt???????? 0.9011982
#x:trt?????? 0.6236769
You didn't say specifically whether you got NA's in example data or your actual data.? I am getting the p-values in R 2.15.
1) I think it should work with binary response variable.
Another example in the same documentation:
?data(ohio)
str(ohio)
#'data.frame':??? 2148 obs. of? 4 variables:
# $ resp : int? 0 0 0 0 0 0 0 0 0 0 ...
# $ id?? : int? 0 0 0 0 1 1 1 1 2 2 ...
# $ age? : int? -2 -1 0 1 -2 -1 0 1 -2 -1 ...
# $ smoke: int? 0 0 0 0 0 0 0 0 0 0 ...
It is not even factors
???? fit <- geese(resp ~ age + smoke + age:smoke, id=id, data=ohio,
????????????????? family=binomial, corstr="exch", scale.fix=TRUE)
?summary(fit)$mean["p"]
?# ?????????????????? p
#(Intercept) 0.00000000
#age???????? 0.01523698
#smoke?????? 0.09478252
#age:smoke?? 0.42234200
# also tested after converting to factors
ohio1<-ohio
ohio1$smoke<-factor(ohio1$smoke)
?ohio1$age<-factor(ohio1$age)
str(ohio1)
#'data.frame':??? 2148 obs. of? 4 variables:
# $ resp : int? 0 0 0 0 0 0 0 0 0 0 ...
# $ id?? : int? 0 0 0 0 1 1 1 1 2 2 ...
# $ age? : Factor w/ 4 levels "-2","-1","0",..: 1 2 3 4 1 2 3 4 1 2 ...
# $ smoke: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
fit1<-geese(resp~age+smoke+age:smoke,id=id,data=ohio1,family=binomial,corstr="exch",scale.fix=TRUE)
?summary(fit1)$mean["p"]
#????????????????????? p
#(Intercept)? 0.00000000
#age-1??????? 0.60555454
#age0???????? 0.45322698
#age1???????? 0.01187725
#smoke1?????? 0.86262269
#age-1:smoke1 0.17239050
#age0:smoke1? 0.32223942
#age1:smoke1? 0.36686706
A.K.
----- Original Message -----
From: rex2013 <usha.nathan at gmail.com>
To: r-help at r-project.org
Cc:
Sent: Monday, January 7, 2013 6:15 AM
Subject: Re: [R] random effects model
Hi A.K
Below is the comment I get, not sure why.
BP.sub3 is the stacked data without the missing values.
BP.geese3 <- geese(HiBP~time*MaternalAge,data=BP.sub3,id=CODEA,
family=binomial, corstr="unstructured", na.action=na.omit)Error in
`contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
? contrasts can be applied only to factors with 2 or more levels
Even though age has 3 levels; time has 14 years & 21 years; HIBP is a
binary response outcome.
2) When you mentioned summary(m1)$mean["p"] what did the p mean? i
used this in one of the gee command, it produced NA as answer?
Many thanks
On Mon, Jan 7, 2013 at 5:26 AM, arun kirshna [via R] <
ml-node+s789695n4654795h72 at n4.nabble.com> wrote:
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 <[hidden email]<http://user/SendEmail.jtp?type=node&node=4654795&i=0>>
To: [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=1>
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] <
[hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=2>>
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:
. 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]]
______________________________________________ [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=3>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. ______________________________________________ [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=4>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. ------------------------------ ? If you reply to this email, your message will be added to the discussion below: http://r.789695.n4.nabble.com/random-effects-model-tp4654346p4654795.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-tp4654346p4654826.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.