[R-sig-ME] lme heteroscedasticity
Joshua Wiley
jwiley.psych at gmail.com
Sun Mar 31 23:19:09 CEST 2013
Dear Iasonas,
Consider this simple example:
mtcars <- mtcars
mtcars$weights <- sample(0:4, size = nrow(mtcars), replace=TRUE)
lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~weights), data = mtcars)
which gives many errors. But...
lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~(weights+1)),
data = mtcars)
and even
lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~(weights+.001)),
data = mtcars)
but not
lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~(weights+0)),
data = mtcars)
Anyway, the short of it is that having 0 weights does not play well
with lme. They do not have to be far away from zero, but it does not
like them to be exactly zero.
In your case, you have count of days absent, with most students having
zero. If that is the baseline, and then there is more variability in
the students who are absent at least some days, since the weights are
a multiplier essentially of the residual variance, consider adding 1,
which then would make the residual variance for no absent, and then it
would increase based on number of days absent.
Does that make sense?
Cheers,
Joshua
On Sun, Mar 31, 2013 at 1:07 AM, Iasonas Lamprianou
<lamprianou at yahoo.com> wrote:
>
>
> Dear all, a few days ago I sent an email to this list and got some responses (I thank you all). This time, I am trying to fit an heteroscedastic model on lme:
>
> model0<-lme(score ~ 1 + overallabs +PARENTALEDUCATIONALLEVEL + PARENTALOCCUPATIONALLEVEL + GENDER + SUBJECT * ETHNICITY + SCHOOL2 ,random=~1|id, weights=varFixed(~overallabs), data=galatia.comb_small.vert)
>
> overallabs is the number of absences from school (which is a heavily skewed variable; most students have zero absences) and this causes some heteroscedasticity in the model. id is the id of the students (each one has scores on 4 subjects which is modelled as fixed effects)
>
>
> However, when I try this model, I get the result:
>
> Error in MEestimate(lmeSt, grps) :
> NA/NaN/Inf in foreign function call (arg 1)
>
> Not sure why this is the case. Anyone any ideas? The internet cannot help me in this case....
>
> I am using the laterst nlme version on linux mint.
>
> Thanks
>
>
> Dr. Iasonas Lamprianou
> Department of Social and Political Sciences
> University of Cyprus
>
>
>>________________________________
>> From: "r-sig-mixed-models-request at r-project.org" <r-sig-mixed-models-request at r-project.org>
>>To: r-sig-mixed-models at r-project.org
>>Sent: Friday, 29 March 2013, 8:23
>>Subject: R-sig-mixed-models Digest, Vol 75, Issue 45
>>
>>Send R-sig-mixed-models mailing list submissions to
>> r-sig-mixed-models at r-project.org
>>
>>To subscribe or unsubscribe via the World Wide Web, visit
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>or, via email, send a message with subject or body 'help' to
>> r-sig-mixed-models-request at r-project.org
>>
>>You can reach the person managing the list at
>> r-sig-mixed-models-owner at r-project.org
>>
>>When replying, please edit your Subject line so it is more specific
>>than "Re: Contents of R-sig-mixed-models digest..."
>>
>>
>>Today's Topics:
>>
>> 1. how to impose same variance for two random effects with fixed
>> correlation patterns? (Gabriel Baud-Bovy)
>> 2. Re: glmmADMB troubles [diagnosed + worked around] (Ross Boylan)
>> 3. glmmADMB truncated distributions (Ross Boylan)
>> 4. Re: glmmADMB truncated distributions (Ben Bolker)
>> 5. confidence interval for the SD of random effects
>> (Iasonas Lamprianou)
>> 6. variance of fixed effects (Iasonas Lamprianou)
>>
>>
>>----------------------------------------------------------------------
>>
>>Message: 1
>>Date: Thu, 28 Mar 2013 21:14:47 +0100
>>From: Gabriel Baud-Bovy <baud-bovy.gabriel at hsr.it>
>>To: r-sig-mixed-models at r-project.org
>>Subject: [R-sig-ME] how to impose same variance for two random effects
>> with fixed correlation patterns?
>>Message-ID: <5154A4B7.2070702 at hsr.it>
>>Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>
>>Dear all,
>>
>>Pinheiro & Bates (1996) give an example of use of pdBlocked to represent
>>a two-level mixed-effects model as a single-level model (p. 162):
>>
>># two-level model
>>lme(yield ~ nitro, data = Oats, random = list(Block=pdIdent(~1),
>>Variety=pdIdent(~1)) )
>>
>># single-level model
>>lme(yield ~ nitro, data = Oats,
>> random = list(Block=pdBlocked(list(pdIdent(~1),pdIdent(~Variety-1))))
>>
>>This yields a diagonal covariance matrix with variance sigma_1^2 in the
>>first block
>>and variance sigma_2^2 in the second block.
>>
>>I would like to have the SAME variance in both blocks (this would be a way
>>to impose that the two variances in the corresponding two-level model
>>are equal).
>>
>>I have looked at pdConstruct.pdBlocked to see whether it was
>>possible to impose this restriction but the code is daunting and
>>I would appreciate any help about how to proceed.
>>
>>Best,
>>
>>Gabriel
>>
>>P.S. Let me know if my question (or the previous one about
>>how to make the covariance matrix for random effects depend
>>on some characteristics of the grouping factor) is more
>>appropriate at R-help .
>>
>>
>>
>>--
>>---------------------------------------------------------------------
>>Gabriel Baud-Bovy tel.: (+39) 02 2643 4839 (office)
>>UHSR University (+39) 02 2643 3429 (laboratory)
>>via Olgettina, 58 (+39) 02 2643 4891 (secretary)
>>20132 Milan, Italy fax: (+39) 02 2643 4892
>>
>>
>>
>>------------------------------
>>
>>Message: 2
>>Date: Thu, 28 Mar 2013 15:22:35 -0700
>>From: Ross Boylan <ross at biostat.ucsf.edu>
>>To: Ben Bolker <bbolker at gmail.com>
>>Cc: r-sig-mixed-models at r-project.org
>>Subject: Re: [R-sig-ME] glmmADMB troubles [diagnosed + worked around]
>>Message-ID: <5154C2AB.5000502 at biostat.ucsf.edu>
>>Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>
>>I stepped through the code and think I found the problem. Since I'm
>>running under emacs/ESS the value of the environment variable SHELL is
>>funny:
>>
>>[context: run_bin does think it's in windows and invokes
>> shell(cmd, invisible = TRUE, intern = !verbose)]
>>
>>The relevant branch of shell() is (note it uses shell as a variable
>>inside the function)
>> if (missing(shell)) {
>> shell <- Sys.getenv("R_SHELL")
>> if (!nzchar(shell))
>> shell <- Sys.getenv("SHELL")
>> if (!nzchar(shell))
>> shell <- Sys.getenv("COMSPEC")
>>R_SHELL is empty
>>SHELL is "C:/Program Files/GNU Emacs 24.2/bin/cmdproxy.exe"
>>I believe that is where the "C:/Program" was coming from in the error
>>messages.
>>COMSPEC is "C:\\Windows\\system32\\cmd.exe".
>>
>>I set R_SHELL to the value of COMSPEC and now the command runs.
>>
>>Ross
>>
>>P.S. shell() in trun invokes system(), which is why you were seeing it
>>in the traceback.
>>
>>On 3/6/2013 6:37 PM, Ross Boylan wrote:
>>> On 3/2/2013 3:04 PM, Ben Bolker wrote:
>>>> On 13-02-28 01:15 PM, Ross Boylan wrote:
>>>>> On 2/27/2013 8:05 PM, Ben Bolker wrote:
>>>>>> Ross Boylan <ross at ...> writes:
>>>> [snip: not including version without NAs removed, since we already
>>>> know what the issue is there]
>>>>
>>>>>> r <- glmmadmb(sexActs~(1|id), sexpartner[!is.na(sexpartner$sexActs),],
>>>>>> debug=TRUE)
>>>>> platform: windows 32
>>>>> executable name: glmmadmb.exe
>>>>> bin_loc:
>>>>> c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe
>>>>>
>>>>>
>>>>> using temp directory
>>>>> C:\Users\rdboylan\AppData\Local\Temp\Rtmpy2JsMY\glmmADMB17085c19f4d
>>>>> creating temp directory
>>>>> changed working directory to
>>>>> C:/Users/rdboylan/AppData/Local/Temp/Rtmpy2JsMY/glmmADMB17085c19f4d
>>>>> Command line:
>>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>>>
>>>>> -maxfn 500 -maxph 5 -noinit -shess
>>>>> Error in system(cmd, intern = intern, wait = wait | intern,
>>>>> show.output.on.console = wait, :
>>>>> 'C:/Program' not found
>>>> I'm a little bit baffled here. What happens if you use save.dir to
>>>> save the input files to a temporary directory and run
>>>>
>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>>
>>>> -maxfn 500 -maxph 5 -noinit -shess
>>>>
>>>> from the command line?
>>>
>>>> r <- glmmadmb(sexActs~(1|id),
>>>> sexpartner[!is.na(sexpartner$sexActs),], debug=TRUE,
>>>> save.dir="I:/LAMOC/Ross/")
>>> platform: windows 32
>>> executable name: glmmadmb.exe
>>> bin_loc:
>>> c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe
>>> changed working directory to I:/LAMOC/Ross
>>> Command line:
>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>> -maxfn 500 -maxph 5 -noinit -shess
>>> Error in system(cmd, intern = intern, wait = wait | intern,
>>> show.output.on.console = wait, :
>>> 'C:/Program' not found
>>>
>>> Then from a Windows Command Prompt (not cygwin)
>>> H:\>
>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32
>>> /glmmadmb.exe" -maxfn 500 -maxph 5 -noinit -shess
>>> Error trying to open data input file glmmadmb.dat
>>> Error trying to read in model data
>>> This is usual caused by a missing DAT file
>>> H:\>I:
>>>
>>> I:\>cd LAMOC/Ross
>>>
>>> I:\LAMOC\Ross>
>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/b
>>> /windows32/glmmadmb.exe" -maxfn 500 -maxph 5 -noinit -shess
>>>
>>> Initial statistics: 1 variables; iteration 0; function evaluation 0;
>>> phase 1
>>> Function value 8.9003579e+04; maximum gradient component mag
>>> -5.3609e+02
>>> Var Value Gradient |Var Value Gradient |Var Value
>>> Gradient
>>>
>>> 1 0.00000 -5.36095e+02 |
>>>
>>> - final statistics:
>>> 1 variables; iteration 8; function evaluation 14
>>> Function value 5.0576e+04; maximum gradient component mag 6.2393e-08
>>> Exit code = 1; converg criter 1.0000e-04
>>> Var Value Gradient |Var Value Gradient |Var Value
>>> Gradient
>>>
>>> 1 112.5711 6.23928e-08 |
>>> etc
>>>
>>> So it seems to work, provided I start in the save.dir. Note that is
>>> the directory that R is running in.
>>> glmmadmb.exe is still running as I hit send.
>>>
>>>>
>>>> What is the result of .Platform (and .Platform$OS in particular)
>>>
>>>> .Platform
>>> $OS.type
>>> [1] "windows"
>>>
>>> $file.sep
>>> [1] "/"
>>>
>>> $dynlib.ext
>>> [1] ".dll"
>>>
>>> $GUI
>>> [1] "RTerm"
>>>
>>> $endian
>>> [1] "little"
>>>
>>> $pkgType
>>> [1] "win.binary"
>>>
>>> $path.sep
>>> [1] ";"
>>>
>>> $r_arch
>>> [1] "i386"
>>>
>>> Ross
>>>
>>>>
>>>> It looks conceivably like R is misdiagnosing your system as *not*
>>>> being
>>>> windows, as that's the only way system() should be running. Are you
>>>> running under Cygwin (you say it's installed below) ... ?
>>>>
>>>>> changed working directory to i:/LAMOC/Ross
>>>>> removed temp directory
>>>>> C:\Users\rdboylan\AppData\Local\Temp\Rtmpy2JsMY\glmmADMB17085c19f4d
>>>>>> glmmADMB:::get_bin_loc()
>>>>> $bin_loc
>>>>> [1]
>>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>>>
>>>>>
>>>>>
>>>>> $platform
>>>>> [1] "windows"
>>>>>
>>>>> P.S. about lme4; I don't have a build environment and so trying the
>>>>> github version will not be my first move.
>>>>> Although perhaps lack of a build environment is why the second version
>>>>> is failing. I do have cygwin installed, althoughI would not expect
>>>>> R to
>>>>> know how to find it.
>>>> lme4 should be installable from lme4.r-forge.r-project.org/repos now,
>>>> as stated in a message earlier today.
>>>>
>>>
>>
>>
>>
>>------------------------------
>>
>>Message: 3
>>Date: Thu, 28 Mar 2013 15:39:07 -0700
>>From: Ross Boylan <ross at biostat.ucsf.edu>
>>To: r-sig-mixed-models at r-project.org
>>Subject: [R-sig-ME] glmmADMB truncated distributions
>>Message-ID: <5154C68B.5090807 at biostat.ucsf.edu>
>>Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>
>>How are they truncated? I expected there might be a way to set the
>>truncation threshhold, but I don't see it in the docs.
>>Ross Boylan
>>
>>
>>
>>------------------------------
>>
>>Message: 4
>>Date: Fri, 29 Mar 2013 03:39:53 +0000 (UTC)
>>From: Ben Bolker <bbolker at gmail.com>
>>To: r-sig-mixed-models at r-project.org
>>Subject: Re: [R-sig-ME] glmmADMB truncated distributions
>>Message-ID: <loom.20130329T043015-425 at post.gmane.org>
>>Content-Type: text/plain; charset=us-ascii
>>
>>Ross Boylan <ross at ...> writes:
>>
>>>
>>> How are they truncated? I expected there might be a way to set the
>>> truncation threshhold, but I don't see it in the docs.
>>> Ross Boylan
>>>
>>
>> They're truncated at zero.
>>
>> One of the problems with developing a flexible system is that everyone
>>comes up with new uses for it: the original intent of the truncated
>>distributions was just to make hurdle models possible.
>>
>> It would be possible but not completely trivial to adapt glmmADMB
>>to handle other truncation points -- you would need a straightforward
>>way to compute the cumulative distribution function of the response up
>>to the truncation point. Probably the best thing for people who want
>>to extend glmmADMB in these various ways is to learn how to use the
>>full version of AD Model Builder ...
>>
>> Ben Bolker
>>
>>
>>
>>------------------------------
>>
>>Message: 5
>>Date: Thu, 28 Mar 2013 23:19:24 -0700 (PDT)
>>From: Iasonas Lamprianou <lamprianou at yahoo.com>
>>To: "r-sig-mixed-models at r-project.org"
>> <r-sig-mixed-models at r-project.org>
>>Subject: [R-sig-ME] confidence interval for the SD of random effects
>>Message-ID:
>> <1364537964.28361.YahooMailNeo at web160104.mail.bf1.yahoo.com>
>>Content-Type: text/plain
>>
>>
>>
>>Dear all, I am running a very simple mixed effects model using lmer. I get the following results:
>>
>>Random effects: Groups Name Variance Std.Dev. id (Intercept) 7.4223 2.7244 Residual 1.9767 1.4060
>>Number of obs: 8080, groups: id, 2020
>>
>>
>>Apparently, a very good question is whethe0r the estimate of the parameters of the random effect are indeed bigger than zero. So, I used the bootstrap method described by Faraway where the distribution of likelihood rations is built between a model with and withot the random effect. I get a probbility p<0.0001, so the random effect estimates are significant. I also confirmed these findings using the library("RLRsim").
>>
>>
>>So far so good, but one of the reviewers asked me to give 95% confidence intervals. I spent three hours on the internet searching for good solutions and decided that I should use the following:
>>
>>sm1 <- mcmcsamp(mlm.1.1g, n=1000, saveb = TRUE)
>>lme4::HPDinterval(sm1)
>>
>>
>>However, the 95% lower and upper bounds do not indclude the estimate!!!! See below:
>>
>>> r$sigma lower upper
>>[1,] 1.977582 2.058084
>>attr(,"Probability")
>>[1] 0.95
>>
>>How is this possible?
>>
>>Also, I assume that the ST shows the ratio of [random effect variance]/[total variance]. In this case, my ICC is 7.42/(7.42+1.97)=0.79 but see the foloowing:
>>
>>> r$ST lower upper
>>[1,] 0.6905202 0.7342866
>>attr(,"Probability")
>>[1] 0.95
>>My empirical ICC is not included in the confidence interval. Apparently, I am missing something but after spending a lot of time, I need to ask you offer me some of your experience.
>>Thank you for your time
>>
>>Jason
>>
>>
>>
>>
>>Dr. Iasonas Lamprianou
>>Department of Social and Political Sciences
>>University of Cyprus
>>
>>
>>
>>
>>>
>> [[alternative HTML version deleted]]
>>
>>
>>
>>------------------------------
>>
>>Message: 6
>>Date: Thu, 28 Mar 2013 23:23:05 -0700 (PDT)
>>From: Iasonas Lamprianou <lamprianou at yahoo.com>
>>To: "r-sig-mixed-models at r-project.org"
>> <r-sig-mixed-models at r-project.org>
>>Subject: [R-sig-ME] variance of fixed effects
>>Message-ID:
>> <1364538185.83291.YahooMailNeo at web160102.mail.bf1.yahoo.com>
>>Content-Type: text/plain
>>
>>
>>Note: I am sending this again because I got an email saying that it was not delivered. In case it was delivered initially, I apologise for the double-posting.
>>
>>
>>
>>Dear all, I am running a very simple mixed effects model using lmer. I get the following results:
>>
>>Random
>>effects: Groups Name Variance Std.Dev. id (Intercept)
>>7.4223 2.7244 Residual 1.9767 1.4060
>>Number of obs: 8080, groups: id, 2020
>>
>>
>>Apparently,
>>a very good question is whethe0r the estimate of the parameters of the
>>random effect are indeed bigger than zero. So, I used the bootstrap
>>method described by Faraway where the distribution of likelihood rations
>>is built between a model with and withot the random effect. I get a
>>probbility p<0.0001, so the random effect estimates are significant. I
>>also confirmed these findings using the library("RLRsim").
>>
>>
>>So
>>far so good, but one of the reviewers asked me to give 95% confidence
>>intervals. I spent three hours on the internet searching for good
>>solutions and decided that I should use the following:
>>
>>sm1 <- mcmcsamp(mlm.1.1g, n=1000, saveb = TRUE)
>>lme4::HPDinterval(sm1)
>>
>>
>>However, the 95% lower and upper bounds do not indclude the estimate!!!! See below:
>>
>>> r$sigma lower upper
>>[1,] 1.977582 2.058084
>>attr(,"Probability")
>>[1] 0.95
>>
>>How is this possible?
>>
>>Also,
>>I assume that the ST shows the ratio of [random effect variance]/[total
>>variance]. In this case, my ICC is 7.42/(7.42+1.97)=0.79 but see the
>>foloowing:
>>
>>> r$ST lower upper
>>[1,] 0.6905202 0.7342866
>>attr(,"Probability")
>>[1] 0.95
>>My
>>empirical ICC is not included in the confidence interval. Apparently, I
>>am missing something but after spending a lot of time, I need to ask
>>you offer me some of your experience.
>>Thank you for your time
>>
>>Jason
>>
>>
>>
>>Dr. Iasonas Lamprianou
>>Department of Social and Political Sciences
>>University of Cyprus
>> [[alternative HTML version deleted]]
>>
>>
>>
>>------------------------------
>>
>>_______________________________________________
>>R-sig-mixed-models mailing list
>>R-sig-mixed-models at r-project.org
>>https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>End of R-sig-mixed-models Digest, Vol 75, Issue 45
>>**************************************************
>>
>>
>>
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://joshuawiley.com/
Senior Analyst - Elkhart Group Ltd.
http://elkhartgroup.com
More information about the R-sig-mixed-models
mailing list