John Sorkin | 26 Jun 05:00 2016

function for over dispersed poisson regression in the setting of a random effects model

Is there a function that will run a model appropriate for over dispersed data (such as a negative binomial or quasipoisson)
with a random effects (or mixed effects) model in R? GLMER will not accept:  
family=quasipoisson(link="log") or
family=negbinomial(link="log") 

I want to run something like the following:
fit0 <- glmer(Fall ~ Group+(1|PID)+offset(log(TimeYrs)),family=quasipoisson(link="log"),data=data)
Thank  you
John
John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric Medicine
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing) 

Confidentiality Statement:
This email message, including any attachments, is for the sole use of the intended recipient(s) and may
contain confidential and privileged information. Any unauthorized use, disclosure or distribution is
prohibited. If you are not the intended recipient, please contact the sender by reply email and destroy
all copies of the original message. 
Gonçalo Ferraz | 25 Jun 16:13 2016
Picon

strange behavior of lchoose in combinatorics problem

Hi,

I am working on interactions between animals, studying whether animal 1 is attracted to animal 2 (or
vice-versa). I looked for the two animals in N=2178 sampling occasions, finding animal 1 a total of N1=165
times, and animal 2 a total of N2=331 times. In J=97 occasions, I saw both animals at the same time. 

The more frequently I see the two animals in the same sampling occasion, the more I will believe that they are
attracted to each other. Therefore, I want to calculate the probability of finding J<=97 when the two
animals are moving around independently of each other. The higher this probability, the stronger the attraction.

Following Veech (Journal of Biogeography 2014, 41: 1029-1035) I compute the log probability of obtaining
a number n of encounters between animals as ‘lpn’ in the function lveech:

lveech <- function(N,N1,N2,n) {
    a <- lchoose(N,n)
    b <- lchoose(N-n,N2-n)
    c <- lchoose(N-N2,N1-n)
    d <- lchoose(N,N2)
    e <- lchoose(N,N1)
    lpn <- (a+b+c)-(d+e)
    return(lpn)
}

I have tested this function with shorter, intuitive numbers, and it works as expected. I use log
probabilities to stay away from intractable numbers.

Next, I use function lveech to obtain a vector ‘lpvec’ that gives me all the log probabilities of
getting n=0,1,2,…,97 simultaneous observations of the two animals:

lpvec <- rep(NA,J+1)
(Continue reading)

Shivi Bhatia | 25 Jun 21:05 2016
Picon

Help with the Cut Function

Dear Team,

Please see the code below:

Age1<- cut(desc$Age, breaks = c(20,30,40,Inf),labels = c("Low","Mid","Top"))
here i am creating three categories as mentioned from the age var from desc
data set.
All the values are set correctly however the values which are below 20 are
set to NA.
Is there anything i am doing incorrect.

Regards, Shivi

	[[alternative HTML version deleted]]

T.Riedle | 25 Jun 15:52 2016
Picon

Error using function seas()

Dear all,

I am trying to run the seas() function. If I run the seas() function as shown below I get following errors.
What is wrong with my code?

> data<-Shiller_data[,2]

> ts<-ts(data,start=c(1871, 01), end=c(2015, 12), frequency=12)

> ts

         Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep     Oct

1871    4.44    4.50    4.61    4.74    4.86    4.82    4.73    4.79    4.84    4.59

1872    4.86    4.88    5.04    5.18    5.18    5.13    5.10    5.04    4.95    4.97

1873    5.11    5.15    5.11    5.04    5.05    4.98    4.97    4.97    4.59    4.19

1874    4.66    4.80    4.73    4.60    4.48    4.46    4.46    4.47    4.54    4.53

1875    4.54    4.53    4.59    4.65    4.47    4.38    4.39    4.41    4.37    4.30

1876    4.46    4.52    4.51    4.34    4.18    4.15    4.10    3.93    3.69    3.67

1877    3.55    3.34    3.17    2.94    2.94    2.73    2.85    3.05    3.24    3.31

1878    3.25    3.18    3.24    3.33    3.34    3.41    3.48    3.45    3.52    3.48

1879    3.58    3.71    3.65    3.77    3.94    3.96    4.04    4.07    4.22    4.68
(Continue reading)

Marius Hofert | 25 Jun 02:40 2016
Picon
Picon

Re: What's box() (exactly) doing?

Hi Jim,

Thanks a lot, exactly what I was looking for.

Cheers,
Marius

On Thu, Jun 23, 2016 at 11:06 PM, Jim Lemon <drjimlemon <at> gmail.com> wrote:
> Hi Marius,
> There are a few things that are happening here. First, the plot area
> is not going to be the same as your x and y limits unless you say so:
>
> # run your first example
> par("usr")
> [1] -0.04  1.04 -0.04  1.04
>
> # but
> plot(NA, type = "n", ann = FALSE, axes = FALSE,
>  xlim = 0:1, ylim = 0:1,xaxs="i",yaxs="i")
> box()
> rect(xleft = 0, ybottom = 0, xright = 1, ytop = 1, col = "grey80")
> par("usr")
> [1] 0 1 0 1
>
> Second, the "rect" function is automatically clipped to the plot area,
> so you may lose a bit at the edges if you don't override this:
>
> par(xpd=TRUE)
> rect(...)
> par(xpd=FALSE)
(Continue reading)

Steven Yen | 24 Jun 22:55 2016
Picon

Reading csv file with missing value

I read a csv file (with read.csv) containing missing values (as shown 
below). Is there a convenient way to set these NA into zeros?
Better yet, is there an option to assign zeros to these blank cells in 
reading the csv file? Thank you!

NA  -1  NA  NA  NA   1  NA
NA  NA  NA  NA  NA  NA  NA
NA  NA  NA  NA  NA  NA  NA
NA  -1  NA  NA  NA  NA  NA

	[[alternative HTML version deleted]]

deva d | 24 Jun 20:51 2016
Picon

R package updating

hi all,

i notice that the R package available at CRAN is a more recent one compared
to what I have.

but i failed to update my machine and unfortunately, R does not work any
more.

can someone kindly suggest what is a way to update the R ver.

is there a good way to retain all R files and update it without upsetting
other settings.

would be happy to hear.

regds,

*....*

*D*


...............

*in search of knowledge, everyday something is added ....*

*in search of wisdom, everyday something is dropped  ... an old Chinese
Proverb*
:::::::::::::::::::::::::

(Continue reading)

deva d | 24 Jun 20:49 2016
Picon

Rattle

hi all,

i am getting stuck while using Rattle - specifically, I am unable to get
the graphics.

i am using R 3.2.5 and some packages do not work in that.

can someone pl help ?

*....*

*Deva*


...............

*in search of knowledge, everyday something is added ....*

*in search of wisdom, everyday something is dropped  ... an old Chinese
Proverb*
:::::::::::::::::::::::::

On Fri, Jun 24, 2016 at 5:34 PM, Bert Gunter <bgunter.4567 <at> gmail.com> wrote:

> This is very basic. Have you gone through any R tutorials? There are
> many good ones on the web. e.g., see here:
> https://www.rstudio.com/online-learning/#R
>
> In any case, you should not expect this list to teach you basic R. You
> *should* expect it to help you learn and improve your own efforts. I
(Continue reading)

fgoetz | 24 Jun 19:45 2016

Fwd: Fwd: RE: Heatmap.2 Breaks argument

Dear Mr. or Mrs.,

Please see the previous messages for information. I had a problem with 
the heatmap.2 breaks argument and was wondering if someone could help me 
with that problem. I could not find information for the other authors of 
heatmap.2 to contact them.
Could you provide help or contact information?

Best regards,

Florian Goetz

Phd candidate

Woods Hole Oceanographic Institution

508 289 3715

-------- Weitergeleitete Nachricht --------
Betreff: 	RE: Heatmap.2 Breaks argument
Datum: 	Fri, 24 Jun 2016 13:16:53 -0400
Von: 	Liaw, Andy <andy_liaw <at> merck.com>
An: 	fgoetz <fgoetz <at> whoi.edu>

Hi Florian,

I wrote the original version of heatmap(), but not heatmap.2().  Also, heatmap() was modified by R Core
before being included into base R, and is now supported by R Core.  Please contact the author of heatmap.2()
first.  Thanks!

(Continue reading)

André Luis Neves | 24 Jun 18:05 2016
Picon
Picon

How to do it in R

Dear all,

I`ve got to calculate the ratio of methanogens to bacteria, but I wouldn`t
like to divide the total copy numbers of methanogens ( on average 10^8) by
bacteria (10^10) because they have different exponents and bases. So, my
idea is to standardize both microorganisms counts to 10^3.

Hypothetical example of what I`d like to do:

Total Methanogens: 28 x 10^3
Total bacteria: 710 x 10^3

Total: 710 + 28= 738 x 10^3

Met/bac Ratio : (28/738)*100 = 3.79%

How could I perform this calculation in R?

Thanks,

--

-- 
Andre

	[[alternative HTML version deleted]]

John Sorkin | 24 Jun 15:45 2016

Add column to the output of summary(glht).


I am trying to make the leap from an R users to an R aficionado . . .

I am trying to understand how add a column to the output of summary (and to understand how summary() works).

I have run a glmer
fit0 <- glmer(Fall ~ Group+(1|PID),family=poisson(link="log"),data=data[data[,"Group"]!=0,])

and I want to perform adjusted multiple comparisons:

SumTukey <- summary(glht(fit0, linfct= mcp(Group="Tukey")))

which gives beautiful output:

> SumTukey

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glmer(formula = Fall ~ Group + (1 | PID), data = data[data[, 
    "Group"] != 0, ], family = poisson(link = "log"))

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0   0.5320     0.5075   1.048    0.717
3 - 1 == 0   0.6554     0.5000   1.311    0.551
4 - 1 == 0   0.9357     0.4655   2.010    0.181
3 - 2 == 0   0.1234     0.4174   0.296    0.991
4 - 2 == 0   0.4037     0.3754   1.075    0.700
(Continue reading)


Gmane