[R-sig-ME] [R] random effects model

Thu Jan 17 03:04:50 CET 2013

HI,

In the same link at the bottom of the page,

"

All is well now after updating all packages with the following:
update.packages()"

It may or may not solve your problem.

I got your attachments.  You should post those questions in (r-sig-mixed-models at r-project.org).  I suggest you to read lme4 book (http://lme4.r-forge.r-project.org/lMMwR/)
#lrgprt.pdf

A.K.

----- Original Message -----
From: rex2013 <usha.nathan at gmail.com>
To: r-help at r-project.org
Cc: 
Sent: Wednesday, January 16, 2013 5:06 AM
Subject: Re: [R] random effects model

Hi

I tried removing the missing values and installing "plyr". Still error
message appears with ggplot2

Btw, did you get the attachments with my earlier mail?

Ta.

On Wed, Jan 16, 2013 at 3:16 AM, arun kirshna [via R] <
ml-node+s789695n4655612h99 at n4.nabble.com> wrote:

>
>
> Hi,
> Check these links:
> http://comments.gmane.org/gmane.comp.lang.r.ggplot2/6527
> https://groups.google.com/forum/#!msg/ggplot2/nfVjxL0DXnY/5zf50zCeZuMJ
> A.K.
>
> ________________________________
> From: Usha Gurunathan <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=0>>
>
> To: arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=1>>
>
> Cc: R help <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=2>>
>
> Sent: Tuesday, January 15, 2013 6:31 AM
> Subject: Re: [R] random effects model
>
>
> Hi AK
>
> Got an error message with
> library(ggplot2) >
> ggplot(BP.stack1,aes(x=factor(HiBP),fill=Obese))+geom_bar(position="fill")
> Error in rename(x, .base_to_ggplot, warn_missing = FALSE) :  could not find
> function "revalue" >
> ggplot(BP.stack1,aes(x=factor(HiBP),fill=Overweight))+geom_bar(position="fill")
> Error in rename(x, .base_to_ggplot, warn_missing = FALSE) :  could not find
> function "revalue"
> I got the dot plot, thanks for that.
>
> I have attached some plots, not sure how to interpret, they had unusual
> patterns.Is it because of missing data? I tried removing the missing data
> too. They still appeared the same. Do I need to transform the data?
>
>
> Thanks in advance.
>
>
>
>
>
> On Tue, Jan 15, 2013 at 8:54 AM, arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=3>>
> wrote:
>
> HI,
>
> >
> >
> >BP_2b<-read.csv("BP_2b.csv",sep="\t")
> >BP_2bNM<-na.omit(BP_2b)
> >
> >BP.stack3 <-
> reshape(BP_2bNM,idvar="CODEA",timevar="time",sep="",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21"),c("hibp14","hibp21")),v.names=c("Obese","Overweight","HiBP"),times=factor(c(1,2)),direction="long")
>
> >library(car)
> >BP.stack3$Obese<- recode(BP.stack3$Obese,"1='Obese';0='Not Obese'")
> >BP.stack3$Overweight<- recode(BP.stack3$Overweight,"1='Overweight';0='Not
> Overweight'")
> >
> >library(ggplot2)
> >ggplot(BP.stack3,aes(x=factor(HiBP),fill=Obese))+geom_bar(position="fill")
>
> >ggplot(BP.stack3,aes(x=factor(HiBP),fill=Overweight))+geom_bar(position="fill")
>
> >
> >You could try lmer() from lme4.
> >library(lme4)
> >fm1<-lmer(HiBP~time+(1|CODEA), family=binomial,data=BP.stack3) #check
> codes, not sure
> >print(dotplot(ranef(fm1,post=TRUE),
> >              scales = list(x = list(relation = "free")))[[1]])
> >qmt1<- qqmath(ranef(fm1, postVar=TRUE))
> >print(qmt1[[1]])
> >
> >
> >A.K.
> >
> >
> >
> >
> >
> >________________________________
> >From: Usha Gurunathan <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=4>>
>
> >To: arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=5>>
>
> >Cc: R help <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=6>>
>
> >Sent: Monday, January 14, 2013 6:32 AM
> >
> >Subject: Re: [R] random effects model
> >
> >
> >Hi AK
> >
> >I have been trying to create some plots. All being categorical variables,
> I am not getting any luck with plots. The few ones that have worked are
> below:
> >
> >barchart(~table(HiBP)|Obese,data=BP.sub3) ## BP.sub3 is the stacked data
> without missing values
> >
> >barchart(~table(HiBP)|Overweight,data=BP.sub3)
> >
> >plot(jitter(hibp14,factor=2)~jitter(Obese14,factor=2),col="gray",cex=0.7,
> data=Copy.of.BP_2)  ## Copy.of.BP_2 is the original wide format
> >
> >## not producing any good plots with mixed models as well.
> >summary(lme.3 <- lme(HiBP~time, data=BP.sub3,random=~1|CODEA,
> na.action=na.omit))
> >anova(lme.3)
> >head(ranef(lme.3))
> >print(plot(ranef(lme.3))) ##
> >
> >Thanks for any help.
> >
> >
> >
> >
> >
> >On Mon, Jan 14, 2013 at 4:33 AM, arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=7>>
> wrote:
> >
> >
> >>
> >>
> >>HI,
> >>
> >>I think I mentioned to you before that when you reshape the
> >>columns excluding the response variable, response variable gets repeated
> >>(in this case hibp14 or hibp21) and creates the error"
> >>
> >>
> >>I run your code, there are obvious problems in the code so I didn't
> reach up to BP.gee
> >>
> >>
> >>BP_2b<-read.csv("BP_2b.csv",sep="\t")
> >>BP.stack3 <-
> reshape(BP_2b,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),times=factor(c(1,2)),direction="long")
>
> >>
> >>
> >>BP.stack3 <-
> transform(BP.stack3,CODEA=factor(CODEA),Sex=factor(Sex,labels=c("Male","Female")),MaternalAge=factor(MaternalAge,labels=c("39years
> or less","40-49 years","50 years or
> older")),Education=factor(Education,labels=c("Primary/special","Started
> secondary","Completed grade10", "Completed grade12",
> "College","University")),Birthplace=factor(Birthplace,labels=c("Australia","Other
> English-speaking","Other")))
> >>
> >> BP.stack3$Sex <-
> factor(BP.stack3$Sex,levels=levels(BP.stack3$Sex)[c(2,1)])
> >>
> >> BP.sub3a <-  subset(BP.stack3,subset=!(is.na(Sex)| is.na(Education)|
> is.na(Birthplace)|is.na(Education)|is.na(hibp14)| is.na(hibp21)))
> >> nrow(BP.sub3a)
> >>#[1] 3364
> >> BP.sub5a <- BP.sub3a[order(BP.sub3a$CODEA),] # your code was BP.sub5a
> <- BP.sub3a[order(BP.sub5a$CODEA),]
> >>
> ^^^^^ was not defined before
> >>#Next line
> >>BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==1&BPsub3$Obese==0]<-
> "Overweight14"  #It should be BP.sub3 and what is BPsub6, it was not
> defined previously.
> >>#Error in BPsub3$Categ[BPsub6$Overweight == 1 & BPsub3$time == 1 &
> BPsub3$Obese ==  :
> >>  #object 'BPsub3' not found
> >>
> >>
> >>
> >>
> >>
> >>
> >>
> >>A.K.
> >>
> >>
> >>________________________________
> >>From: Usha Gurunathan <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=8>>
>
> >>To: arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=9>>
>
> >>
> >>Sent: Sunday, January 13, 2013 1:51 AM
> >>
> >>Subject: Re: [R] random effects model
> >>
> >>
> >>
> >>HI AK
> >>
> >>Thanks a lot  for explaining that.
> >>
> >>1. With the chi sq. ( in order to find out if the diffce is significant
> between groups) do I have create a separate excel file and make a
> dataframe.How do I go about it?
> >>
> >>I have resent a mail to Jun Yan at a difft email ad( first add.didn't
> work, mail not delivered).
> >>
> >>2. With my previous query ( reg. Obese/Overweight/ Normal at age 14 Vs
> change of blood pressure status at 21), even though I had compromised
> without the age-specific regression, but I am still keen to explore why the
> age-specific regression didn't work despite it looking okay. I have given
> below the syntax. If you get time, could you kindly look at it and see if
> it could work by any chance? Apologies for persisting with this query.
> >>
> >>
> >>BP.stack3 <-
> >>reshape(Copy.of.BP_2,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),times=factor(c(1,2)),direction="long
>
> >>BP.stack3
> >>head(BP.stack3)
> >>tail(BP.stack3)
> >>
> >> names(BP.stack3)[c(2,3,4,5,6,7)] <-
> >>c("Sex","MaternalAge","Education","Birthplace","AggScore","IntScore")
> >>
> >>BP.stack3 <-
> >>transform(BP.stack3,CODEA=factor(CODEA),Sex=factor(Sex,labels=c("Male","Female")),MaternalAge=factor(MaternalAge,labels=c("39years
>
> >>or less","40-49 years","50 years or
> >>older")),Education=factor(Education,labels=c("Primary/special","Started
> >>secondary","Completed grade10", "Completed grade12",
> >>"College","University")),Birthplace=factor(Birthplace,labels=c("Australia","Other
>
> >>English-speaking","Other")))
> >>
> >>table(BP.stack3$Sex)
> >>BP.stack3$Sex <-
> >>factor(BP.stack3$Sex,levels=levels(BP.stack3$Sex)[c(2,1)])
> >>
> >>levels(BP.stack3$Sex)
> >>BP.sub3a <-  subset(BP.stack3,subset=!(is.na(Sex)| is.na(Education)|
> is.na(Birthplace)|is.na(Education)|is.na(hibp14)| is.na(hibp21)))
> >>summary(BP.sub3a)
> >>BP.sub5a <- BP.sub3a[order(BP.sub5a$CODEA),]
> >> BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==1&BPsub3$Obese==0]
> >><- "Overweight14"
> >>BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==2&BPsub3$Obese==0]
> >><- "Overweight21"
> >>BPsub3$Categ[BPsub3$Obese==1&BPsub3$time==1&BPsub3$Overweight==0|BPsub3$Obese==1&BPsub3$time==1&BPsub3$Overweight==1
>
> >>] <- "Obese14"
> >>BPsub3$Categ[BPsub3$Obese==0&BPsub3$time==1&BPsub3
> BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==1&BPsub3$Obese==0]
> >><- "Overweight14"$Overweight==0]
> >>
> >><- "Normal14"
> >>BPsub3$Categ[BPsub3$Obese==0&BPsub3$time==2&BPsub3$Overweight==0]
> >><- "Normal21"
> >>BPsub3$Categ[BPsub3$Obese==1&BPsub3$time==2&BPsub3$Overweight==0|BPsub3$Obese==1&BPsub3$time==2&BPsub3$Overweight==1]
>
> >><- "Obese21"
> >>
> >>
> >>
> >>BPsub3$Categ <- factor(BPsub3$Categ)
> >>BPsub3$time <- factor(BPsub3$time)
> >>summary(BPsub3$Categ)
> >>BPsub7 <- subset(BPsub6,subset=!(is.na(Categ)))
> >>BPsub7$time <-
> >>recode(BPsub7$time,"1=14;2=21")
> >>BPsub7$hibp14 <- factor(BPsub7$hibp14)
> >>levels(BPsub7$hibp14)
> >>levels(BPsub7$Categ)
> >>names(BPsub7)
> >>head(BPsub7)    ### this was looking quite okay.
> >>
> >>tail(BPsub7)
> >>str(BPsub7)
> >>
> >>library(gee)
> >>
> >>BP.gee <- geese(hibp14~ time*Categ,
> >>data=BPsub7,id=CODEA,family=binomial,
> >>corstr="exchangeable",na.action=na.omit)
> >>
> >>Thanks again.
> >>
> >>
> >>
> >>On Sun, Jan 13, 2013 at 1:22 PM, arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=10>>
> wrote:
> >>
> >>HI,
> >>>
> >>>table(BP_2b$Sex) #original dataset
> >>>#   1    2
> >>>#3232 3028
> >>> nrow(BP_2b)
> >>>#[1] 6898
> >>> nrow(BP_2bSexNoMV)
> >>>#[1] 6260
> >>> 6898-6260
> >>>#[1] 638 #these rows were removed from the BP_2b to create BP_2bSexNoMV
> >>>BP_2bSexMale<-BP_2bSexNoMV[BP_2bSexNoMV$Sex=="Male",]
> >>> nrow(BP_2bSexMale)
> >>>#[1] 3232
> >>> nrow(BP_2bSexMale[!complete.cases(BP_2bSexMale),]) #Missing rows with
> Male
> >>>#[1] 2407
> >>> nrow(BP_2bSexMale[complete.cases(BP_2bSexMale),]) #Non missing rows
> with Male
> >>>#[1] 825
> >>>
> >>>
> >>>You did the chisquare test on the new dataset with 6260 rows, right.
> >>>I removed those 638 rows because these doesn't belong to either male or
> female, but you want the % of missing value per male or female.  So, I
> thought this will bias the results.  If you want to include the missing
> values, you could do it, but I don't know where you would put that missing
> values as it cannot be classified as belonging specifically to males or
> females.  I hope you understand it.
> >>>
> >>>Sometimes, the maintainer's respond a bit slow.  You have to sent an
> email reminding him again.
> >>>
> >>>Regarding the vmv package, you could email Waqas Ahmed Malik ([hidden
> email] <http://user/SendEmail.jtp?type=node&node=4655612&i=11>) regarding
> options for changing the title and the the font etc.
> >>>You could also use this link (
> http://www.r-bloggers.com/visualizing-missing-data-2/ ) to plot missing
> value (?plot.missing()).  I never used that package, but you could try.
> Looks like it gives more information.
> >>>
> >>>A.K.
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>________________________________
> >>>From: Usha Gurunathan <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=12>>
>
> >>>To: arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=13>>
>
> >>>Sent: Saturday, January 12, 2013 9:05 PM
> >>>
> >>>Subject: Re: [R] random effects model
> >>>
> >>>
> >>>Hi A.K
> >>>
> >>>So it is number of females missing/total female participants enrolled:
> 72.65%
> >>>Number of females missing/total (of males+ females)  participants
> enrolled : 35.14%
> >>>
> >>>The total no. with the master data: Males: 3232, females: 3028 ( I got
> this before removing any missing values)
> >>>
> >>>with table(Copy.of.BP_2$ Sex)  ## BP
> >>>
> >>>
> >>>If I were to write a table (  and do a chi sq. later),
> >>>
> >>>as Gender            Study                    Non study/missing
> Total
> >>>      Male              825 (25.53%)             2407 (74.47%)
> 3232 (100%)
> >>>    Female           828 (27.35%)             2200 (72.65%)       3028
> ( 100%)
> >>>     Total              1653                          4607
>          6260
> >>>
> >>>
> >>>The problem is when I did
> >>>>colSums(is.na(Copy.of.BP_2), the sex category showed N=638.
> >>>
> >>>I cannot understand the discrepancy.Also, when you have mentioned to
> remove NA, is that not a missing value that needs to be included in the
> total number missing. I am a bit confused. Can you help?
> >>>
> >>>## I tried sending email to gee pack maintainer at the ID with R site,
> mail didn't go through??
> >>>
> >>>Many thanks
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>On Sun, Jan 13, 2013 at 9:17 AM, arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=14>>
> wrote:
> >>>
> >>>Hi,
> >>>>Yes, you are right.  72.655222% was those missing among females.
>  35.14377% of values in females are missing from among the whole dataset
> (combined total of Males+Females data after removing the NAs from the
> variable "Sex").
> >>>>
> >>>>A.K.
> >>>>
> >>>>
> >>>>
> >>>>________________________________
> >>>>From: Usha Gurunathan <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=15>>
>
> >>>>To: arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=16>>
>
> >>>>Cc: R help <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=17>>
>
> >>>>Sent: Saturday, January 12, 2013 5:59 PM
> >>>>
> >>>>Subject: Re: [R] random effects model
> >>>>
> >>>>
> >>>>
> >>>>Hi AK
> >>>>That works. I was trying to get  similar results from any other
> package. Being a beginner, I was not sure how to modify the syntax to get
> my output.
> >>>>
> >>>>lapply(split(BP_2bSexNoMV,BP_
> >>>>2bSexNoMV$Sex),function(x)
> (nrow(x[!complete.cases(x[,-2]),])/nrow(x))*100) #gives the percentage of
> rows of missing #values from the overall rows for Males and Females
> >>>>#$Female
> >>>>#[1] 72.65522
> >>>>#
> >>>>#$Male
> >>>>#[1] 74.47401
> >>>>
> >>>>#iF you want the percentage from the total number rows in Males and
> Females (without NA's in the the Sex column)
> >>>> lapply(split(BP_2bSexNoMV,BP_2bSexNoMV$Sex),function(x)
> (nrow(x[!complete.cases(x[,-2]),])/nrow(BP_2bSexNoMV))*100)
> >>>>#$Female
> >>>>#[1] 35.14377
> >>>>#
> >>>>#$Male
> >>>>#[1] 38.45048
> >>>>
> >>>>How do I interpret the above 2 difft results? 72.66% of values were
> missing among female participants?? Can you pl. clarify.
> >>>>
> >>>>Many thanks.
> >>>>
> >>>>
> >>>>On Sun, Jan 13, 2013 at 3:28 AM, arun <[hidden email]<http://user/SendEmail.jtp?type=node&node=4655612&i=18>>
> wrote:
> >>>>
> >>>>lapply(split(BP_2bSexNoMV,BP_2bSexNoMV$Sex),function(x)
> (nrow(x[!complete.cases(x[,-2]),])/nrow(x))*100) #gives the percentage of
> rows of missing #values from the overall rows for Males and Females
> >>>>>#$Female
> >>>>>#[1] 72.65522
> >>>>>#
> >>>>>#$Male
> >>>>>#[1] 74.47401
> >>>>>
> >>>>>#iF you want the percentage from the total number rows in Males and
> Females (without NA's in the the Sex column)
> >>>>> lapply(split(BP_2bSexNoMV,BP_2bSexNoMV$Sex),function(x)
> (nrow(x[!complete.cases(x[,-2]),])/nrow(BP_2bSexNoMV))*100)
> >>>>>#$Female
> >>>>>#[1] 35.14377
> >>>>>#
> >>>>>#$Male
> >>>>>#[1] 38.45048
> >>>>
> >>>
> >>
> >
>
> ______________________________________________
> [hidden email] <http://user/SendEmail.jtp?type=node&node=4655612&i=19>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-tp4654346p4655612.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-tp4654346p4655704.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.