Skip to content
Prev 314797 / 398503 Next

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:

            
--
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.