[R] list contingency tables

[li li]

Hi all,
  I am trying to list all the 4 by 2 tables with some fixed margins.
  For example, consider 4 by 2 tables with row margins 1,2,2,1 and
column margins 3,3. I was able to do it using the code below. However,
as seen below, I had to first count the total number of tables with
the specific row margins and column margins in order to create space
to store the tables.
Is there a way to skip the step of counting the number of tables?
  Also, wanted to avoid for loops as much as possible since it can be
extremely slow and inefficient.
   Thanks so much in advance for you insight and help.
   Hanna



> library(gtools)
> A <- permutations(n=4,r=2,v=0:3, repeats.allowed=TRUE)
> B <- apply(A, 1, sum)
> rmg <- c(1,2,2,1)
> cmg <- c(3,3)
> m1 <- t(A[which(B==1),])
> m2 <- t(A[which(B==2),])
> m3 <- t(A[which(B==2),])
>
> ##count number of tables with row margins 1,2,2,1 and column margins 3,3.
> num <- 0
> for (i in 1:ncol(m1)){
+ for (j in 1:ncol(m2)){
+ for (k in 1:ncol(m3)){
+ M <- t(cbind(m1[,i], m2[,j], m3[,k]))
+ M1 <- rbind(M, cmg-apply(M,2,sum))
+ num <- num+(sum(M1[4,] < 0) == 0)
+ }}}
>
>
> #create space to store the tables
> C <- array(NA, dim=c(4,2,num))
>
> # list all the tables with fixed margins
> num <- 0
> for (i in 1:ncol(m1)){
+ for (j in 1:ncol(m2)){
+ for (k in 1:ncol(m3)){
+ M <- t(cbind(m1[,i], m2[,j], m3[,k]))
+ M1 <- rbind(M,cmg-apply(M,2,sum))
+ if (sum(M1[4,] < 0) == 0) {
+ num <- num+1
+C[,,num] <- M1
+ }
+ }}}
>
> C
, , 1

 [,1] [,2]
[1,]01
[2,]02
[3,]20
[4,]10

, , 2

 [,1] [,2]
[1,]01
[2,]11
[3,]11
[4,]10

, , 3

 [,1] [,2]
[1,]01
[2,]11
[3,]20
[4,]01

, , 4

 [,1] [,2]
[1,]01
[2,]20
[3,]02
[4,]10

, , 5

 [,1] [,2]
[1,]01
[2,]20
[3,]11
[4,]01

, , 6

 [,1] [,2]
[1,]10
[2,]02
[3,]11
[4,]10

, , 7

 [,1] [,2]
[1,]10
[2,]02
[3,]20
[4,]01

, , 8

 [,1] [,2]
[1,]10
[2,]11
[3,]02
[4,]10

, , 9

 [,1] [,2]
[1,]10
[2,]11
[3,]11
[4,]01

, , 10

 [,1] [,2]
[1,]10
[2,]20
[3,]02
[4,]01

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Identify row indices corresponding to each distinct row of a matrix

[li li]

Thanks. It makes sense.

Jeff Newmiller  于2018年11月8日周四 下午8:05写道：

> The duplicated function returns TRUE for rows that have already
> appeared... exactly one of the rows is not represented in the output of
> duplicated. For the intended purpose of removing duplicates this behavior
> is ideal. I have no idea what your intended purpose is, since every row has
> duplicates elsewhere in the matrix. If you really want every set identified
> this way then a loop/apply seems inevitable (most opportunities for
> optimization come about by not visiting every combination).
>
> Cm <- as.matrix( C )
> D <- which( !duplicated( Cm, MARGIN=1 ) )
> nCm <- nrow( Cm )
> F <- lapply( D, function(d) {
>idxrep <- rep( d, nCm )
>which( 0 == unname( rowSums( Cm[idxrep,] != Cm ) ) )
>   } )
>
>
> On November 8, 2018 1:42:40 PM PST, li li  wrote:
> >Thanks to all the reply. I will try to use plain text in the future.
> >One question regarding using "which( ! duplicated( m, MARGIN=1 ) )".
> >This seems to return the fist row indices corresponding to the distinct
> >rows but it does not give all the row indices
> >corresponding to each of the distinct rows. For example, in the my
> >example
> >below, rows 1, 13 15 are all (1,9).
> >Thanks.
> >  Hanna
> >> A <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),8,2)
> >> B <- rbind(A,A,A)
> >> C <- as.data.frame(B[sample(nrow(B)),])
> >> C
> >   V1 V2
> >1   1  9
> >2   2 10
> >3   3 11
> >4   5 13
> >5   7 15
> >6   6 14
> >7   4 12
> >8   3 11
> >9   8 16
> >10  5 13
> >11  7 15
> >12  2 10
> >13  1  9
> >14  8 16
> >15  1  9
> >16  3 11
> >17  7 15
> >18  4 12
> >19  2 10
> >20  6 14
> >21  4 12
> >22  8 16
> >23  5 13
> >24  6 14
> >> T <- unique(C)
> >> T
> >  V1 V2
> >1  1  9
> >2  2 10
> >3  3 11
> >4  5 13
> >5  7 15
> >6  6 14
> >7  4 12
> >9  8 16
> >>
> >> i <- 1
> >> which(C[,1]==T[i,1]& C[,2]==T[i,2])
> >[1]  1 13 15
> >
> >
> >Bert Gunter  于2018年11月8日周四 上午10:43写道：
> >
> >> Yes -- much better than mine. I didn't know about the MARGIN argument
> >of
> >> duplicated().
> >>
> >> -- Bert
> >>
> >>
> >> On Wed, Nov 7, 2018 at 10:32 PM Jeff Newmiller
> >
> >> wrote:
> >>
> >>> Perhaps
> >>>
> >>> which( ! duplicated( m, MARGIN=1 ) )
> >>>
> >>> ? (untested)
> >>>
> >>> On November 7, 2018 9:20:57 PM PST, Bert Gunter
> >
> >>> wrote:
> >>> >A mess -- due to your continued use of html formatting.
> >>> >
> >>> >But something like this may do what you want (hard to tell with the
> >>> >mess):
> >>> >
> >>> >> m <- matrix(1:16,nrow=8)[rep(1:8,2),]
> >>> >> m
> >>> >  [,1] [,2]
> >>> > [1,]19
> >>> > [2,]2   10
> >>> > [3,]3   11
> >>> > [4,]4   12
> >>> > [5,]5   13
> >>> > [6,]6   14
> >>> > [7,]7   15
> >>> > [8,]8   16
> >>> > [9,]19
> >>> >[10,]2   10
> >>> >[11,]3   11
> >>> >[12,]4   12
> >>> >[13,]5   13
> >>> >[14,]6   14
> >>> >[15,]7   15
> >>> >[16,]8   16
> >>> >> vec <- apply(m,1,paste,collapse="-") ## converts rows into
> >character
> >>> >vector
> >>> >> vec
> >>> >[1] "1-9"  "2-10" "3-11" "4-12" "5-13" "6-14" "7-15" "8-16" "1-9"
> >>> >"2-10"
> >>> >"3-11" "4-12" "5-13" "6-14"
> >>> >[15] "7-15" "8-16"
> >>> >> ## Then maybe:
> >>> >> tapply(seq_along(vec),vec, I)
> >>> >$`1-9`
> >>> >[1] 1 9
> >>> >
> >>> >$`2-10`
> >>> >[1]  2 10
> >>> >
> >>> >$`3-11`
> >>> >[1]  3 11
> >>> >
> >>> >$`4-12`
> >>> >[1]  4 12
> >>> >
> >>> >$`5-13`
> >>> >

Re: [R] Identify row indices corresponding to each distinct row of a matrix

[li li]

Thanks to all the reply. I will try to use plain text in the future.
One question regarding using "which( ! duplicated( m, MARGIN=1 ) )".
This seems to return the fist row indices corresponding to the distinct
rows but it does not give all the row indices
corresponding to each of the distinct rows. For example, in the my example
below, rows 1, 13 15 are all (1,9).
Thanks.
  Hanna
> A <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),8,2)
> B <- rbind(A,A,A)
> C <- as.data.frame(B[sample(nrow(B)),])
> C
   V1 V2
1   1  9
2   2 10
3   3 11
4   5 13
5   7 15
6   6 14
7   4 12
8   3 11
9   8 16
10  5 13
11  7 15
12  2 10
13  1  9
14  8 16
15  1  9
16  3 11
17  7 15
18  4 12
19  2 10
20  6 14
21  4 12
22  8 16
23  5 13
24  6 14
> T <- unique(C)
> T
  V1 V2
1  1  9
2  2 10
3  3 11
4  5 13
5  7 15
6  6 14
7  4 12
9  8 16
>
> i <- 1
> which(C[,1]==T[i,1]& C[,2]==T[i,2])
[1]  1 13 15


Bert Gunter  于2018年11月8日周四 上午10:43写道：

> Yes -- much better than mine. I didn't know about the MARGIN argument of
> duplicated().
>
> -- Bert
>
>
> On Wed, Nov 7, 2018 at 10:32 PM Jeff Newmiller 
> wrote:
>
>> Perhaps
>>
>> which( ! duplicated( m, MARGIN=1 ) )
>>
>> ? (untested)
>>
>> On November 7, 2018 9:20:57 PM PST, Bert Gunter 
>> wrote:
>> >A mess -- due to your continued use of html formatting.
>> >
>> >But something like this may do what you want (hard to tell with the
>> >mess):
>> >
>> >> m <- matrix(1:16,nrow=8)[rep(1:8,2),]
>> >> m
>> >  [,1] [,2]
>> > [1,]19
>> > [2,]2   10
>> > [3,]3   11
>> > [4,]4   12
>> > [5,]5   13
>> > [6,]6   14
>> > [7,]7   15
>> > [8,]8   16
>> > [9,]19
>> >[10,]2   10
>> >[11,]3   11
>> >[12,]4   12
>> >[13,]5   13
>> >[14,]6   14
>> >[15,]7   15
>> >[16,]8   16
>> >> vec <- apply(m,1,paste,collapse="-") ## converts rows into character
>> >vector
>> >> vec
>> >[1] "1-9"  "2-10" "3-11" "4-12" "5-13" "6-14" "7-15" "8-16" "1-9"
>> >"2-10"
>> >"3-11" "4-12" "5-13" "6-14"
>> >[15] "7-15" "8-16"
>> >> ## Then maybe:
>> >> tapply(seq_along(vec),vec, I)
>> >$`1-9`
>> >[1] 1 9
>> >
>> >$`2-10`
>> >[1]  2 10
>> >
>> >$`3-11`
>> >[1]  3 11
>> >
>> >$`4-12`
>> >[1]  4 12
>> >
>> >$`5-13`
>> >[1]  5 13
>> >
>> >$`6-14`
>> >[1]  6 14
>> >
>> >$`7-15`
>> >[1]  7 15
>> >
>> >$`8-16`
>> >[1]  8 16
>> >
>> >> ## gives the row numbers for each unique row
>> >
>> >There may well be slicker ways to do this -- if this is actually what
>> >you
>> >want to do.
>> >
>> >-- Bert
>> >
>> >
>> >
>> >On Wed, Nov 7, 2018 at 7:56 PM li li  wrote:
>> >
>> >> Hi all,
>> >>I use the following example to illustrate my question. As you can
>> >see,
>> >> in matrix C some rows are repeated and I would like to find the
>> >indices of
>> >> the rows corresponding to each of the distinct rows.
>> >>   For example, for the row c(1,9), I have used the "which" function
>> >to
>> >> identify the row indices corresponding to c(1,9). Using this
>> >approach, in
>> >> order to cover all distinct rows, I need to use a for loop.
>> >>I am wondering whether there is an easier way where a for loop can
>> >be
>> >> avoided?
>> >>Thanks very much!
>> >>   Hanna
>> >>
>> >>
>> >>
>> >> > A <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),8,2)> B <-
>> >> rbind(A,A,A)> C <- as.data.frame(B[sample(nrow(B)),])> C   V1 V2
>> >> 1   1  9
>> >> 2   2 10
>> >> 3   3 11
>> >> 4   5 13
>> >> 5   7 15
>> >> 6   6 14
>> >> 7   4 12
>> >> 8   3 11
>> >> 9   8 16
>> >> 10  5 13
>> >> 11  7 15
>> >> 12  2 10
>> >> 13  1  9
>> >> 14  8 16
>> >> 15  1  9
>> >> 16  3 11
>&g

[R] Identify row indices corresponding to each distinct row of a matrix

[li li]

Hi all,
   I use the following example to illustrate my question. As you can see,
in matrix C some rows are repeated and I would like to find the indices of
the rows corresponding to each of the distinct rows.
  For example, for the row c(1,9), I have used the "which" function to
identify the row indices corresponding to c(1,9). Using this approach, in
order to cover all distinct rows, I need to use a for loop.
   I am wondering whether there is an easier way where a for loop can be
avoided?
   Thanks very much!
  Hanna



> A <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),8,2)> B <- 
> rbind(A,A,A)> C <- as.data.frame(B[sample(nrow(B)),])> C   V1 V2
1   1  9
2   2 10
3   3 11
4   5 13
5   7 15
6   6 14
7   4 12
8   3 11
9   8 16
10  5 13
11  7 15
12  2 10
13  1  9
14  8 16
15  1  9
16  3 11
17  7 15
18  4 12
19  2 10
20  6 14
21  4 12
22  8 16
23  5 13
24  6 14> T <- unique(C)> T  V1 V2
1  1  9
2  2 10
3  3 11
4  5 13
5  7 15
6  6 14
7  4 12
9  8 16> > i <- 1> which(C[,1]==T[i,1]&
C[,2]==T[i,2])[1]  1 13 15

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] question with integrate function

[li li]

Oh ok. Thanks very much. I will have to restrict to a shorter interval.
Hanna

2018-02-06 14:33 GMT-05:00 Göran Broström :

> Hi Hanna,
>
> your function is essentially zero outside a short interval around 9. And
> the help page states: "If the function is approximately constant (in
> particular, zero) over nearly all its range it is possible that the result
> and error estimate may be seriously wrong."
>
> You could try to integrate over a finite interval, say (7, 12).
>
> Göran Broström
>
>
> On 2018-02-06 19:40, li li wrote:
>
>> Sorry. I meant in the previous email that the function h() is a monotone
>> decreasing function. Thanks very much.
>>
>> 2018-02-06 13:32 GMT-05:00 li li :
>>
>> Hi all,
>>>The function h below is a function of c and it should be a monotone
>>> increasing function since the integrand is nonnegative and integral is
>>> taken from c to infinity. However, as we can see from the plot, it is not
>>> shown to be monotone. Something wrong with the usage of integrate
>>> function?
>>> Thanks so much for your help.
>>>  Hanna
>>>
>>>
>>>
>>> h <- function(c){
>>>  g <- function(x){pnorm(x-8.8, mean=0.4, sd=0.3,
>>> lower.tail=TRUE)*dnorm(x, mean=9,sd=0.15)}
>>>  integrate(g, lower=c, upper=Inf)$value}
>>>
>>> xx <- seq(-20,20,by=0.001)
>>> y <- xx
>>> for (i in 1:length(xx)){y[i] <- h(xx[i])}
>>> plot(xx, y)
>>>
>>>
>>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] question with integrate function

[li li]

Sorry. I meant in the previous email that the function h() is a monotone
decreasing function. Thanks very much.

2018-02-06 13:32 GMT-05:00 li li :

> Hi all,
>   The function h below is a function of c and it should be a monotone
> increasing function since the integrand is nonnegative and integral is
> taken from c to infinity. However, as we can see from the plot, it is not
> shown to be monotone. Something wrong with the usage of integrate function?
> Thanks so much for your help.
> Hanna
>
>
>
> h <- function(c){
> g <- function(x){pnorm(x-8.8, mean=0.4, sd=0.3,
> lower.tail=TRUE)*dnorm(x, mean=9,sd=0.15)}
> integrate(g, lower=c, upper=Inf)$value}
>
> xx <- seq(-20,20,by=0.001)
> y <- xx
> for (i in 1:length(xx)){y[i] <- h(xx[i])}
> plot(xx, y)
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] question with integrate function

[li li]

Hi all,
  The function h below is a function of c and it should be a monotone
increasing function since the integrand is nonnegative and integral is
taken from c to infinity. However, as we can see from the plot, it is not
shown to be monotone. Something wrong with the usage of integrate function?
Thanks so much for your help.
Hanna



h <- function(c){
g <- function(x){pnorm(x-8.8, mean=0.4, sd=0.3,
lower.tail=TRUE)*dnorm(x, mean=9,sd=0.15)}
integrate(g, lower=c, upper=Inf)$value}

xx <- seq(-20,20,by=0.001)
y <- xx
for (i in 1:length(xx)){y[i] <- h(xx[i])}
plot(xx, y)

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] An iterative function

[li li]

Dear all,
  The function f() below is a function of m1 and m2, both of which are
matrices with 3 rows. The function works sequentially one row after
another.
So altogether there are three stages. I am trying to update the coding to
write a generic function that will work for arbitrary k stages.
  I am hoping to get some suggestion and help. Thanks so much!
 Hanna




 ##x, y are two vectors of the same length
f0 <- function(x, y){
  I <- which((x-y)^2>0.5)
  if (length(I)==0){
a <- 0; b <-0; c<- 0
  } else {
  a <- min(I)
  b <- x[a]
  c <- y[a]}
  return(list(a=a, b=b, c=c))
}

##both m1 and m2 are matrix with 3 rows and same number of columns
f <- function(m1, m2){
  n <- dim(m1)[2]
  tmp1 <- f0(m1[1,], m2[1,])
  S2 <- which(m1[1,] > tmp1$a)
  if (length(S2) == 0){
t1 <- c(tmp1$b, 0, 0)
t2 <- c(tmp1$c, 0, 0)} else {
tmp2 <- f0(m1[2,S2], m2[2, (n-length(S2)+1):n])
S3 <- S2[which(m2[2, S2] > tmp2$a)]
if (length(S3) == 0) {
  t1 <- c(tmp1$b, tmp2$b, 0)
  t2 <- c(tmp1$c, tmp2$c, 0)} else {
tmp3 <- f0(m1[3,S3], m2[3,  (n-length(S3)+1):n])
t1 <- c(tmp1$b, tmp2$b, tmp3$c)
t2 <- c(tmp1$c, tmp2$c, tmp3$b)
  }}
  return(list(t1=t1, t2=t2))
}

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] problem with "unique" function

[li li]

Thank you Peter. Rounding fixed the problem.


2017-07-28 12:15 GMT-04:00 peter dalgaard :

> Most likely, previous computations have ended up giving slightly different
> values of say 0.1. A pragmatic way out is to round to, say, 5 digits
> before applying unique. In this particular case, it seems like all numbers
> are multiples of 1/30, so another idea could be to multiply by 30, round,
> and divide by 30.
>
> -pd
>
> > On 28 Jul 2017, at 17:17 , li li  wrote:
> >
> > I have the joint distribution of three discrete random variables z1, z2
> and
> > z3 which is captured by "z"
> > and "prob" as described below.
> >
> > For example, the probability for z1=0.46667, z2=-1 and z3=-1 is
> 2.752e-13.
> > Also, the probability adds up to 1.
> >
> >> head(z)   z1  z2  z3
> > [1,] -0.46667 -1. -1.
> > [2,] -0.3 -0.9333 -0.9333
> > [3,] -0.2 -0.8667 -0.8667
> > [4,] -0.06667 -0.8000 -0.8000
> > [5,]  0.06667 -0.7333 -0.7333
> > [6,]  0.2 -0.6667 -0.6667> prob[1:5][1] 2.752e-13 3.210e-12
> > 1.348e-11 2.656e-11 2.656e-11> sum(prob)[1] 1
> >
> >
> > I want to put the distribution into a joint probability table. I use the
> > following code. But the probability no longer adds up to 1.
> >
> >> z1 <- sort(unique(z[,1])); z2 <- sort(unique(z[,2]));  z3 <-
> sort(unique(z[,3]))> P <- array(0, dim=c(length(z1),length(z2),length(z3)),
> dimnames=list(A=z1, H=z2, M=z3))> > for (i in 1:(dim(z)[1])){+ ind <-
> z[i,]+ P[dimnames(P)$A==ind[1], dimnames(P)$H==ind[2],
> dimnames(P)$M==ind[3]] <- prob[i] + }> sum(P)[1] 37.5
> >
> >
> > The problem is when we look z1 as below, there are lot of repeated
> values.
> >
> >> unique(z[,1]) [1] -0.46667 -0.3 -0.2 -0.06667  0.06667
> 0.2  0.3  0.46667 -0.5 -0.4
> > [11] -0.26667 -0.1  0.0  0.1  0.26667  0.4  0.5
> > -0.6 -0.06667  0.06667
> > [21]  0.6 -0.7 -0.1  0.1  0.7 -0.7 -0.6
> > -0.2  0.2  0.6
> > [31]  0.7 -0.8 -0.5 -0.26667  0.26667  0.5  0.8
> > -0.86667 -0.46667 -0.3
> > [41]  0.3  0.46667  0.86667 -0.9 -0.4  0.4  0.9
> > -1.0 -0.46667 -0.3
> > [51]  0.3  0.46667  1.0 -0.5 -0.26667  0.26667  0.5
> > -0.2  0.2 -0.7
> > [61] -0.1  0.1  0.7 -0.7 -0.06667  0.06667  0.7
> > -0.2  0.2 -0.1
> > [71]  0.1 -0.06667  0.06667 -0.06667  0.06667
> >
> >
> >
> > Is there a way to fix this? Any idea and suggestions? Thanks very much!!
> >
> >   Hanna
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd@cbs.dk  Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] problem with "unique" function

[li li]

I have the joint distribution of three discrete random variables z1, z2 and
z3 which is captured by "z"
and "prob" as described below.

For example, the probability for z1=0.46667, z2=-1 and z3=-1 is 2.752e-13.
Also, the probability adds up to 1.

> head(z)   z1  z2  z3
[1,] -0.46667 -1. -1.
[2,] -0.3 -0.9333 -0.9333
[3,] -0.2 -0.8667 -0.8667
[4,] -0.06667 -0.8000 -0.8000
[5,]  0.06667 -0.7333 -0.7333
[6,]  0.2 -0.6667 -0.6667> prob[1:5][1] 2.752e-13 3.210e-12
1.348e-11 2.656e-11 2.656e-11> sum(prob)[1] 1


I want to put the distribution into a joint probability table. I use the
following code. But the probability no longer adds up to 1.

> z1 <- sort(unique(z[,1])); z2 <- sort(unique(z[,2]));  z3 <- 
> sort(unique(z[,3]))> P <- array(0, dim=c(length(z1),length(z2),length(z3)), 
> dimnames=list(A=z1, H=z2, M=z3))> > for (i in 1:(dim(z)[1])){+ ind <- 
> z[i,]+ P[dimnames(P)$A==ind[1], dimnames(P)$H==ind[2], 
> dimnames(P)$M==ind[3]] <- prob[i] + }> sum(P)[1] 37.5


The problem is when we look z1 as below, there are lot of repeated values.

> unique(z[,1]) [1] -0.46667 -0.3 -0.2 -0.06667  0.06667  0.2  
> 0.3  0.46667 -0.5 -0.4
[11] -0.26667 -0.1  0.0  0.1  0.26667  0.4  0.5
-0.6 -0.06667  0.06667
[21]  0.6 -0.7 -0.1  0.1  0.7 -0.7 -0.6
-0.2  0.2  0.6
[31]  0.7 -0.8 -0.5 -0.26667  0.26667  0.5  0.8
-0.86667 -0.46667 -0.3
[41]  0.3  0.46667  0.86667 -0.9 -0.4  0.4  0.9
-1.0 -0.46667 -0.3
[51]  0.3  0.46667  1.0 -0.5 -0.26667  0.26667  0.5
-0.2  0.2 -0.7
[61] -0.1  0.1  0.7 -0.7 -0.06667  0.06667  0.7
-0.2  0.2 -0.1
[71]  0.1 -0.06667  0.06667 -0.06667  0.06667



Is there a way to fix this? Any idea and suggestions? Thanks very much!!

   Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Adding zeros in each dimension of an array

[li li]

Thanks! That is the method I am using now but I was checking whether there
is a simpler option. I will look into the abind package too.
   Hanna

2017-06-07 12:25 GMT-04:00 Jeff Newmiller :

> Please read
> https://www.lifewire.com/how-to-send-a-message-in-plain-text
> -from-gmail-1171963
>
> Re your question: easy is in the eye of the beholder.
>
> a <- array( 1:12, dim = c( 2, 2, 3 ) )
> b <- array( 0, dim = dim( a ) + 1 )
> b[ 1:2, 1:2, 1:3 ] <- a
>
> If you want to do a lot of this, it could be wrapped for convenience
> though implementing that last expression without hardcoding the sequences
> may not be obvious:
>
> expandArray <- function( a ) {
>   b <- array( 0, dim = dim( a ) + 1 )
>   do.call( `[<-`, c( list( b ), lapply( dim( a ), seq.int ), list( a ) ) )
> }
> expandArray( a )
>
>
> On Wed, 7 Jun 2017, li li wrote:
>
> For a data frame, we can add an additional row or column easily. For
>> example, we can add an additional row of zero and an additional row of
>> column as follows.
>>
>> Is there an easy and similar way to add zeros in each dimension? For
>> example, for
>> array(1:12, dim=c(2,2,3))?
>>
>> Thanks for your help!!
>>   Hanna
>>
>>
>>
>> x <- as.data.frame(matrix(1:20,4,5))> x[5,] <- 0> x[,6] <- 0> x  V1 V2
>>> V3 V4 V5 V6
>>>
>> 1  1  5  9 13 17  0
>> 2  2  6 10 14 18  0
>> 3  3  7 11 15 19  0
>> 4  4  8 12 16 20  0
>> 5  0  0  0  0  0  0
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
> 
> ---
> Jeff NewmillerThe .   .  Go Live...
> DCN:Basics: ##.#.   ##.#.  Live
> Go...
>   Live:   OO#.. Dead: OO#..  Playing
> Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
> /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
> 
> ---
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Adding zeros in each dimension of an array

[li li]

For a data frame, we can add an additional row or column easily. For
example, we can add an additional row of zero and an additional row of
column as follows.

Is there an easy and similar way to add zeros in each dimension? For
example, for
array(1:12, dim=c(2,2,3))?

Thanks for your help!!
   Hanna



> x <- as.data.frame(matrix(1:20,4,5))> x[5,] <- 0> x[,6] <- 0> x  V1 V2 V3 V4 
> V5 V6
1  1  5  9 13 17  0
2  2  6 10 14 18  0
3  3  7 11 15 19  0
4  4  8 12 16 20  0
5  0  0  0  0  0  0

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] An R question

[li li]

Thank you!

2017-06-07 10:38 GMT-04:00 Ivan Calandra :

> Hi,
>
> Check the FAQ 7.31
> https://cran.rstudio.com/doc/FAQ/R-FAQ.html#Why-doesn_0027t-
> R-think-these-numbers-are-equal_003f
>
> And read the posting guide too...
> https://www.r-project.org/posting-guide.html
>
> HTH,
> Ivan
>
> --
> Dr. Ivan Calandra
> TraCEr, Laboratory for Traceology and Controlled Experiments
> MONREPOS Archaeological Research Centre and
> Museum for Human Behavioural Evolution
> Schloss Monrepos
> 56567 Neuwied, Germany
> +49 (0) 2631 9772-243
> https://www.researchgate.net/profile/Ivan_Calandra
>
>
> On 07/06/2017 16:32, li li wrote:
>
>> Hi all,
>>In checking my R codes, I encountered the following problem. Is there a
>> way to fix this?
>> I tried to specify options(digits=). I did not fix the problem.
>>Thanks so much for your help!
>> Hanna
>>
>>
>> cdf(pmass)[2,2]==pcum[2,2][1] FALSE> cdf(pmass)[2,2][1] 0.758>
>>> pcum[2,2][1] 0.758
>>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posti
> ng-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] An R question

[li li]

Hi all,
  In checking my R codes, I encountered the following problem. Is there a
way to fix this?
I tried to specify options(digits=). I did not fix the problem.
  Thanks so much for your help!
   Hanna


> cdf(pmass)[2,2]==pcum[2,2][1] FALSE> cdf(pmass)[2,2][1] 0.758> 
> pcum[2,2][1] 0.758

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] subletting an array according to dimnames

[li li]

That works. Thank you!

2017-06-02 8:00 GMT-04:00 Adams, Jean :

> Have you tried P2["20", "10", "0"] ?
>
> Jean
>
> On Thu, Jun 1, 2017 at 3:10 PM, li li  wrote:
>
>> Hi all,
>>   I have a three dimensional array with the corresponding dimension names.
>> I would like to subset the array according to the dimension names. For
>> example,
>> suppose I want to extract the values corresponding to A=20, B=10, C=0.  I
>> know I
>> can do:
>>  P2[dimnames(P2)$A==20, dimnames(P2)$B==10, dimnames(P2)$C==0]
>>
>> But is there a better way for doing this? Thanks for your help!
>>   Hanna
>>
>> > dimnames(P2)
>>
>> $A
>>
>> [1] "20" "25" "30" "35" "40"
>>
>>
>> $B
>>
>> [1] "5"  "10" "15" "20" "25" "30" "35" "40"
>>
>>
>> $C
>>
>> [1] "0"  "5" "10" "15" "20" "25" "30" "35"
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] subletting an array according to dimnames

[li li]

Hi all,
  I have a three dimensional array with the corresponding dimension names.
I would like to subset the array according to the dimension names. For
example,
suppose I want to extract the values corresponding to A=20, B=10, C=0.  I
know I
can do:
 P2[dimnames(P2)$A==20, dimnames(P2)$B==10, dimnames(P2)$C==0]

But is there a better way for doing this? Thanks for your help!
  Hanna

> dimnames(P2)

$A

[1] "20" "25" "30" "35" "40"


$B

[1] "5"  "10" "15" "20" "25" "30" "35" "40"


$C

[1] "0"  "5" "10" "15" "20" "25" "30" "35"

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] [FORGED] Re: Question on function "scatterplot3d"

[li li]

Thanks so much for the help!!
   Hanna


2017-06-01 11:02 GMT-04:00 Uwe Ligges :

>
>
> On 01.06.2017 10:03, Rolf Turner wrote:
>
>> On 01/06/17 19:54, Uwe Ligges wrote:
>>
>>> A design flaw, whether the labels are cut depends somewhat on the sizce
>>> of the device, hence there is the argument
>>>
>>> y.margin.add
>>>
>>> add additional space between tick mark labels and axis label of the y
>>> axis
>>>
>>> for working around that limittation that can be set to some positive
>>> value
>>>
>>
>>
>> This seems to be addressing Hannah's (li li's) original enquiry, not my
>> follow-up in which I worried about the position, along the y-axis, of the
>> y-axis label.
>>
>
> Ah, that is intended as a smart way of rotating it along the axis is not
> easy (if not impossible) with the bas egraohics system.
>
> Best,
> Uwe Ligges
>
>
>
>
>> Or am I misunderstanding/missing something?
>>
>> cheers,
>>
>> Rolf
>>
>>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posti
> ng-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Question on function "scatterplot3d"

[li li]

Hi all,
  I have a question with regard to making plots using function
"scatterplot3d".
Please see the example below. It looks like, for y axis, the tickmark text
was cutoff.
The number "10" does not show up completely. I tried to work with par(mpg).
It does not
seem to work. Hope to get some advice here. Thanks much!
   Hanna




C <- runif(30)
B <- rep(1:3, each=10)
A <- rep(1:10,3)
scatterplot3d(B,A,C, type = "h", lwd = 1, pch = 16, color="red",  main =
"",
  grid=TRUE,  col.grid="lightgreen",
  xlab="x", ylab="y", zlab="z")

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] creat contingency tables with fixed row and column margins

[li li]

I meant that the post in your link above was NOT asked by me.

2017-05-27 16:53 GMT-04:00 li li :

> Hi Dave,
>   Thanks for your reply but the post in your link above was certainly
> asked by me.
>Hanna
>
> 2017-05-27 16:04 GMT-04:00 David Winsemius :
>
>>
>> > On May 27, 2017, at 12:49 PM, li li  wrote:
>> >
>> > Hi all,
>> >  Is there an R function that can be used to enumerate all the
>> contingency
>> > tables with fixed row and column margins. For example, can we list all
>> 3 by
>> > 3 tables with row margins 3,6,6 and column margins 5,5,5.
>>
>> The Posting Guide suggests that you describe the results of your efforts
>> at searching:
>>
>> https://www.google.com/search?q=use+r%3Alang+to+enumerate+co
>> ntingency+tables+with+fixed+margins&oq=use+r%3Alang+to+enume
>> rate+contingency+tables+with+fixed+margins&aqs=chrome..
>> 69i57.35014j0j4&sourceid=chrome&ie=UTF-8
>>
>> >   Thanks very much!
>> >   Hanna
>> >
>> >   [[alternative HTML version deleted]]
>>
>> And (yet again) you are asked to post in plain text.
>> >
>> > __
>> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide http://www.R-project.org/posti
>> ng-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>>
>> David Winsemius
>> Alameda, CA, USA
>>
>>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] creat contingency tables with fixed row and column margins

[li li]

Hi Dave,
  Thanks for your reply but the post in your link above was certainly asked
by me.
   Hanna

2017-05-27 16:04 GMT-04:00 David Winsemius :

>
> > On May 27, 2017, at 12:49 PM, li li  wrote:
> >
> > Hi all,
> >  Is there an R function that can be used to enumerate all the contingency
> > tables with fixed row and column margins. For example, can we list all 3
> by
> > 3 tables with row margins 3,6,6 and column margins 5,5,5.
>
> The Posting Guide suggests that you describe the results of your efforts
> at searching:
>
> https://www.google.com/search?q=use+r%3Alang+to+enumerate+
> contingency+tables+with+fixed+margins&oq=use+r%3Alang+to+
> enumerate+contingency+tables+with+fixed+margins&aqs=chrome.
> .69i57.35014j0j4&sourceid=chrome&ie=UTF-8
>
> >   Thanks very much!
> >   Hanna
> >
> >   [[alternative HTML version deleted]]
>
> And (yet again) you are asked to post in plain text.
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>
> David Winsemius
> Alameda, CA, USA
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] creat contingency tables with fixed row and column margins

[li li]

Hi all,
  Is there an R function that can be used to enumerate all the contingency
tables with fixed row and column margins. For example, can we list all 3 by
3 tables with row margins 3,6,6 and column margins 5,5,5.
   Thanks very much!
   Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Finding Infimum in R

[li li]

Yes. If the function f takes the value zero at some discontinuity point,
then the code gives the inf of the set I described.
Otherwise, it is an approximation since we need to worry about numerical
accuracy.

2017-04-10 14:08 GMT-04:00 Bert Gunter :

> Well, I haven't checked carefully, but  of course this does not find infs
> or sups at all, just mins or maxes in the sample, which are not the same.
>
> You'll have to do your own full testing  and debugging, however. I do not
> provide such a service.
>
> -- Bert
>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Finding Infimum in R

[li li]

Hi Burt,
Yes, the function is monotone increasing and points of discontinuity
are all known.
They are all numbers between 0 and 1.  Thanks very much!
   Hanna


2017-04-09 16:55 GMT-04:00 Bert Gunter :

> Details matter!
>
> 1. Are the points of discontinuity known? This is critical.
>
> 2. Can we assume monotonic increasing, as is shown?
>
>
> -- Bert
>
>
>
>
> Bert Gunter
>
> "The trouble with having an open mind is that people keep coming along
> and sticking things into it."
> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>
>
> On Sun, Apr 9, 2017 at 1:28 PM, li li  wrote:
> > Dear all,
> >   For a piecewise function F similar to the attached graph, I would like
> to
> > find
> > inf{x| F(x) >=0}.
> >
> >
> >  I tried to uniroot. It does not seem to work. Any suggestions?
> >  Thank you very much!!
> > Hanna
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Looking for

[li li]

Dear all,
  For a piecewise function F similar to the attached graph, I would like to
find
inf{x| F(x) >=0}.


 I tried to uniroot. It does not seem to work. Any suggestions?
 Thank you very much!!
Hanna


F.pdf
Description: Adobe PDF document
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] standard error for regression coefficients corresponding to factor levels

[li li]

Hi all,
  I have the following data called "data1". After fitting the ancova model
with different slopes and intercepts for each region, I calculated the
regression coefficients and the corresponding standard error. The standard
error (for intercept or for slope) are all the same for different regions.
Is there something wrong?
  I know the SE is related to (X^T X)^-1, where X is design matrix. So does
this happen whenever each factor level has the same set of values for
"week"?
 Thanks.
 Hanna



> mod <- lm(response ~ region*week, data1)> tmp <- coef(summary(mod))> res <- 
> matrix(NA, 5,4)> res[1,1:2] <- tmp[1,1:2]> res[2:5,1] <- tmp[1,1]+tmp[2:5,1]> 
> res[2:5,2] <- sqrt(tmp[2:5,2]^2-tmp[1,2]^2)> res[1,3:4] <- tmp[6,1:2]> 
> res[2:5,3] <- tmp[6,1]+tmp[7:10,1]> res[2:5,4] <- 
> sqrt(tmp[7:10,2]^2-tmp[6,2]^2)

> colnames(res) <- c("intercept", "intercept SE", "slope", "slope SE")> 
> rownames(res) <- letters[1:5]> res   intercept intercept SEslope   
> slope SE
a 0.18404464   0.08976301 -0.018629310 0.01385073
b 0.17605666   0.08976301 -0.022393789 0.01385073
c 0.16754130   0.08976301 -0.022367770 0.01385073
d 0.12554452   0.08976301 -0.017464385 0.01385073
e 0.06153256   0.08976301  0.007714685 0.01385073







> data1week region response
5  3  c  0.057325067
6  6  c  0.066723632
7  9  c -0.025317808
12 3  d  0.024692613
13 6  d  0.021761492
14 9  d -0.099820335
19 3  c  0.119559235
20 6  c -0.054456186
21 9  c  0.078811180
26 3  d  0.091667189
27 6  d -0.053400777
28 9  d  0.090754363
33 3  c  0.163818085
34 6  c  0.008959741
35 9  c -0.115410852
40 3  d  0.193920693
41 6  d -0.087738914
42 9  d  0.004987542
47 3  a  0.121332285
48 6  a -0.020202707
49 9  a  0.037295785
54 3  b  0.214304603
55 6  b -0.052346480
56 9  b  0.082501222
61 3  a  0.053540767
62 6  a -0.019182819
63 9  a -0.057629113
68 3  b  0.068592791
69 6  b -0.123298216
70 9  b -0.230671818
75 3  a  0.330741562
76 6  a  0.013902905
77 9  a  0.190620360
82 3  b  0.151002874
83 6  b  0.086177696
84 9  b  0.178982656
89 3  e  0.062974799
90 6  e  0.062035391
91 9  e  0.206200831
96 3  e  0.123102197
97 6  e  0.040181790
98 9  e  0.121332285
1033  e  0.147557564
1046  e  0.062035391
1059  e  0.144965770

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Test individual slope for each factor level in ANCOVA

[li li]

Hi John. Thanks much for your help. It is great to know this.
  Hanna

2017-03-16 8:02 GMT-04:00 Fox, John :

> Dear Hanna,
>
> You can test the slope in each non-reference group as a linear hypothesis.
> You didn’t make the data available for your example, so here’s an example
> using the linearHypothesis() function in the car package with the Moore
> data set in the same package:
>
> - - - snip - - -
>
> > library(car)
> > mod <- lm(conformity ~ fscore*partner.status, data=Moore)
> > summary(mod)
>
> Call:
> lm(formula = conformity ~ fscore * partner.status, data = Moore)
>
> Residuals:
> Min  1Q  Median  3Q Max
> -7.5296 -2.5984 -0.4473  2.0994 12.4704
>
> Coefficients:
>   Estimate Std. Error t value Pr(>|t|)
> (Intercept)   20.793483.26273   6.373 1.27e-07 ***
> fscore-0.151100.07171  -2.107  0.04127 *
> partner.statuslow-15.534084.40045  -3.530  0.00104 **
> fscore:partner.statuslow   0.261100.09700   2.692  0.01024 *
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Residual standard error: 4.562 on 41 degrees of freedom
> Multiple R-squared:  0.2942,Adjusted R-squared:  0.2426
> F-statistic: 5.698 on 3 and 41 DF,  p-value: 0.002347
>
> > linearHypothesis(mod, "fscore + fscore:partner.statuslow")
> Linear hypothesis test
>
> Hypothesis:
> fscore  + fscore:partner.statuslow = 0
>
> Model 1: restricted model
> Model 2: conformity ~ fscore * partner.status
>
>   Res.DfRSS Df Sum of Sq  F  Pr(>F)
> 1 42 912.45
> 2 41 853.42  159.037 2.8363 0.09976 .
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> - - - snip - - -
>
> In this case, there are just two levels for partner.status, but for a
> multi-level factor you can simply perform more than one test.
>
>
> I hope this helps,
>
>  John
>
> -----
> John Fox, Professor
> McMaster University
> Hamilton, Ontario, Canada
> Web: http://socserv.mcmaster.ca/jfox/
>
>
>
>
> On 2017-03-15, 9:43 PM, "R-help on behalf of li li"
>  wrote:
>
> >Hi all,
> >   Consider the data set where there are a continuous response variable, a
> >continuous predictor "weeks" and a categorical variable "region" with five
> >levels "a", "b", "c",
> >"d", "e".
> >  I fit the ANCOVA model as follows. Here the reference level is region
> >"a"
> >and there are 4 dummy variables. The interaction terms (in red below)
> >represent the slope
> >difference between each region and  the baseline region "a" and the
> >corresponding p-value is for testing whether this slope difference is
> >zero.
> >Is there a way to directly test whether the slope corresponding to each
> >individual factor level is 0 or not, instead of testing the slope
> >difference from the baseline level?
> >  Thanks very much.
> >  Hanna
> >
> >
> >
> >
> >
> >
> >> mod <- lm(response ~ weeks*region,data)> summary(mod)
> >Call:
> >lm(formula = response ~ weeks * region, data = data)
> >
> >Residuals:
> > Min   1Q   Median   3Q  Max
> >-0.19228 -0.07433 -0.01283  0.04439  0.24544
> >
> >Coefficients:
> >Estimate Std. Error t value Pr(>|t|)
> >(Intercept)1.2105556  0.0954567  12.682  1.2e-14 ***
> >weeks -0.021  0.0147293  -1.4480.156
> >regionb   -0.0257778  0.1349962  -0.1910.850
> >regionc   -0.034  0.1349962  -0.2550.800
> >regiond   -0.075  0.1349962  -0.5590.580
> >regione   -0.148  0.1349962  -1.0980.280weeks:regionb
> >-0.0007222  0.0208304  -0.0350.973
> >weeks:regionc -0.0017778  0.0208304  -0.0850.932
> >weeks:regiond  0.003  0.0208304   0.1440.886
> >weeks:regione  0.0301667  0.0208304   1.4480.156---
> >Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> >Residual standard error: 0.1082 on 35 degrees of freedom
> >Multiple R-squared:  0.2678,   Adjusted R-squared:  0.07946
> >F-statistic: 1.422 on 9 and 35 DF,  p-value: 0.2165
> >
> >   [[alternative HTML version deleted]]
> >
> >__
> >R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >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.
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Test individual slope for each factor level in ANCOVA

[li li]

Hi all,
   Consider the data set where there are a continuous response variable, a
continuous predictor "weeks" and a categorical variable "region" with five
levels "a", "b", "c",
"d", "e".
  I fit the ANCOVA model as follows. Here the reference level is region "a"
and there are 4 dummy variables. The interaction terms (in red below)
represent the slope
difference between each region and  the baseline region "a" and the
corresponding p-value is for testing whether this slope difference is zero.
Is there a way to directly test whether the slope corresponding to each
individual factor level is 0 or not, instead of testing the slope
difference from the baseline level?
  Thanks very much.
  Hanna






> mod <- lm(response ~ weeks*region,data)> summary(mod)
Call:
lm(formula = response ~ weeks * region, data = data)

Residuals:
 Min   1Q   Median   3Q  Max
-0.19228 -0.07433 -0.01283  0.04439  0.24544

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)1.2105556  0.0954567  12.682  1.2e-14 ***
weeks -0.021  0.0147293  -1.4480.156
regionb   -0.0257778  0.1349962  -0.1910.850
regionc   -0.034  0.1349962  -0.2550.800
regiond   -0.075  0.1349962  -0.5590.580
regione   -0.148  0.1349962  -1.0980.280weeks:regionb
-0.0007222  0.0208304  -0.0350.973
weeks:regionc -0.0017778  0.0208304  -0.0850.932
weeks:regiond  0.003  0.0208304   0.1440.886
weeks:regione  0.0301667  0.0208304   1.4480.156---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1082 on 35 degrees of freedom
Multiple R-squared:  0.2678,Adjusted R-squared:  0.07946
F-statistic: 1.422 on 9 and 35 DF,  p-value: 0.2165

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] nested structure for Ancova

[li li]

Thanks so much, Richard. That  works.
 Hanna

2017-03-12 14:31 GMT-04:00 Richard M. Heiberger :

> I think you need either
>
> mod4 <- lm( y ~ -1 + area / (month*group), data=two)
>
> mod5 <- lm( y ~ area / (month*group), data=two)
>
> With either of those,  area:month and area:group and Residuals add up.
>
>
> On Sun, Mar 12, 2017 at 10:39 li li  wrote:
>
>> Hi All,
>>   I have a dataset which contains 4 variables: area, group, time, y,
>> Area is a factor that has two levels A and B, group is a factor that is
>> nested within area. There are four groups within each area.
>> y is the response variable, and time refers to different days.
>>
>>   Below is the how data looks like.
>>
>>   First, I fit separate ancova for each area (area A and B), For each
>> area,
>> I obtained different regression lines for each group within that area.
>>
>>   mod1 <- lm(y ~ month * group, data=two[two$area=="A"]
>>
>>   mod2 <- lm(y ~ month * group, data=two[two$area=="B"]
>>
>>  I want to fit the model at one time by using the nested structure
>>
>>  mod3 <-  lm(y ~ -1+ month * (area/group), data=two)
>>
>> I get different results using mod 3 from using mod1 and mod2.
>>
>>   Can someone give some suggestion on this.  How can I specify the model
>> in
>> R to simultaneously fit the model to get the same results as from mod1 and
>> mod2.
>>
>> Thanks  very much!
>> Hanna
>>
>> > head(two[two$area=="A",])   group time y area
>> 79 13 -1.394327A
>> 80 16 -1.435485A
>> 81 19 -1.406497A
>> 82 1   12 -1.265848A
>> 83 10 -1.316768A
>> 84 16 -1.431292A> head(two[two$area=="B",])  group time
>>  y area
>> 1 10 -2.145581B
>> 2 10 -1.910543B
>> 3 10 -2.128632B
>> 4 13 -2.079442B
>> 5 13 -2.273026B
>> 6 16 -2.312635B
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] nested structure for Ancova

[li li]

Hi All,
  I have a dataset which contains 4 variables: area, group, time, y,
Area is a factor that has two levels A and B, group is a factor that is
nested within area. There are four groups within each area.
y is the response variable, and time refers to different days.

  Below is the how data looks like.

  First, I fit separate ancova for each area (area A and B), For each area,
I obtained different regression lines for each group within that area.

  mod1 <- lm(y ~ month * group, data=two[two$area=="A"]

  mod2 <- lm(y ~ month * group, data=two[two$area=="B"]

 I want to fit the model at one time by using the nested structure

 mod3 <-  lm(y ~ -1+ month * (area/group), data=two)

I get different results using mod 3 from using mod1 and mod2.

  Can someone give some suggestion on this.  How can I specify the model in
R to simultaneously fit the model to get the same results as from mod1 and
mod2.

Thanks  very much!
Hanna

> head(two[two$area=="A",])   group time y area
79 13 -1.394327A
80 16 -1.435485A
81 19 -1.406497A
82 1   12 -1.265848A
83 10 -1.316768A
84 16 -1.431292A> head(two[two$area=="B",])  group time
 y area
1 10 -2.145581B
2 10 -1.910543B
3 10 -2.128632B
4 13 -2.079442B
5 13 -2.273026B
6 16 -2.312635B

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] weighted sum of independent chi square random variables

[li li]

Hi all,
   Is there a function in R that can calculate the quantiles or percentiles
for the weighted sum of independent chi square random variables. I found a
few functions for calculating probability distribution function for such
random variables (e.g. pchisqsum..), but can not find any function for
finding quantiles.
   Thanks in advance.
Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] ancova for dependent response

[li li]

Hi all,
  Is there an R function which can handles dependent response in Analysis
of covariance model. The dependence structure is known and is there to
account for it in ANCOVA analysis in R?
  Thanks.
  Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] drm function in drc package for fitting dose response curve

[li li]

Hi all,
  When using drm function in drc package, we can fit several does response
curves simultaneously or fitting each curve separately. If we fit
simultaneous curves without any constraining condition, for example,
parallelism condition, will we get the same results as the results obtained
by fitting separate curves, particularly parameter estimates?
If the results are different, what is the root differences, underlying
model and so on?
Thanks much.
   Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] plot.drm in "drc" package

[li li]

Thanks Petr for the reply.
When I run "plot(mod, type="all",log="x")", the tickmarks of the x axix
include (0.02, 0.1, 1,10).
But the log(Dose) should be less than 4.
Anyway, I think there is something missing in the plot.drm function.

   Hanna

2016-08-31 7:41 GMT-04:00 PIKAL Petr :

> Hi
>
> Thanks for code.
> see in line
>
> > -Original Message-
> > From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of li li
> > Sent: Tuesday, August 30, 2016 5:07 PM
> > To: r-help 
> > Subject: [R] plot.drm in "drc" package
> >
> > Hi all,
> >   I am trying to use the drc package to fit 4L curve.
> > I am so confused about the plot.drm function in the drc package.
> > Particularly, I am confused about the scale of the xaxis in the plots
> generated
> > using the plot.drm function. See the example below:
> >
> > ## generate data and fit the model
> > dose <- rep(50*2^(-(0:11)),3)
> > dose
> > d <- 100
> > c <- 1
> > b <- 1
> > e <- 1.6
> > y <- rnorm(length(dose))+ c+ (d-c)/(1+exp(b*(log(dose)-log(e
> > library(drc)
> > mod <- drm(y~dose, fct = LL.4())
> > summary(mod)
> >
> > Now I plot the data and the fitted curve with the plot.drm using the code
> > below and get the figure 1 below.
> >
> >
> > ##obtaining figure 1
> > plot(mod, type="all",log="x")
> >
> >
> > Next I plot the raw data and add the curve by extracting the estimate of
> the
> > parameters.
> >
> > ##extract parameters
> > para <- mod$fit$par
> > bhat<- para[1]
> > chat <- para[2]
> > dhat <- para[3]
> > ehat <- para[4]
> >
> > ##plot figure 2
> > plot(log(dose),y)
> > points(log(50*2^(-(0:11))),  chat +
> > (dhat-chat)/(1+exp(bhat*(log(50*2^(-(0:11)))-log(ehat, type="l")
> >
> > My question is regarding the figure 1 generated by the plot.drm.
> > The x axis is the not the log scale of the doses. I checked the package
> manual,
>  ^^^
> Well I am either missing something obvios but when I tried
>
> plot(mod, type="all",log="x")
>
> I got x axis with log scaling. Actually the result is the same as
>
> plot(mod, type="all")
>
> If you want x axis to be in original range you need to put empty string
>
> plot(mod, type="all",log="")
>
> Cheers
> Petr
>
> > it says the default is log base 10. But it is not true in this case.
> > Does some have some insight on the correct usage of the plot.drm
> function.
> >
> > Thanks much in advance.
> >   Hanna
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>
> 
> Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou
> určeny pouze jeho adresátům.
> Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě
> neprodleně jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie
> vymažte ze svého systému.
> Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email
> jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
> Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi
> či zpožděním přenosu e-mailu.
>
> V případě, že je tento e-mail součástí obchodního jednání:
> - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření
> smlouvy, a to z jakéhokoliv důvodu i bez uvedení důvodu.
> - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout;
> Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany
> příjemce s dodatkem či odchylkou.
> - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve
> výslovným dosažením shody na všech jejích náležitostech.
> - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za
> společnost žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn
> nebo písemně pověřen a takové pověření nebo plná moc byly adresátovi tohoto
> emailu případně osobě, kterou adresát zastupuje, předloženy nebo jejich
> existence je adresátovi či osobě jím zastoupené známá.
>
> This e-mail and any documents attached to it may be confidential and are
> inte

Re: [R] plot with different symbols and colors according to the factor levels

[li li]

Thank you all.
   Hanna

2016-08-30 16:55 GMT-04:00 Paulo Moniz :

> Obter o Outlook para Android <https://aka.ms/ghei36>
>
>
>
> On Tue, Aug 30, 2016 at 2:41 PM -0300, "Clint Bowman" 
> wrote:
>
> Hanna,
>
> lili<-read.table("lili.txt",header=T)  # don't forget to label the row
> number if it's in your data
>
> with(lili,plot(y,conc,pch=sample,col=sample))
>
> Clint
>
>
> Clint BowmanINTERNET:   cl...@ecy.wa.gov
> Air Quality Modeler INTERNET:   cl...@math.utah.edu
> Department of Ecology   VOICE:  (360) 407-6815
> PO Box 47600FAX:(360) 407-7534
> Olympia, WA 98504-7600
>
>  USPS:   PO Box 47600, Olympia, WA 98504-7600
>  Parcels:300 Desmond Drive, Lacey, WA 98503-1274
>
> On Tue, 30 Aug 2016, li li wrote:
>
> > Hi all,
> > I have the following data. I want to plot the data (y ~ conc)
> > with different symbols and colors corresponding to different levels of
> the
> > factor sample.
> > I could create a column with color and pch and then do the plot, but I am
> > sure there are much better ways.
> > Can anyone make suggestions?
> >  Hanna
> >
> >
> >
> >   y conc sample
> > 1  33 20.0  1
> > 2  33  5.0  1
> > 3  35  1.25000  1
> > 4  43  0.31250  1
> > 5  58  0.078125000  1
> > 6  54  0.019531250  1
> > 7  57  0.004882812  1
> > 8  57  0.001220703  1
> > 9  32 20.0  1
> > 10 32  5.0  1
> > 11 34  1.25000  1
> > 12 52  0.31250  1
> > 13 57  0.078125000  1
> > 14 58  0.019531250  1
> > 15 59  0.004882812  1
> > 16 50  0.001220703  1
> > 17 34 20.0  2
> > 18 34  5.0  2
> > 19 38  1.25000  2
> > 20 53  0.31250  2
> > 21 57  0.078125000  2
> > 22 57  0.019531250  2
> > 23 57  0.004882812  2
> > 24 52  0.001220703  2
> > 25 34 20.0  2
> > 26 33  5.0  2
> > 27 36  1.25000  2
> > 28 48  0.31250  2
> > 29 58  0.078125000  2
> > 30 57  0.019531250  2
> > 31 58  0.004882812  2
> > 32 53  0.001220703  2
> > 33 34 20.0  2
> > 34 35  5.0  2
> > 35 37  1.25000  2
> > 36 49  0.31250  2
> > 37 55  0.078125000  2
> > 38 59  0.019531250  2
> > 39 57  0.004882812  2
> > 40 54  0.001220703  2
> > 41 36 20.0  3
> > 42 33  5.0  3
> > 43 36  1.25000  3
> > 44 51  0.31250  3
> > 45 57  0.078125000  3
> > 46 57  0.019531250  3
> > 47 59  0.004882812  3
> > 48 56  0.001220703  3
> > 49 33 20.0  3
> > 50 32  5.0  3
> > 51 35  1.25000  3
> > 52 47  0.31250  3
> > 53 57  0.078125000  3
> > 54 56  0.019531250  3
> > 55 57  0.004882812  3
> > 56 53  0.001220703  3
> > 57 33 20.0  3
> > 58 34  5.0  3
> > 59 38  1.25000  3
> > 60 52  0.31250  3
> > 61 56  0.078125000  3
> > 62 61  0.019531250  3
> > 63 56  0.004882812  3
> > 64 55  0.001220703  3
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
> >
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] plot with different symbols and colors according to the factor levels

[li li]

Hi all,
 I have the following data. I want to plot the data (y ~ conc)
with different symbols and colors corresponding to different levels of the
factor sample.
I could create a column with color and pch and then do the plot, but I am
sure there are much better ways.
Can anyone make suggestions?
  Hanna



   y conc sample
1  33 20.0  1
2  33  5.0  1
3  35  1.25000  1
4  43  0.31250  1
5  58  0.078125000  1
6  54  0.019531250  1
7  57  0.004882812  1
8  57  0.001220703  1
9  32 20.0  1
10 32  5.0  1
11 34  1.25000  1
12 52  0.31250  1
13 57  0.078125000  1
14 58  0.019531250  1
15 59  0.004882812  1
16 50  0.001220703  1
17 34 20.0  2
18 34  5.0  2
19 38  1.25000  2
20 53  0.31250  2
21 57  0.078125000  2
22 57  0.019531250  2
23 57  0.004882812  2
24 52  0.001220703  2
25 34 20.0  2
26 33  5.0  2
27 36  1.25000  2
28 48  0.31250  2
29 58  0.078125000  2
30 57  0.019531250  2
31 58  0.004882812  2
32 53  0.001220703  2
33 34 20.0  2
34 35  5.0  2
35 37  1.25000  2
36 49  0.31250  2
37 55  0.078125000  2
38 59  0.019531250  2
39 57  0.004882812  2
40 54  0.001220703  2
41 36 20.0  3
42 33  5.0  3
43 36  1.25000  3
44 51  0.31250  3
45 57  0.078125000  3
46 57  0.019531250  3
47 59  0.004882812  3
48 56  0.001220703  3
49 33 20.0  3
50 32  5.0  3
51 35  1.25000  3
52 47  0.31250  3
53 57  0.078125000  3
54 56  0.019531250  3
55 57  0.004882812  3
56 53  0.001220703  3
57 33 20.0  3
58 34  5.0  3
59 38  1.25000  3
60 52  0.31250  3
61 56  0.078125000  3
62 61  0.019531250  3
63 56  0.004882812  3
64 55  0.001220703  3

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] plot.drm in "drc" package

[li li]

Hi all,
  I am trying to use the drc package to fit 4L curve.
I am so confused about the plot.drm function in the drc package.
Particularly, I am confused about the scale of the xaxis in the plots
generated using the plot.drm function. See the example below:

## generate data and fit the model
dose <- rep(50*2^(-(0:11)),3)
dose
d <- 100
c <- 1
b <- 1
e <- 1.6
y <- rnorm(length(dose))+ c+ (d-c)/(1+exp(b*(log(dose)-log(e
library(drc)
mod <- drm(y~dose, fct = LL.4())
summary(mod)

Now I plot the data and the fitted curve with the plot.drm using the code
below and get the figure 1 below.


##obtaining figure 1
plot(mod, type="all",log="x")


Next I plot the raw data and add the curve by extracting the estimate of
the parameters.

##extract parameters
para <- mod$fit$par
bhat<- para[1]
chat <- para[2]
dhat <- para[3]
ehat <- para[4]

##plot figure 2
plot(log(dose),y)
points(log(50*2^(-(0:11))),  chat +
(dhat-chat)/(1+exp(bhat*(log(50*2^(-(0:11)))-log(ehat, type="l")

My question is regarding the figure 1 generated by the plot.drm.
The x axis is the not the log scale of the doses. I checked the package
manual, it says the default is log base 10. But it is not true in this case.
Does some have some insight on the correct usage of the plot.drm function.

Thanks much in advance.
  Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Help on improving an algorithm

[li li]

Thanks very much Don! You alternative method does make it much faster.
However I think I made a mistake in my previous coding. Sorry.

Instead of checking whether each row sum is larger than 8, what I wanted is
to
check whether each column sum is lager than 8. See the code below.

So in this case, I should still try to use expand.grid? How to use
expand.grid if each element is a row instead of a single number or
character?

Thanks for you help.


> dim(tmp_a);dim(tmp_b); dim(tmp_c); dim(tmp_d)[1] 252   6
[1] 126   6
[1] 126   6
[1] 126   6


p <- 0
for (i in 1:dim(tmp_a)[1]){
for (j in 1:dim(tmp_b)[1]){
for (k in 1:dim(tmp_c)[1]){
for (l in 1:dim(tmp_c)[1]){
conti <- rbind(tmp_a[i,], tmp_b[j,], tmp_c[k,],tmp_d[l,])

* if (sum(apply(conti,2,sum)>8)==0)*p <- p+1
print(p)


2016-07-28 12:16 GMT-04:00 MacQueen, Don :

> This is a good example to illustrate that R is a vectorized language, and
> programming concepts that would work in fortran, C, or any of many other
> languages are not always effective in R.
>
> The vectorized method below has produced the same value for 'p' in every
> example I¹ve tried, and is much faster.
>
> Note that the desired calculation looks at every combination of one row
> from each of four matrices, i.e., 126^4 combinations (a large number).
> This suggests the use of expand.grid().
>
> I don't see any need to actually construct the 'conti' matrix. Only its
> row sums are needed, and those can be calculated directly from the rows of
> the input matrices.
>
> The vectorized method makes using 126 rows tractable. It still takes a
> while on my machine, but it's tolerable. The original method takes longer
> than I'm willing to wait. I used a smaller number of rows for testing.
>
> ## number of rows in each input matrix
> nr <- 126
> nr <- 11
> ## number of columns in each input matrix
> nc <- 6
>
> ## simulate input data, in order to create
> ## a reproducible example
> xr <- 0:2
> tmpa <- matrix( sample(xr, nr*nc, replace=TRUE) , ncol=nc)
> tmpb <- matrix( sample(xr, nr*nc, replace=TRUE) , ncol=nc)
> tmpc <- matrix( sample(xr, nr*nc, replace=TRUE) , ncol=nc)
> tmpd <- matrix( sample(xr, nr*nc, replace=TRUE) , ncol=nc)
>
> cat('method 2 (vectorized)\n')
> t2 <- system.time({
>   rsa <- rowSums(tmpa)
>   rsb <- rowSums(tmpb)
>   rsc <- rowSums(tmpc)
>   rsd <- rowSums(tmpd)
>
>   bigr <- expand.grid(rsa, rsb, rsc, rsd)
>
>   p2 <- sum( rowSums( bigr > 8) == 0)
>   print(p2)
> })
>
> cat('original method (looping)\n')
>
> torig <- system.time({
>   p <- 0
>   for (i in 1:nrow(tmpa)) {
> for (j in 1:nrow(tmpb)) {
>   for (k in 1:nrow(tmpc)) {
> for (l in 1:nrow(tmpd)) {
>   conti <- rbind(tmpa[i,], tmpb[j,], tmpc[k,],tmpd[l,])
>   if (sum(apply(conti,1,sum)>8)==0) p <- p+1
>   ##print(p)
> }
>   }
> }
>   }
>
> print(p)
> })
>
> cat('\n')
> print(t2)
> print(torig)
>
>
> A couple of minor points:
>   dim(tmp_a)[1] can be replaced by nrow(tmp_a)
>   in the original code, apply() can be replaced by rowSums()
> It might be faster if rowSums() was replaced by .rowSums(); I haven't
> tried that.
>
> -Don
>
>
> --
> Don MacQueen
>
> Lawrence Livermore National Laboratory
> 7000 East Ave., L-627
> Livermore, CA 94550
> 925-423-1062
>
>
>
>
>
> On 7/27/16, 1:41 PM, "R-help on behalf of li li"
>  wrote:
>
> >Hi all,
> > I have four matrix tmp_a, tmp_b, tmp_c and tmp_d whose dimensions are
> >shown as below.
> >I want to take one row from each of these matrices and then put the four
> >selected rows into a matrix.
> >I want to count the number of such matrices for which the vector of row
> >sum
> >is less than or equal to (8,8,8,8,8) (componentwise).
> >Below is the code I use now. Is there a way to make this more efficient
> >and
> >faster?
> > Thanks for the help in advance.
> >
> >> dim(tmp_a);dim(tmp_b); dim(tmp_c); dim(tmp_d)[1] 252   6
> >[1] 126   6
> >[1] 126   6
> >[1] 126   6
> >
> >
> >p <- 0
> >for (i in 1:dim(tmp_a)[1]){
> >for (j in 1:dim(tmp_b)[1]){
> >for (k in 1:dim(tmp_c)[1]){
> >for (l in 1:dim(tmp_c)[1]){
> >conti <- rbind(tmp_a[i,], tmp_b[j,], tmp_c[k,],tmp_d[l,])
> >if (sum(apply(conti,1,sum)>8)==0)
> >p <- p+1
> >print(p)
> >
> >

[R] Help on improving an algorithm

[li li]

Hi all,
 I have four matrix tmp_a, tmp_b, tmp_c and tmp_d whose dimensions are
shown as below.
I want to take one row from each of these matrices and then put the four
selected rows into a matrix.
I want to count the number of such matrices for which the vector of row sum
is less than or equal to (8,8,8,8,8) (componentwise).
Below is the code I use now. Is there a way to make this more efficient and
faster?
 Thanks for the help in advance.

> dim(tmp_a);dim(tmp_b); dim(tmp_c); dim(tmp_d)[1] 252   6
[1] 126   6
[1] 126   6
[1] 126   6


p <- 0
for (i in 1:dim(tmp_a)[1]){
for (j in 1:dim(tmp_b)[1]){
for (k in 1:dim(tmp_c)[1]){
for (l in 1:dim(tmp_c)[1]){
conti <- rbind(tmp_a[i,], tmp_b[j,], tmp_c[k,],tmp_d[l,])
if (sum(apply(conti,1,sum)>8)==0)
p <- p+1
print(p)

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Fixed Effects in lme function

[li li]

 Dear all,
For the data below, I would like to fit a model with common random
slope and  common random intercept as shown below. I am interested in
obtaining
separate fixed effect estimates (intercept and slope and corresponding
hypothesis test)
for each method. Instead of performing the analysis method by method, can I
use the syntax as below, specifically, the part "fixed = response ~
method/time"?
I know this is legitimate specification if there is no random effects
involved.
   Thanks so much in advance!
  Hanna

lme(fixed= response ~ method/time, random=~ 1+time | lot, data=dat,
weights= varIdent(form=~1|method),   control = lmeControl(opt = "optim"),
na.action = na.exclude)



  response individual time method
1102.9  30  3
2103.0  33  3
3103.0  36  3
4102.8  39  3
5102.2  3   12  3
6102.5  3   15  3
7103.0  3   18  3
8102.0  3   24  3
9102.8  10  3
10   102.7  13  3
11   103.0  16  3
12   102.2  19  3
13   103.0  1   12  3
14   102.8  1   15  3
15   102.8  1   18  3
16   102.9  1   24  3
17   102.2  20  3
18   102.6  23  3
19   103.4  26  3
20   102.3  29  3
21   101.3  2   12  3
22   102.1  2   15  3
23   102.1  2   18  3
24   102.2  2   24  3
25   102.7  40  3
26   102.3  43  3
27   102.6  46  3
28   102.7  49  3
29   102.8  4   12  3
30   102.5  50  3
31   102.4  53  3
32   102.1  56  3
33   102.3  60  3
34   102.3  63  3
35   101.9  70  3
36   102.0  73  3
37   107.4  30  1
38   101.3  3   12  1
3992.8  3   15  1
4073.7  3   18  1
41   104.7  3   24  1
4292.6  10  1
43   101.9  1   12  1
44   106.3  1   15  1
45   104.1  1   18  1
4695.6  1   24  1
4779.8  20  1
4889.7  2   12  1
4997.0  2   15  1
50   108.4  2   18  1
51   103.5  2   24  1
5296.4  40  1
5389.3  4   12  1
54   112.6  50  1
5593.3  60  1
5699.6  70  1
57   109.5  30  2
5898.5  3   12  2
59   103.5  3   24  2
60   113.5  10  2
6194.5  1   12  2
6288.5  1   24  2
6399.5  20  2
6497.5  2   12  2
6598.5  2   24  2
66   103.5  40  2
6789.5  50  2
6887.5  60  2
6982.5  70  2

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] [R-sig-ME] multcomp package

[li li]

0.279301  0.7810
method2:time -0.12616  0.360585 57 -0.349877  0.7277
method3:time -0.08010  0.251105 57 -0.318999  0.7509
 Correlation:
 (Intr) methd2 methd3 time   mthd2:
method2  -0.726
method3  -0.999  0.726
time -0.779  0.566  0.779
method2:time  0.542 -0.712 -0.542 -0.696
method3:time  0.778 -0.566 -0.779 -0.999  0.696

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max
-2.67575293 -0.51633192  0.06742723  0.59706762  2.81061874

Number of Observations: 69
Number of Groups: 7










2016-06-06 11:21 GMT-04:00 Gabriel Baud-Bovy :

> On 06/06/2016 4:57 PM, li li wrote:
>
> Hi all,
>   After fitting a random slope and random intercept model using lme
> function, I want
> to test whether each of the fixed slopes is equal to zero (The output of
> model is below).
> Can this be done (testing each individual slope) using multcomp package?
>
> I don't understand what you mean by testing individual slope ?
>
> For the fixed effects, you might test whether there is a method, time or
> interaction
> effect using one of the methods described below
>
> For the randome effects, according to your model specification, the time
> dependency might vary for each individual. the sd for the time
> (0.0001841179)
> is small.  You might want to test whether to include randome slope  by
> doing a LRT
> between a model with it and another model without.
>
> Whay not include method in the random effects ?
>
> Gabriel
>
>   Thanks much for the help.
>Hanna
>
> To get p values:
>
>
> http://stats.stackexchange.com/questions/118416/getting-p-value-with-mixed-effect-with-lme4-package
>
>
> http://mindingthebrain.blogspot.it/2014/02/three-ways-to-get-parameter-specific-p.html
>
> Using lmerTest package
> https://cran.r-project.org/web/packages/lmerTest/lmerTest.pdf
>
> or using mixed in afex package
> http://rpackages.ianhowson.com/cran/afex/man/mixed.html
>
> both use pbkrtest packages
> https://cran.r-project.org/web/packages/pbkrtest/pbkrtest.pdf
>
> a faq
> http://glmm.wikidot.com/faq
>
>
>
>
> summary(mod1)Linear mixed-effects model fit by REML
>
>  Data: one
>AIC  BIClogLik
>   304.4703 330.1879 -140.2352
>
> Random effects:
>  Formula: ~1 + time | individual
>  Structure: General positive-definite, Log-Cholesky parametrization
> StdDev   Corr
> (Intercept) 0.2487869075 (Intr)
> time0.0001841179 -0.056
> Residual0.3718305953
>
> Variance function:
>  Structure: Different standard deviations per stratum
>  Formula: ~1 | method
>  Parameter estimates:
>312
>  1.0 26.59750 24.74476
> Fixed effects: reponse ~ method * time
> Value Std.Error DF   t-value p-value(Intercept)
> 96.65395  3.528586 57 27.391694  0.
> method2   1.17851  4.856026 57  0.242689  0.8091
> method3   5.87505  3.528617 57  1.664973  0.1014time
> 0.07010  0.250983 57  0.279301  0.7810
> method2:time -0.12616  0.360585 57 -0.349877  0.7277
> method3:time -0.08010  0.251105 57 -0.318999  0.7509
>  Correlation:
>  (Intr) methd2 methd3 time   mthd2:
> method2  -0.726
> method3  -0.999  0.726
> time -0.779  0.566  0.779
> method2:time  0.542 -0.712 -0.542 -0.696
> method3:time  0.778 -0.566 -0.779 -0.999  0.696
>
> Standardized Within-Group Residuals:
> Min  Q1 Med  Q3 Max
> -2.67575293 -0.51633192  0.06742723  0.59706762  2.81061874
>
> Number of Observations: 69
> Number of Groups: 7 >
>
>   [[alternative HTML version deleted]]
>
> ___r-sig-mixed-mod...@r-project.org
>  mailing listhttps://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> .
>
>
>
>
> --
> -
> 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
> -
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] multcomp package

[li li]

Hi all,
  After fitting a random slope and random intercept model using lme
function, I want
to test whether each of the fixed slopes is equal to zero (The output of
model is below).
Can this be done (testing each individual slope) using multcomp package?
  Thanks much for the help.
   Hanna

> summary(mod1)Linear mixed-effects model fit by REML
 Data: one
   AIC  BIClogLik
  304.4703 330.1879 -140.2352

Random effects:
 Formula: ~1 + time | individual
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 0.2487869075 (Intr)
time0.0001841179 -0.056
Residual0.3718305953

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | method
 Parameter estimates:
   312
 1.0 26.59750 24.74476
Fixed effects: reponse ~ method * time
Value Std.Error DF   t-value p-value(Intercept)
96.65395  3.528586 57 27.391694  0.
method2   1.17851  4.856026 57  0.242689  0.8091
method3   5.87505  3.528617 57  1.664973  0.1014time
0.07010  0.250983 57  0.279301  0.7810
method2:time -0.12616  0.360585 57 -0.349877  0.7277
method3:time -0.08010  0.251105 57 -0.318999  0.7509
 Correlation:
 (Intr) methd2 methd3 time   mthd2:
method2  -0.726
method3  -0.999  0.726
time -0.779  0.566  0.779
method2:time  0.542 -0.712 -0.542 -0.696
method3:time  0.778 -0.566 -0.779 -0.999  0.696

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max
-2.67575293 -0.51633192  0.06742723  0.59706762  2.81061874

Number of Observations: 69
Number of Groups: 7 >

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Fwd: model specification using lme

[li li]

Thanks Thierry for the reply. I think I now have a better understanding for
the specification of the random effects when using lme function.
Are my interpretations below correct?

random=~ 1 | individual   (same random intercept no random slope)

random=~ 1 +method| individual(same random intercept and same random
slope)

random=~ 1 +method:time| individual(same random intercept and different
random slope for different method)
random=~ 1 +method + method:time| individual(different random intercept
and different random slope for different method

 The summary results from the lme function shows whether the slopes for the
three methods are equal (parallelism). I also wanted to test the hypotheses
that each of the fixed slopes (corresponding to the three methods) equals
0, can I use multicomp package for that purpose? I am confused on how to
make correct specifications in glht function to test these hypotheses.

Hanna


> summary(mod1)
Linear mixed-effects model fit by REML
 Data: one
   AIC  BIClogLik
  304.4703 330.1879 -140.2352

Random effects:
 Formula: ~1 + time | individual
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 0.2487869075 (Intr)
time0.0001841179 -0.056
Residual0.3718305953

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | method
 Parameter estimates:
   312
 1.0 26.59750 24.74476
Fixed effects: reponse ~ method * time
Value Std.Error DF   t-value p-value(Intercept)
96.65395  3.528586 57 27.391694  0.
method2   1.17851  4.856026 57  0.242689  0.8091
method3   5.87505  3.528617 57  1.664973  0.1014
*time  0.07010  0.250983 57  0.279301  0.7810 method2:time
-0.12616  0.360585 57 -0.349877  0.7277 method3:time -0.08010  0.251105 57
-0.318999  0.7509*
 Correlation:
 (Intr) methd2 methd3 time   mthd2:
method2  -0.726
method3  -0.999  0.726
time -0.779  0.566  0.779
method2:time  0.542 -0.712 -0.542 -0.696
method3:time  0.778 -0.566 -0.779 -0.999  0.696

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max
-2.67575293 -0.51633192  0.06742723  0.59706762  2.81061874

Number of Observations: 69
Number of Groups: 7 >



-- Forwarded message --
From: Thierry Onkelinx 
Date: 2016-05-30 4:40 GMT-04:00
Subject: Re: [R] model specification using lme
To: li li 
Cc: r-help 



Dear Hanna,

None of the models are correct is you want the same random intercept for
the different methods but different random slope per method.

You can random = ~ 1 + time:method | individual

The easiest way to get alpha_0 and tau_i is to apply post-hoc contrasts.
That is fairly easy to do with the multcomp package.

alpha_0 = (m1 + m2 + m3) / 3
m1 = intercept
m2 = intercept + method2
m3 = intercept + method3
hence alpha_0 = intercept + method2/3 + method3/3

m1 = alpha_0 + tau_1
tau_1 = intercept - method2/3 - method3/3

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2016-05-29 21:23 GMT+02:00 li li :


> Hi all,
>   For the following data, I consider the following random intercept and
> random slope model. Denote as y_ijk the response value from *j*th
> individual within *i*th method at time point *k*. Assume the following
> model for y_ijk:
>
>   y_ijk= (alpha_0+ tau_i +a_j(i))+(beta_i+b_j(i)) T_k + e_ijk
>
>
> Here alpha_0 is the grand mean;
>   tau_i is the fixed effect for ith method;
>   a_j(i) is random intercept corresponding to the *j*th individual
> within *i*th method, assumed to be common for all three methods;
>   beta_i is the fixed slope corresponding to the ith method;
>   b_j(i) is the random slope corresponding to jth individual for
> the ith method, assumed to be different for different methods;
>   T_k is the time corresponding to y_ijk;
>   e_ijk is the residual.
>
> For this model, I consider the three specification using  the lme function
> as follows:
>
>
> mod1 <- lme(fixed= reponse ~ method*time, random=~ 1 +time | individual,
> data=one, weights= varIdent(form=~1|method),
> control = lmeControl(opt = "optim"))
>
> mod2 <- lme(fixed= reponse 

[R] model specification using lme

[li li]

Hi all,
  For the following data, I consider the following random intercept and
random slope model. Denote as y_ijk the response value from *j*th
individual within *i*th method at time point *k*. Assume the following
model for y_ijk:

  y_ijk= (alpha_0+ tau_i +a_j(i))+(beta_i+b_j(i)) T_k + e_ijk


Here alpha_0 is the grand mean;
  tau_i is the fixed effect for ith method;
  a_j(i) is random intercept corresponding to the *j*th individual
within *i*th method, assumed to be common for all three methods;
  beta_i is the fixed slope corresponding to the ith method;
  b_j(i) is the random slope corresponding to jth individual for
the ith method, assumed to be different for different methods;
  T_k is the time corresponding to y_ijk;
  e_ijk is the residual.

For this model, I consider the three specification using  the lme function
as follows:


mod1 <- lme(fixed= reponse ~ method*time, random=~ 1 +time | individual,
data=one, weights= varIdent(form=~1|method),
control = lmeControl(opt = "optim"))

mod2 <- lme(fixed= reponse ~ method*time, random=~ 0 +time | individual,
data=one, weights= varIdent(form=~1|method),
control = lmeControl(opt = "optim"))

mod3 <- lme(fixed= reponse ~ method*time, random=~ method +time |
individual, data=one, weights= varIdent(form=~1|method),
control = lmeControl(opt = "optim"))

I think mod1 is the correct one. However, I am kind of confused with the
right usage of lme function. Can someone familiar with this give some help
here?

Another question is regarding the fixed effect   tau_1, tau_2 and tau_3
(corresponding to the three methods). One main question I am interested in
is whether each of them are statistically different from zero. In the
summary results below (shaded part), it looks that the result for method 2
and 3 are given with reference to method 1). Is there a way to obtain
specific result separately for alpha_0 (the overall mean) and also tau_1,
tau_2 and tau3?

Thanks very much for the help!
   Hanna

> summary(mod1)Linear mixed-effects model fit by REML
 Data: one
   AIC  BIClogLik
  304.4703 330.1879 -140.2352

Random effects:
 Formula: ~1 + time | individual
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 0.2487869075 (Intr)
time0.0001841179 -0.056
Residual0.3718305953

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | method
 Parameter estimates:
   312
 1.0 26.59750 24.74476
Fixed effects: reponse ~ method * time
Value Std.Error DF   t-value p-value(Intercept)
96.65395  3.528586 57 27.391694  0.
method2   1.17851  4.856026 57  0.242689  0.8091
method3   5.87505  3.528617 57  1.664973  0.1014time
0.07010  0.250983 57  0.279301  0.7810
method2:time -0.12616  0.360585 57 -0.349877  0.7277
method3:time -0.08010  0.251105 57 -0.318999  0.7509
 Correlation:
 (Intr) methd2 methd3 time   mthd2:
method2  -0.726
method3  -0.999  0.726
time -0.779  0.566  0.779
method2:time  0.542 -0.712 -0.542 -0.696
method3:time  0.778 -0.566 -0.779 -0.999  0.696

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max
-2.67575293 -0.51633192  0.06742723  0.59706762  2.81061874

Number of Observations: 69
Number of Groups: 7 >





> one   response individual time method
1102.9  30  3
2103.0  33  3
3103.0  36  3
4102.8  39  3
5102.2  3   12  3
6102.5  3   15  3
7103.0  3   18  3
8102.0  3   24  3
9102.8  10  3
10   102.7  13  3
11   103.0  16  3
12   102.2  19  3
13   103.0  1   12  3
14   102.8  1   15  3
15   102.8  1   18  3
16   102.9  1   24  3
17   102.2  20  3
18   102.6  23  3
19   103.4  26  3
20   102.3  29  3
21   101.3  2   12  3
22   102.1  2   15  3
23   102.1  2   18  3
24   102.2  2   24  3
25   102.7  40  3
26   102.3  43  3
27   102.6  46  3
28   102.7  49  3
29   102.8  4   12  3
30   102.5  50  3
31   102.4  53  3
32   102.1  56  3
33   102.3  60  3
34   102.3  63  3
35   101.9  70  3
36   102.0  73  3
37   107.4  30  1
38   101.3  3   12  1
3992.8  3   15  1
4073.7  3   18  1
41   104.7  3   24  1
4292.6  10  1
43   101.9  1   12  1
44   106.3  1   15  1
45   104.1  1  

[R] read multiple sheets of excel data into R

[li li]

Hi all,
  I tried to use the package "XLConnect" to read excel data into R.  I got
the following error message:

Error : .onLoad failed in loadNamespace() for 'rJava', details:
  call: fun(libname, pkgname)
  error: No CurrentVersion entry in Software/JavaSoft registry! Try
re-installing Java and make sure R and Java have matching
architectures.


  I tried read.xls and got the following error  message:

> library(gdata)> one <- read.xls ("one.xlsx", sheet=1)Error in 
> findPerl(verbose = verbose) :
  perl executable not found. Use perl= argument to specify the correct
path.Error in file.exists(tfn) : invalid 'file' argument


  Can anyone give me some input on this?
  Thanks.
Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Error because of large dimension

[li li]

Thanks all for the reply. I think I need to think of other ways to approach
the problem.
Hanna

2016-01-24 17:45 GMT-05:00 William Dunlap :

> > 28 PiB.  Storing such a large matrix even on file is not possible.
>
> The ads for Amazon Red Shift say it is possible.  E.g.,
>   Amazon Redshift is a fast, fully managed, petabyte-scale data
>   warehouse that makes it simple and cost-effective to analyze
>   all your data using your existing business intelligence tools.
>   Start small for $0.25 per hour with no commitments and scale
>   to petabytes for $1,000 per terabyte per year, less than a tenth
>   the cost of traditional solutions. Customers typically see 3x
>   compression, reducing their costs to $333 per uncompressed
>   terabyte per year.
>
> Cost may be an issue:
>
> 28 petabytes * 1024 petabytes/terabyte * $333 terabyte/year ~= $9.5
> million/year
> or $26 thousand/day.
>
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
> On Sun, Jan 24, 2016 at 1:29 PM, Henrik Bengtsson <
> henrik.bengts...@gmail.com> wrote:
>
>> FYI, the matrix you tried to allocate would hold
>> (3195*1290*495*35*35*35*15) * 3 = 3.936248e+15 values.  Each value
>> would occupy 8 bytes of memory (for the double data type).  In other
>> words, in order to keep this data matrix in memory you would require a
>> computer with at least 3.148998e+16 bytes of RAM, i.e. 29327331 GiB =
>> 28640 TiB = 28 PiB.  Storing such a large matrix even on file is not
>> possible.
>>
>> In other words, you need to figure out how to approach your original
>> problem in a different way.
>>
>> /Henrik
>>
>> On Sun, Jan 24, 2016 at 8:46 AM, li li  wrote:
>> > Hi all,
>> >   I am doing some calculation with very large dimension. I need to
>> create a
>> > matrix
>> > with three columns and a very large number of rows
>> > (3195*1290*495*35*35*35*15=1.312083e+15) i
>> > n order to allocate calculation result from a for loop.
>> > R does not allow me to create such a matrix because of the large
>> dimension
>> > (see below). Is there a way to go around this?
>> >   Thanks very much!!
>> >  Hanna
>> >
>> >
>> >> matrix(0, 3195*1290*495*35*35*35*15, 3)
>> > Error in matrix(0, 3195 * 1290 * 495 * 35 * 35 * 35 * 15, 3) :
>> >   invalid 'nrow' value (too large or NA)
>> > In addition: Warning message:
>> > In matrix(0, 3195 * 1290 * 495 * 35 * 35 * 35 * 15, 3) :
>> >   NAs introduced by coercion
>> >>
>> >
>> > [[alternative HTML version deleted]]
>> >
>> > __
>> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> > 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.
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Error because of large dimension

[li li]

Hi all,
  I am doing some calculation with very large dimension. I need to create a
matrix
with three columns and a very large number of rows
(3195*1290*495*35*35*35*15=1.312083e+15) i
n order to allocate calculation result from a for loop.
R does not allow me to create such a matrix because of the large dimension
(see below). Is there a way to go around this?
  Thanks very much!!
 Hanna


> matrix(0, 3195*1290*495*35*35*35*15, 3)
Error in matrix(0, 3195 * 1290 * 495 * 35 * 35 * 35 * 15, 3) :
  invalid 'nrow' value (too large or NA)
In addition: Warning message:
In matrix(0, 3195 * 1290 * 495 * 35 * 35 * 35 * 15, 3) :
  NAs introduced by coercion
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Logical operator in R

[li li]

I see. Thanks!

2016-01-22 10:57 GMT-05:00 Rmh :

> FAQ 7.31
>
> in this case subtract the two numbers and see that
> they differ by about 1e-16
>
> Sent from my iPhone
>
> > On Jan 22, 2016, at 10:46, li li  wrote:
> >
> > Hi all,
> >  I encountered the following strange phenomenon.
> > For some reason, the obs_p[1] and res1$st_p[89] have
> > the same value but when I run "==", it returns FALSE.
> > Can anyone help give some explanation on this?
> >  Thanks very much!
> >Hanna
> >
> >> obs_p[1]
> > [1] 0.002201438
> >> res1$st_p[89]
> > [1] 0.002201438
> >> res1$st_p[89]==obs_p[1]
> > [1] FALSE
> >> res1$st_p[89] > [1] FALSE
> >> res1$st_p[89]>obs_p[1]
> > [1] TRUE
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Logical operator in R

[li li]

Hi all,
  I encountered the following strange phenomenon.
For some reason, the obs_p[1] and res1$st_p[89] have
the same value but when I run "==", it returns FALSE.
Can anyone help give some explanation on this?
  Thanks very much!
Hanna

> obs_p[1]
[1] 0.002201438
> res1$st_p[89]
[1] 0.002201438
> res1$st_p[89]==obs_p[1]
[1] FALSE
> res1$st_p[89] res1$st_p[89]>obs_p[1]
[1] TRUE

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] permTREND function in package "perm"

[li li]

Hi all,
  I am trying to figure out how to use the function permTREND correctly in
package "perm". There does not seem to be any examples given for this
function. In the help file, it says the following:

## Default S3 method:
permTREND(x, y, alternative = c("two.sided", "less", "greater"), exact =
NULL, method = NULL, methodRule = methodRuleTREND1,
control=permControl(),...)

## S3 method for class 'formula'
permTREND(formula,data,subset,na.action,...)

where,

x

numeric vector of response scores for the first group
y

numeric vector of either response scores for the second group (for permTS)
or trend scores for each observation (for permTREND)
g

a factor or character vector denoting group membership

But I am still not sure. What does it mean by S3  method for formula?
also, trend score in the explanation for "y" mean what?

Take the data below for example, x should be the count value for the row
corresponding to "Yes", and "y" should be the dose levels "0, 0.15. 0.5,
1.5, 5"?  Need some help on this. Thanks!
Hanna


> addmargins(dat)

dose 0 dose 0.15 dose 0.5 dose 1.5 dose 5 Sum

yes  4 345  8  24

no   4 543  0  16

Sum  8 888  8  40

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Fwd: exact trend test (enumerate all possible contingency tables with fixed row and column margins)

[li li]

As a follow up, I found out the proc freq in SAS can perform the exact
permutation trend test.

  Hanna

2016-01-08 9:30 GMT-05:00 li li :

> Thanks Bert.
>
>
> 2016-01-07 13:39 GMT-05:00 Bert Gunter :
>
>> Sorry -- neglected to reply to the list. -- Bert
>>
>>
>>
>> -- Forwarded message --
>> From: Bert Gunter 
>> Date: Thu, Jan 7, 2016 at 10:38 AM
>> Subject: Re: [R] exact trend test (enumerate all possible contingency
>> tables with fixed row and column margins)
>> To: li li 
>>
>>
>> I do not know whether there is any package to do what you want.
>>
>> I **do** know that the algorithms required to do this are very
>> sophisticated and that with more than a few groups, all possible
>> enumerations are out of the question so that approximating shortcuts
>> must be used. See http://www.cytel.com/software-solutions/statxact for
>> some background.
>>
>> I **suspect** that you have no need to do what you have requested and
>> **suggest** that you consult a local statistician or
>> stats.stackexchange.com for another approach to whatever your
>> underlying issue is.
>>
>> Cheers,
>> Bert
>> Bert Gunter
>>
>> "The trouble with having an open mind is that people keep coming along
>> and sticking things into it."
>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>>
>>
>> On Thu, Jan 7, 2016 at 10:18 AM, li li  wrote:
>> > I did check the coin package before. I did not see a function in that
>> > package that can be used to list all the possible contingency tables
>> with
>> > fixed margins.
>> > Of course I googled "exact trend test using R". There is not enough help
>> > there.
>> > For up to three groups, I can easily enumerate all the contingency table
>> > with fixed margins, but with 5 groups it is not that easy.
>> > But as mentioned before, this is done implicitly and routinely in
>> > fisher.test function in R. So if anyone who have done this in R before,
>> > please help.
>> > Thanks.
>> >Hanna
>> >
>> >
>> > 2016-01-07 12:20 GMT-05:00 Michael Dewey :
>> >
>> >> You received a number of suggestions about where to look and packages
>> that
>> >> might be suitable. Did you do that? If you did which ones did you look
>> at
>> >> and why did you reject them?
>> >>
>> >>
>> >> On 07/01/2016 16:29, li li wrote:
>> >>
>> >>> Thanks for all the reply. Below is the data in a better format.
>> >>>
>> >>> addmargins(dat)
>> >>>>
>> >>>
>> >>>  dose 0 dose 0.15 dose 0.5 dose 1.5 dose 5 Sum
>> >>>
>> >>> yes  4 345  8  24
>> >>>
>> >>> no   4 543  0  16
>> >>>
>> >>> Sum      8 888  8  40
>> >>>
>> >>> I think it is easier and better that I rephrase my question. I would
>> like
>> >>> to enumerate all possible
>> >>> contingency tables with the row margins and column margins fixed as
>> in the
>> >>> above table. Yes. In fisher's exact test, this should have been done
>> >>> internally. But I need explicitly find all such tables. Need some
>> help on
>> >>> this and thanks very much in advance.
>> >>>
>> >>>  Hanna
>> >>>
>> >>>
>> >>> 2016-01-07 7:15 GMT-05:00 peter dalgaard :
>> >>>
>> >>>
>> >>>> On 07 Jan 2016, at 08:31 , David Winsemius 
>> >>>> wrote:
>> >>>>
>> >>>>
>> >>>>>> On Jan 6, 2016, at 8:16 PM, li li  wrote:
>> >>>>>>
>> >>>>>> Hi all,
>> >>>>>> Is there an R function that does exact randomization trend test?
>> >>>>>> For example, consider the 2 by 5 contingency table below:
>> >>>>>>
>> >>>>>>dose0dose 0.15dose 0.5dose 1.5dose 5
>> >>>>>>  row
>> >>>>>> margin
>> >>>>>> Yes  43  4   5
>> >>>>>> 8   24
>> >>>

Re: [R] Fwd: exact trend test (enumerate all possible contingency tables with fixed row and column margins)

[li li]

Thanks Bert.


2016-01-07 13:39 GMT-05:00 Bert Gunter :

> Sorry -- neglected to reply to the list. -- Bert
>
>
>
> -- Forwarded message --
> From: Bert Gunter 
> Date: Thu, Jan 7, 2016 at 10:38 AM
> Subject: Re: [R] exact trend test (enumerate all possible contingency
> tables with fixed row and column margins)
> To: li li 
>
>
> I do not know whether there is any package to do what you want.
>
> I **do** know that the algorithms required to do this are very
> sophisticated and that with more than a few groups, all possible
> enumerations are out of the question so that approximating shortcuts
> must be used. See http://www.cytel.com/software-solutions/statxact for
> some background.
>
> I **suspect** that you have no need to do what you have requested and
> **suggest** that you consult a local statistician or
> stats.stackexchange.com for another approach to whatever your
> underlying issue is.
>
> Cheers,
> Bert
> Bert Gunter
>
> "The trouble with having an open mind is that people keep coming along
> and sticking things into it."
> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>
>
> On Thu, Jan 7, 2016 at 10:18 AM, li li  wrote:
> > I did check the coin package before. I did not see a function in that
> > package that can be used to list all the possible contingency tables with
> > fixed margins.
> > Of course I googled "exact trend test using R". There is not enough help
> > there.
> > For up to three groups, I can easily enumerate all the contingency table
> > with fixed margins, but with 5 groups it is not that easy.
> > But as mentioned before, this is done implicitly and routinely in
> > fisher.test function in R. So if anyone who have done this in R before,
> > please help.
> > Thanks.
> >Hanna
> >
> >
> > 2016-01-07 12:20 GMT-05:00 Michael Dewey :
> >
> >> You received a number of suggestions about where to look and packages
> that
> >> might be suitable. Did you do that? If you did which ones did you look
> at
> >> and why did you reject them?
> >>
> >>
> >> On 07/01/2016 16:29, li li wrote:
> >>
> >>> Thanks for all the reply. Below is the data in a better format.
> >>>
> >>> addmargins(dat)
> >>>>
> >>>
> >>>  dose 0 dose 0.15 dose 0.5 dose 1.5 dose 5 Sum
> >>>
> >>> yes  4 345  8  24
> >>>
> >>> no   4 543  0  16
> >>>
> >>> Sum  8 888  8  40
> >>>
> >>> I think it is easier and better that I rephrase my question. I would
> like
> >>> to enumerate all possible
> >>> contingency tables with the row margins and column margins fixed as in
> the
> >>> above table. Yes. In fisher's exact test, this should have been done
> >>> internally. But I need explicitly find all such tables. Need some help
> on
> >>> this and thanks very much in advance.
> >>>
> >>>  Hanna
> >>>
> >>>
> >>> 2016-01-07 7:15 GMT-05:00 peter dalgaard :
> >>>
> >>>
> >>>> On 07 Jan 2016, at 08:31 , David Winsemius 
> >>>> wrote:
> >>>>
> >>>>
> >>>>>> On Jan 6, 2016, at 8:16 PM, li li  wrote:
> >>>>>>
> >>>>>> Hi all,
> >>>>>> Is there an R function that does exact randomization trend test?
> >>>>>> For example, consider the 2 by 5 contingency table below:
> >>>>>>
> >>>>>>dose0dose 0.15dose 0.5dose 1.5dose 5
> >>>>>>  row
> >>>>>> margin
> >>>>>> Yes  43  4   5
> >>>>>> 8   24
> >>>>>> No  45   4   3
> >>>>>>   0  16
> >>>>>> col sum88   8   8
> >>>>>>   8   40
> >>>>>>
> >>>>>
> >>>>> Your data presentation has been distorted by your failure to post in
> >>>>>
> >>>> plain text. Surely you have been asked in the past to correct this
> issue?
> >>>>
> >>

Re: [R] exact trend test (enumerate all possible contingency tables with fixed row and column margins)

[li li]

I did check the coin package before. I did not see a function in that
package that can be used to list all the possible contingency tables with
fixed margins.
Of course I googled "exact trend test using R". There is not enough help
there.
For up to three groups, I can easily enumerate all the contingency table
with fixed margins, but with 5 groups it is not that easy.
But as mentioned before, this is done implicitly and routinely in
fisher.test function in R. So if anyone who have done this in R before,
please help.
Thanks.
   Hanna


2016-01-07 12:20 GMT-05:00 Michael Dewey :

> You received a number of suggestions about where to look and packages that
> might be suitable. Did you do that? If you did which ones did you look at
> and why did you reject them?
>
>
> On 07/01/2016 16:29, li li wrote:
>
>> Thanks for all the reply. Below is the data in a better format.
>>
>> addmargins(dat)
>>>
>>
>>  dose 0 dose 0.15 dose 0.5 dose 1.5 dose 5 Sum
>>
>> yes  4 345  8  24
>>
>> no   4 543  0  16
>>
>> Sum  8 888  8  40
>>
>> I think it is easier and better that I rephrase my question. I would like
>> to enumerate all possible
>> contingency tables with the row margins and column margins fixed as in the
>> above table. Yes. In fisher's exact test, this should have been done
>> internally. But I need explicitly find all such tables. Need some help on
>> this and thanks very much in advance.
>>
>>  Hanna
>>
>>
>> 2016-01-07 7:15 GMT-05:00 peter dalgaard :
>>
>>
>>> On 07 Jan 2016, at 08:31 , David Winsemius 
>>> wrote:
>>>
>>>
>>>>> On Jan 6, 2016, at 8:16 PM, li li  wrote:
>>>>>
>>>>> Hi all,
>>>>> Is there an R function that does exact randomization trend test?
>>>>> For example, consider the 2 by 5 contingency table below:
>>>>>
>>>>>dose0dose 0.15dose 0.5dose 1.5dose 5
>>>>>  row
>>>>> margin
>>>>> Yes  43  4   5
>>>>> 8   24
>>>>> No  45   4   3
>>>>>   0  16
>>>>> col sum88   8   8
>>>>>   8   40
>>>>>
>>>>
>>>> Your data presentation has been distorted by your failure to post in
>>>>
>>> plain text. Surely you have been asked in the past to correct this issue?
>>>
>>>>
>>>>
>>>>> To do the exact trend test, we need to enumerate all the contingency
>>>>>
>>>> table
>>>
>>>> with the
>>>>> row and column margins fixed.
>>>>>
>>>>
>>>> Er, how should that be done? A trend test? What is described above would
>>>>
>>> be a general test of no association rather than a trend test. Please use
>>> clear language and be as specific as possible if you choose to respond.
>>>
>>>>
>>>> Find the probability corresponding to
>>>>> obtaining
>>>>> the corresponding contingency tables based on the multivariate
>>>>> hypergeometric distribution. Finally the pvalue is obtained by adding
>>>>> relevant probabilities.
>>>>>
>>>>
>>>> If there is a trend under consideration, then I do not understand such a
>>>>
>>> trend would be modeled under a hypergeometric distribution? A
>>> hypergeometic
>>> distribution would suggest no trend, at least to my current
>>> understanding.
>>>
>>> I'd expect that there is such a beast as a noncentral multivariate
>>> hypergeometric (for the 2x2 case that is what we use to get the CI for
>>> the
>>> odds ratio), but usually, one just wants the null distribution of the
>>> test
>>> statistic.
>>>
>>>
>>>
>>>>
>>>>> Is there an R function that does this? if not, I am wondering whether
>>>>>
>>>> it is
>>>
>>>> possible to
>>>>> enumerate all possible contingency tables that has column sun and row
>>>>>
>>>> sum
>>>
>>>> fixed?
>>>>>
>>>>
>>&g

[R] exact trend test

[li li]

Hi all,
  Is there an R function that does exact randomization trend test?
  For example, consider the 2 by 5 contingency table below:

dose0dose 0.15dose 0.5dose 1.5dose 5   row
margin
 Yes  43  4   5
 8   24
  No  45   4   3
   0  16
col sum88   8   8
   8   40

To do the exact trend test, we need to enumerate all the contingency table
with the
row and column margins fixed. Find the probability corresponding to
obtaining
the corresponding contingency tables based on the multivariate
hypergeometric distribution. Finally the pvalue is obtained by adding
relevant probabilities.

Is there an R function that does this? if not, I am wondering whether it is
possible to
enumerate all possible contingency tables that has column sun and row sum
fixed?

Thanks very much!!

   Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Restricted 4-PL curves using drc package

[li li]

The option 2 works. Thanks.
  Hanna

2015-12-09 16:24 GMT-05:00 Bert Gunter :

> 1. Don't know, but, assuming I understand you correctly, probably not.
>
> 2 Note (subject to my understanding again) that you can rephrase your
> query as fitting a 4 parameter logistic to all the data together with
> a grouping factor to indicate the separate curves, allowing two of the
> parameters (e.g. ED50 and slope) to vary by group. You can then use
> standard nonlinear fitting (e.g. via the optimx package) to fit the
> model. Note that you will have to specify 2c +2 starting values for c
> curves. With c "large" (?whatever that means), this may make
> convergence tricky.
>
> If you care to reply, please do so to the list, not to me.
>
> Cheers,
> Bert
>
>
>
>
> Bert Gunter
>
> "The trouble with having an open mind is that people keep coming along
> and sticking things into it."
> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>
>
> On Wed, Dec 9, 2015 at 12:27 PM, li li  wrote:
> > Hi all,
> >   In drc package, is there a function which can be used to fit restricted
> > 4PL curves? For example, we restrict two 4PL curves have the same lower
> and
> > upper asymptotes?
> >   Thanks for the help in advance!
> > Hanna
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Restricted 4-PL curves using drc package

[li li]

Hi all,
  In drc package, is there a function which can be used to fit restricted
4PL curves? For example, we restrict two 4PL curves have the same lower and
upper asymptotes?
  Thanks for the help in advance!
Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] change the x axis tickmarks when using plot function in drc package

[li li]

Thanks for your help Giorgio! Both options worked.

2015-12-09 11:46 GMT-05:00 Giorgio Garziano :

> Looking at the source code of the package drc, there is something that may
> somehow explain what
> you are experiencing:
>
> file: plot.drc.R, function addAxes(), lines 543-626
>
> ceilingxTicks <- ceiling(log10(xaxisTicks[-1]))
> ...
> xaxisTicks <- c(xaxisTicks[1], 10^(unique(ceilingxTicks)))
> 
> xLabels <- as.character(xaxisTicks)
>
>
> I may suggest two options:
>
>
> 1.  provide the x labels at plot() call time:
>
>
>
>plot(mod, type="all", log="x", xtlab = c(-2, -1, 0, 1, 2),
> xlab="log(dose)")
>
>
> 2.  try to use the option logDose in drm():
>
>   mod <- drm(y~log(dose), fct = LL.4(), logDose=10)
>   plot(mod, type="all")
>
> Not sure if that second option fits your needs.
>
>
> --
> GG
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] change the x axis tickmarks when using plot function in drc package

[li li]

Thanks for the reply but that does not seem to work. With that the plot is
on log scale for both the response and the dose levels, but the tickmarks
are still on the original scale.


2015-12-08 22:22 GMT-05:00 Jim Lemon :

> Hi Hannah,
> Try this:
>
> plot(mod, type="all",log="xy")
>
> Jim
>
>
> On Wed, Dec 9, 2015 at 1:57 PM, li li  wrote:
>
>> Hi all,
>>   When plotting the dose response curve using plot function as in the
>> example codes below, the scale of the axis should really be on the log
>> scale of the dose given the shape of the graph. But as you can see, the
>> tickmarks of the returned graph represent the original scale. How can I
>> change the tick marks to the scale that the correspond to log(dose)?
>>   Thanks much in advance.
>> Hanna
>>
>>
>> dose <- rep(50*2^(-(0:11)),3)
>> dose
>> d <- 100
>> c <- 1
>> b <- 1
>> e <- 1.6
>> y <- rnorm(length(dose))+ c+ (d-c)/(1+exp(b*(log(dose)-log(e
>> library(drc)
>> mod <- drm(y~dose, fct = LL.4())
>> summary(mod)
>> plot(mod, type="all")
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] change the x axis tickmarks when using plot function in drc package

[li li]

Hi all,
  When plotting the dose response curve using plot function as in the
example codes below, the scale of the axis should really be on the log
scale of the dose given the shape of the graph. But as you can see, the
tickmarks of the returned graph represent the original scale. How can I
change the tick marks to the scale that the correspond to log(dose)?
  Thanks much in advance.
Hanna


dose <- rep(50*2^(-(0:11)),3)
dose
d <- 100
c <- 1
b <- 1
e <- 1.6
y <- rnorm(length(dose))+ c+ (d-c)/(1+exp(b*(log(dose)-log(e
library(drc)
mod <- drm(y~dose, fct = LL.4())
summary(mod)
plot(mod, type="all")

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] [FORGED] Error in loading the drc package

[li li]

Thanks.

2015-12-08 15:20 GMT-05:00 Rolf Turner :

>
> Please keep communications on-list.  Others may have relevant comments and
> suggestions to make.
>
>
> On 09/12/15 09:00, li li wrote:
>
>> Thanks for the reply. So a newer version of R can solve the problem? But
>> I was able to successfully load the package yesterday.
>> Thanks.
>>
>
> I have no access to your system so I cannot advise.  Apparently
> *something* changed with your system between sometime yesterday and
> whenever you got the error.  Only you are in a position to know what
> changed.
>
> However that is probably irrelevant.  Just
>
>   * update R
>   * update (re-install) the drc package
>
> and things should work without errors or warnings.
>
> cheers,
>
> Rolf Turner
>
> --
> Technical Editor ANZJS
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Error in loading the drc package

[li li]

Hi all,
  When trying to load the drc package. I got the following error. Any
suggestions?
  Thanks.
 Hanna


> install.packages("drc", dependencies=TRUE)
--- Please select a CRAN mirror for use in this session ---
trying URL 'https://cran.fhcrc.org/bin/windows/contrib/3.1/drc_2.5-12.zip'
Content type 'application/zip' length 502776 bytes (490 Kb)
opened URL
downloaded 490 Kb
package ‘drc’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\Temp\Rtmpysroid\downloaded_packages

> library(drc)
Error in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck
= vI[[j]]) :
  namespace ‘Matrix’ 1.1-3 is already loaded, but >= 1.1.5 is required
In addition: Warning message:
package ‘drc’ was built under R version 3.1.3
Error: package or namespace load failed for ‘drc’
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Plot in drc package

[li li]

Hi all,
  I have the following data and I fit a log logistic model using the drm
function in DRC package.
I saved the fitted model in an object called "mod". (See below)
   I understand that we can use the plot function to plot the regression
curve with the original data points using the code "plot(mod, type="all",
col="blue")". I have two questions regarding the generated plot:

1. the concentration stop at 10 in the graph. How can I show data for
all concentration values? (in the original data set, the largest
concentration is 20).
2. The plot looks like on the log (conc) scale (If I simply plot resp
against conc, the shape of the data is not sigmoid). However, the x label
automatically and the scale of the x axis come out as conc. It seems to be
it should be on the log scale instead.

Thanks very much in advance!!
   Hanna

 > dat
   resp conc
1   -68 20.0
2   -74 20.0
3   -75 20.0
4   -47 10.0
5   -53 10.0
6   -58 10.0
7 8  5.0
8   -25  5.0
9   -18  5.0
10   75  2.5
11   89  2.5
12  116  2.5
13  346  1.25000
14  471  1.25000
15  515  1.25000
16 1003  0.62500
17 1033  0.62500
18  866  0.62500
19 1384  0.31250
20 1431  0.31250
21 1460  0.31250
22 1853  0.15625
23 1918  0.15625
24 1858  0.15625
25 1995  0.078125000
26 2034  0.078125000
27 1991  0.078125000
28 2065  0.039062500
29 2044  0.039062500
30 2090  0.039062500
31 2148  0.019531250
32 2048  0.019531250
33 2003  0.019531250
34 2075  0.009765625
35 2068  0.009765625
36 2056  0.009765625
> mod <- drm(resp ~ conc, fct=LL.4(), data=dat)
> summary(mod)
Model fitted: Log-logistic (ED50 as parameter) (4 parms)
Parameter estimates:
 Estimate  Std. Error t-value p-value
b:(Intercept)1.5367830.060292   25.489167   0e+00
c:(Intercept)  -86.952455   18.941887   -4.590485   1e-04
d:(Intercept) 2089.228320   16.413634  127.286154   0e+00
e:(Intercept)0.5819790.017057   34.119620   0e+00
Residual standard error:
 48.61747 (32 degrees of freedom)

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] question with drm function in "drc package"

[li li]

Thanks. I see what went wrong.


2015-12-07 16:33 GMT-05:00 Jeff Newmiller :

> The fine manual for the drm function mentions a parameter "logDose" ... I
> expect that since you did not specify it that the drm function is helpfully
> attempting to take the log of your data for you.
>
> I highly recommend reading up on making reproducible examples so your
> helpers can experiment with your problem directly instead of guessing by
> reading only. Of course, you would then want to post using plain text as
> the Posting Guide indicates so the HTML would not corrupt your code.
>
>
> http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example
> --
> Sent from my phone. Please excuse my brevity.
>
> On December 7, 2015 1:05:32 PM PST, li li  wrote:
>>
>> Hi Jeff,
>>   I am not sure where the log of 0 was taken?
>>   Thanks.
>> Li
>>
>> 2015-12-07 15:55 GMT-05:00 Jeff Newmiller :
>>
>>> Log of zero is infinity. Don't do that.
>>> --
>>> Sent from my phone. Please excuse my brevity.
>>>
>>> On December 7, 2015 12:15:48 PM PST, li li  wrote:
>>>
>>>> Hi all,
>>>>   I am trying to use the drm function in drc package to fit a 4 PL or 3PL
>>>> curve for an assay response. Please see the listed data below. When I do
>>>> the curve fitting, it returns the following error message. Anyone who
>>>> familiar with this have any input on what went wrong?
>>>>Thanks so much in advance!!
>>>>   Hanna
>>>>
>>>>  dat
>>>>>
>>>>   values   log_dose
>>>> 1   -68.1125  3.9120230
>>>> 2   -73.8955  3.9120230
>>>> 3   -75.1235  3.9120230
>>>> 4   -47.3205  3.2188758
>>>> 5   -52.9835  3.2188758
>>>> 6   -58.0075  3.2188758
>>>> 7 8.2515  2.5257286
>>>> 8   -24.5925  2.5257286
>>>> 9   -18.1425  2.5257286
>>>> 10   75.3375  1.8325815
>>>> 11   89.3755  1.8325815
>>>> 12  115.9685  1.8325815
>>>> 13  345.7675  1.1394343
>>>> 14  470.8125
>>>> 1.1394343
>>>> 15  514.8985  1.1394343
>>>> 16
>>>> 1003.2235  0.4462871
>>>> 17 1033.0345  0.4462871
>>>> 18  866.1365  0.4462871
>>>> 19 1383.6525 -0.2468601
>>>> 20 1431.3245 -0.2468601
>>>> 21 1459.8025 -0.2468601
>>>> 22 1852.5795 -0.9400073
>>>> 23 1917.8015 -0.9400073
>>>> 24 1858.0875 -0.9400073
>>>> 25 1995.1185 -1.6331544
>>>> 26 2033.9455 -1.6331544
>>>> 27 1991.0405 -1.6331544
>>>> 28 2064.5855 -2.3263016
>>>> 29 2043.5195 -2.3263016
>>>> 30 2089.8525 -2.3263016
>>>> 31 2147.8445 -3.0194488
>>>> 32 2047.7905 -3.0194488
>>>> 33 2002.5375 -3.0194488
>>>> 34 2075.2665 -3.7125960
>>>> 35 2068.3545 -3.7125960
>>>> 36 2055.9605 -3.7125960
>>>>
>>>>>  mod1 <- drm(values~log_dose, fct=LL.4(), data=dat)
>>>>>
>>>> Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control
>>>> = list(maxit = maxIt,  :
>>>>   initial value in 'vmmin' is not
>>>> finite
>>>> Error in drmOpt(opfct, opdfct1, startVecSc,
>>>> optMethod, constrained,
>>>> warnVal,  :
>>>>   Convergence failed
>>>>
>>>>>
>>>>>
>>>>  [[alternative HTML version deleted]]
>>>>
>>>> --
>>>>
>>>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>> 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.
>>>>
>>>>
>>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] question with drm function in "drc package"

[li li]

Hi all,
  I am trying to use the drm function in drc package to fit a 4 PL or 3PL
curve for an assay response. Please see the listed data below. When I do
the curve fitting, it returns the following error message. Anyone who
familiar with this have any input on what went wrong?
   Thanks so much in advance!!
  Hanna

> dat
  values   log_dose
1   -68.1125  3.9120230
2   -73.8955  3.9120230
3   -75.1235  3.9120230
4   -47.3205  3.2188758
5   -52.9835  3.2188758
6   -58.0075  3.2188758
7 8.2515  2.5257286
8   -24.5925  2.5257286
9   -18.1425  2.5257286
10   75.3375  1.8325815
11   89.3755  1.8325815
12  115.9685  1.8325815
13  345.7675  1.1394343
14  470.8125  1.1394343
15  514.8985  1.1394343
16 1003.2235  0.4462871
17 1033.0345  0.4462871
18  866.1365  0.4462871
19 1383.6525 -0.2468601
20 1431.3245 -0.2468601
21 1459.8025 -0.2468601
22 1852.5795 -0.9400073
23 1917.8015 -0.9400073
24 1858.0875 -0.9400073
25 1995.1185 -1.6331544
26 2033.9455 -1.6331544
27 1991.0405 -1.6331544
28 2064.5855 -2.3263016
29 2043.5195 -2.3263016
30 2089.8525 -2.3263016
31 2147.8445 -3.0194488
32 2047.7905 -3.0194488
33 2002.5375 -3.0194488
34 2075.2665 -3.7125960
35 2068.3545 -3.7125960
36 2055.9605 -3.7125960
> mod1 <- drm(values~log_dose, fct=LL.4(), data=dat)
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control
= list(maxit = maxIt,  :
  initial value in 'vmmin' is not finite
Error in drmOpt(opfct, opdfct1, startVecSc, optMethod, constrained,
warnVal,  :
  Convergence failed
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] [FORGED] Re: correlated binomial random variables

[li li]

Thanks Dennis and Rolf. Yes. Simulation is one way. I think
correlation does not determine the joint distribution so it will not
be unique. Under specific settings, the joint probability of X, Y can
be calculated. For example, let X=X_0+X_1 and Y=X_0+X_2, with X_0
being Binomial(n_0, p) and X_1, and X_2 are both Binomial(n, p). X_0,
X_1, and X_2 are all independent. Then X, Y are correlated and P(X <=
t, Y <= t) can be exactly calculated.

Thanks!
   Hanna

2015-10-05 18:00 GMT-04:00, Rolf Turner :
> On 06/10/15 04:43, li li wrote:
>> Hi all,
>> Using the "bindata" package, it is possible to gerenerate
>> correlated binomial random variables both with the same number of
>> trials, say n. I am wondering whether there is an R function to
>> calculate the joint probability distribution of the correlated
>> binomial random variables. Say if X is binomial (n, p1) and Y is
>> binomial (n, p2) and the correlation between X and Y is rho and we
>> want to calculate
>> P(X <= c, Y <= c).
>
> (1) The use of correlation in the context of binary or binomial variates
> makes little or no sense, it seems to me.  Correlation is basically
> useful for quantifying linear relationships between continuous variates.
> Linear relationships between count variates are of at best limited
> interest.
>
> (2) I suspect that the correlation does not determine a unique joint
> distribution of X and Y.  If my suspicion is correct then there is not a
> unique (well-defined) answer to the question "What is
> Pr(X <= x, Y <= y)?"
>
> cheers,
>
> Rolf Turner
>
> --
> Technical Editor ANZJS
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] correlated binomial random variables

[li li]

Hi all,
   Using the "bindata" package, it is possible to gerenerate
correlated binomial random variables both with the same number of
trials, say n. I am wondering whether there is an R function to
calculate the joint probability distribution of the correlated
binomial random variables. Say if X is binomial (n, p1) and Y is
binomial (n, p2) and the correlation between X and Y is rho and we
want to calculate
P(X <= c, Y <= c).
   Thanks so much!
 Hanna

2015-10-05 11:27 GMT-04:00, li li :
>
>

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] correlated binomial random variables

[li li]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] cumulative distribtuion function for multinomial distribution in R

[li li]

Hi all,
  In R, is there a function for the cumulative distribution function
for multinomial distribution? I only see pmultinom and rmultinom which
are the prabability mass function and the function for generating
multinomial random variables respectively.
  Thanks!
Hanna

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Order of boxplots

[li li]

Thank you! That worked.


2015-09-15 2:50 GMT-04:00 peter dalgaard :

>
> > On 15 Sep 2015, at 04:31 , li li  wrote:
> >
> > Hi Jeff,
> >  Thanks for replying. I actually tried "ordered(tmp$type, levels=c("c",
> > "b", "a")."
> > But I think only the order of the letters on x axis changed but the order
>
> You _think_ ??? Documentation, please...
>
> The boxplots certainly move if I do
>
> plot(result ~ type, tmp)
> plot(result ~ factor(type, levels=c("c","b","a")), tmp)
>
>
>
> > of the boxplot did not. So there is some problem there. I also tried
> > as.factor(tmp$type); levels(tmp$type)=c("c", "b", "a") and got the same
>
> That changes the level _names_: 1st group name becomes "c" instead of "a";
> you want the 3rd group to become the 1st but still be called "c".
>
> -pd
>
>
> > thing.
> >Thanks.
> >  Li
> >
> > 2015-09-14 21:44 GMT-04:00 Jeff Newmiller :
> >
> >> Make your factor variable deliberately. That is, specify the levels
> >> parameter with the values in order when you create the factor.
> >>
> ---
> >> Jeff NewmillerThe .   .  Go
> Live...
> >> DCN:Basics: ##.#.   ##.#.  Live
> >> Go...
> >>  Live:   OO#.. Dead: OO#..  Playing
> >> Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
> >> /Software/Embedded Controllers)   .OO#.   .OO#.
> rocks...1k
> >>
> ---
> >> Sent from my phone. Please excuse my brevity.
> >>
> >> On September 14, 2015 6:23:50 PM PDT, li li 
> wrote:
> >>> Hi all,
> >>> I have the following data "tmp" and I want to plot  boxplots for
> >>> each level of the factor "type" and the order the factor should be c,
> >>> b ,a. In other words, the boxplot corresponding to the level "c"
> >>> should be the first and so on.
> >>> Any suggestions?
> >>>  Li
> >>>
> >>>> tmp
> >>>  result type
> >>> 1 101a
> >>> 2 101a
> >>> 3 101a
> >>> 4 101a
> >>> 5 101a
> >>> 6 101a
> >>> 7 100a
> >>> 8 106b
> >>> 9  91b
> >>> 10 78b
> >>> 11 95b
> >>> 12111b
> >>> 13 92b
> >>> 14 98b
> >>> 15108c
> >>> 16112c
> >>> 17 98c
> >>> 18102c
> >>> 19 88c
> >>> 20 86c
> >>> 21 81c
> >>>
> >>> __
> >>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >>> 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.
> >>
> >>
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd@cbs.dk  Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Order of boxplots

[li li]

Hi Jeff,
  Thanks for replying. I actually tried "ordered(tmp$type, levels=c("c",
"b", "a")."
But I think only the order of the letters on x axis changed but the order
of the boxplot did not. So there is some problem there. I also tried
as.factor(tmp$type); levels(tmp$type)=c("c", "b", "a") and got the same
thing.
Thanks.
  Li

2015-09-14 21:44 GMT-04:00 Jeff Newmiller :

> Make your factor variable deliberately. That is, specify the levels
> parameter with the values in order when you create the factor.
> ---
> Jeff NewmillerThe .   .  Go Live...
> DCN:Basics: ##.#.   ##.#.  Live
> Go...
>   Live:   OO#.. Dead: OO#..  Playing
> Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
> /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
> -------
> Sent from my phone. Please excuse my brevity.
>
> On September 14, 2015 6:23:50 PM PDT, li li  wrote:
> >Hi all,
> >  I have the following data "tmp" and I want to plot  boxplots for
> >each level of the factor "type" and the order the factor should be c,
> >b ,a. In other words, the boxplot corresponding to the level "c"
> >should be the first and so on.
> >Any suggestions?
> >   Li
> >
> >> tmp
> >   result type
> >1 101a
> >2 101a
> >3 101a
> >4 101a
> >5 101a
> >6 101a
> >7 100a
> >8 106b
> >9  91b
> >10 78b
> >11 95b
> >12111b
> >13 92b
> >14 98b
> >15108c
> >16112c
> >17 98c
> >18102c
> >19 88c
> >20 86c
> >21 81c
> >
> >__
> >R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >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.
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Order of boxplots

[li li]

Hi all,
  I have the following data "tmp" and I want to plot  boxplots for
each level of the factor "type" and the order the factor should be c,
b ,a. In other words, the boxplot corresponding to the level "c"
should be the first and so on.
Any suggestions?
   Li

> tmp
   result type
1 101a
2 101a
3 101a
4 101a
5 101a
6 101a
7 100a
8 106b
9  91b
10 78b
11 95b
12111b
13 92b
14 98b
15108c
16112c
17 98c
18102c
19 88c
20 86c
21 81c

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Adding a second Y axis on a dotplot

[li li]

Hi all,
  I plotted a dotplot based on the data below and code below. I would
like to add another yaxis on the right with a different col, different
tickmarks and a different label. Can anyone give some help?Thanks very
much!!
 Hanna


> tmp1
   result  lot  trt trtsymb trtcol
1  98 lot1 trt1   1   blue
2  99 lot2 trt1   1   blue
3  98 lot3 trt1   1   blue
4 100 lot4 trt1   1   blue
5 100 lot5 trt1   1   blue
6 101 lot6 trt1   1   blue
7 101 lot7 trt1   1   blue
8  99 lot8 trt1   1   blue
9 100 lot9 trt1   1   blue
10 94 lot1 trt2  16red
11105 lot2 trt2  16red
12 87 lot3 trt2  16red
13119 lot4 trt2  16red
14 96 lot5 trt2  16red
15113 lot6 trt2  16red
16106 lot7 trt2  16red
17 71 lot8 trt2  16red
18 95 lot9 trt2  16red


library(lattice)
dotplot(result ~ lot, tmp1, cex=1.1,  ylab = "values", xlab="lot",
jitter.y = F, aspect=1.0,
pch=tmp1$trtsymb, col=tmp1$trtcol, scales=list(rot=30),
main="", key = list(text = list(labels = c("trt1", "trt2"),cex=c(0.9,0.9)),
points = list(pch =c(1,12), col =c("blue", "red")),
space = "right"))

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Different random intercepts but same random slope for groups

[li li]

Thanks all for the reply,

I think what Bert specified is what I wanted. Thanks very much.
So this model allows different random intercept term but the same
random slope term for the three methods.

I have an additional question. I would like to require differnt
residual variance also for the three groups. Is that possible?

Thanks!!

2015-06-09 17:25 GMT-04:00, Bert Gunter :
> Thierry:
>
> I don't think so. It looks to me like her syntax/understanding is confused.
> I think the call should be:
>
> mod2 <- lmer(result  ~ group*time+(group + time|lot), na.action=na.omit,
> data=alldata)
>
> Her request for "the same random slope for each group" -- I assume it's for
> time -- means to me that the time slope will vary "randomly" by lot only,
> the slope would be the same for all groups within the lot.
>
> Of course, I may be wrong also. If so, I suggest that she follow the
> posting guide and post at least head(alldata) using dput() to enable folks
> to understand the structure of her data. And only on r-sig-mixed-models --
> crossposting is frowned upon here and the mixed models list is the best bet
> for this sort of question anyway.
>
> As always, corrections and criticism welcome.
>
> Cheers,
> Bert
>
> Bert Gunter
>
> "Data is not information. Information is not knowledge. And knowledge is
> certainly not wisdom."
>-- Clifford Stoll
>
> On Tue, Jun 9, 2015 at 1:49 PM, Thierry Onkelinx 
> wrote:
>
>> Your model is too complex for the data. This gives you two options: a)
>> simplify the model and b) get more data.
>>
>> Best regards,
>>
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>> and
>> Forest
>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>> Kliniekstraat 25
>> 1070 Anderlecht
>> Belgium
>>
>> To call in the statistician after the experiment is done may be no more
>> than asking him to perform a post-mortem examination: he may be able to
>> say
>> what the experiment died of. ~ Sir Ronald Aylmer Fisher
>> The plural of anecdote is not data. ~ Roger Brinner
>> The combination of some data and an aching desire for an answer does not
>> ensure that a reasonable answer can be extracted from a given body of
>> data.
>> ~ John Tukey
>>
>> 2015-06-09 21:57 GMT+02:00 li li :
>>
>> > Hi all,
>> >   I'd like to fit a random intercept and random slope model. In my
>> > data, there are three groups. I want to have different random
>> > intercept for each group but the same random slope effect for all
>> > three groups. I used the following R command.
>> > However, there seems to be some problem. Any suggestions?
>> >
>> >
>> >
>> > mod2 <- lmer(result  ~ group*time+(0+group1+ group2 +
>> > group3+time|lot), na.action=na.omit, data=alldata)
>> >
>> > > summary(mod2)
>> > Model is not identifiable...
>> > summary from lme4 is returned
>> > some computational error has occurred in lmerTest
>> > Linear mixed model fit by REML ['merModLmerTest']
>> > Formula: result ~ group * time + (0 + group1 + group2 + group3 + time |
>> > lot)
>> >Data: alldata
>> >
>> > REML criterion at convergence: 807.9
>> >
>> > Scaled residuals:
>> > Min  1Q  Median  3Q Max
>> > -3.0112 -0.3364  0.0425  0.2903  3.2017
>> >
>> > Random effects:
>> >  Groups   Name Variance Std.Dev. Corr
>> >  lot  group1   0.0 0.000
>> >   group2   86.20156 9.284  NaN
>> >   group3 55.91479 7.478  NaN  0.06
>> >   time  0.02855 0.169  NaN -0.99  0.10
>> >  Residual  39.91968 6.318
>> > Number of obs: 119, groups:  lot, 15
>> >
>> > Fixed effects:
>> > Estimate Std. Error t value
>> > (Intercept) 100.1566 2.5108   39.89
>> > group  group2-2.9707 3.7490   -0.79
>> > group  group3   -0.0717 2.8144   -0.03
>> > time -0.1346 0.1780   -0.76
>> > group  group2 :time   0.1450 0.29390.49
>> > group  group3:time0.1663 0.21520.77
>> >
>> > Warning messages:
>> > 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
>> :
>> >   Model failed to converge with max|grad| = 0.147314 (tol = 0.002,
>> > comp

[R] Different random intercepts but same random slope for groups

[li li]

Hi all,
  I'd like to fit a random intercept and random slope model. In my
data, there are three groups. I want to have different random
intercept for each group but the same random slope effect for all
three groups. I used the following R command.
However, there seems to be some problem. Any suggestions?



mod2 <- lmer(result  ~ group*time+(0+group1+ group2 +
group3+time|lot), na.action=na.omit, data=alldata)

> summary(mod2)
Model is not identifiable...
summary from lme4 is returned
some computational error has occurred in lmerTest
Linear mixed model fit by REML ['merModLmerTest']
Formula: result ~ group * time + (0 + group1 + group2 + group3 + time |
lot)
   Data: alldata

REML criterion at convergence: 807.9

Scaled residuals:
Min  1Q  Median  3Q Max
-3.0112 -0.3364  0.0425  0.2903  3.2017

Random effects:
 Groups   Name Variance Std.Dev. Corr
 lot  group1   0.0 0.000
  group2   86.20156 9.284  NaN
  group3 55.91479 7.478  NaN  0.06
  time  0.02855 0.169  NaN -0.99  0.10
 Residual  39.91968 6.318
Number of obs: 119, groups:  lot, 15

Fixed effects:
Estimate Std. Error t value
(Intercept) 100.1566 2.5108   39.89
group  group2-2.9707 3.7490   -0.79
group  group3   -0.0717 2.8144   -0.03
time -0.1346 0.1780   -0.76
group  group2 :time   0.1450 0.29390.49
group  group3:time0.1663 0.21520.77

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.147314 (tol = 0.002, component 2)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 2 negative eigenvalues

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Warning message when using lmer function

[li li]

I got the following warning message when using the lmer function.
Does anyone know what is the implication? Thanks!

Warning message:
In anova(model, ddf = "lme4") : bytecode version mismatch; using eval

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Confidence interval for the mean in the random intercept and random slople model

[li li]

Hi all,
  I am fitting a random slope and random intercept model usign lme
fucntion as shown below. Type is factor with two levels. I would like
to to find a confidence interval for mean of this model. Note that the
variance we use in finding the confidence interval should include the
variariance component from random effect.  Any suggestions?


## using lme function
> mod_lme <- lme(ti  ~ type*months, random=~ 1+months|lot, na.action=na.omit,
+ data=one, control = lmeControl(opt = "optim"))
> summary(mod_lme)
Linear mixed-effects model fit by REML
 Data: one
AIC   BIC   logLik
  -82.60042 -70.15763 49.30021

Random effects:
 Formula: ~1 + months | lot
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 8.907584e-03 (Intr)
months  6.039781e-05 -0.096
Residual4.471243e-02

Fixed effects: ti ~ type * months
 Value   Std.Error DF   t-value p-value
(Intercept) 0.25831245 0.016891587 31 15.292373  0.
type0.13502089 0.026676101  4  5.061493  0.0072
months  0.00804790 0.001218941 31  6.602368  0.
type:months -0.00693679 0.002981859 31 -2.326329  0.0267
 Correlation:
   (Intr) typPPQ months
type   -0.633
months -0.785  0.497
type:months  0.321 -0.762 -0.409

Standardized Within-Group Residuals:
  MinQ1   MedQ3   Max
-2.162856e+00 -1.962972e-01 -2.771184e-05  3.749035e-01  2.088392e+00

Number of Observations: 39
Number of Groups: 6

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] How to add legend to a dotplot

[li li]

Thank you! It worked.
   Hanna

2015-06-03 3:49 GMT-04:00, Duncan Mackay :
> Hi
> If you use the lattice package and not resort to other packages read the ?
> xyplot carefully
> There are several ways to put a key in the graph including auto.key if you
> want something simple ie default
> The examples at the bottom give some ideas
>
> If you want something special you can use key = list(...) eg
>
>key  = list(text = list(labels = LETTERS[1:3],   cex  =
> c(1,1,1) ),
>title = "Farm", cex.title = 1,
>points = list(pch = c(20,20,3),
>  col = c("black","grey","grey")),
>lines = list(lwd = c(2,1,2),
> lty = c(1,1,1),
> col = c("black","black","grey80")),
> space = "bottom"),
> which gives 3 columns: text points and lines.
>
> Regards
>
> Duncan
>
> Duncan Mackay
> Department of Agronomy and Soil Science
> University of New England
> Armidale NSW 2351
> Email: home: mac...@northnet.com.au
>
>
>
>
> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of li li
> Sent: Wednesday, 3 June 2015 13:07
> To: r-help
> Subject: [R] How to add legend to a dotplot
>
> Hi all,
>   I wanted to add the legend to a dotplot using legend funciton. If
> does not seem to be working? Anyone have any suggestions?
>   Thanks!
>   Hanna
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] How to add legend to a dotplot

[li li]

Hi all,
  I wanted to add the legend to a dotplot using legend funciton. If
does not seem to be working? Anyone have any suggestions?
  Thanks!
  Hanna

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] different results from lme and lmer function

[li li]

Hi all,
  I am fitting a random slope and random intercept model using R. I
used both lme and lmer funciton for the same model. However I got
different results as shown below (different variance component
estimates and so on). I think that is really confusing. They should
produce close results. Anyone has any thoughts or suggestions. Also,
which one should be comparable to sas results?
 Thanks!
  Hanna

## using lme function
> mod_lme <- lme(ti  ~ type*months, random=~ 1+months|lot, na.action=na.omit,
+ data=one, control = lmeControl(opt = "optim"))
> summary(mod_lme)
Linear mixed-effects model fit by REML
 Data: one
AIC   BIC   logLik
  -82.60042 -70.15763 49.30021

Random effects:
 Formula: ~1 + months | lot
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 8.907584e-03 (Intr)
months  6.039781e-05 -0.096
Residual4.471243e-02

Fixed effects: ti ~ type * months
 Value   Std.Error DF   t-value p-value
(Intercept) 0.25831245 0.016891587 31 15.292373  0.
type0.13502089 0.026676101  4  5.061493  0.0072
months  0.00804790 0.001218941 31  6.602368  0.
type:months -0.00693679 0.002981859 31 -2.326329  0.0267
 Correlation:
   (Intr) typPPQ months
type   -0.633
months -0.785  0.497
type:months  0.321 -0.762 -0.409

Standardized Within-Group Residuals:
  MinQ1   MedQ3   Max
-2.162856e+00 -1.962972e-01 -2.771184e-05  3.749035e-01  2.088392e+00

Number of Observations: 39
Number of Groups: 6




###Using lmer function
> mod_lmer <-lmer(ti  ~ type*months+(1+months|lot), na.action=na.omit, data=one)
> summary(mod_lmer)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to
  degrees of freedom [merModLmerTest]
Formula: ti ~ type * months + (1 + months | lot)
   Data: one

REML criterion at convergence: -98.8

Scaled residuals:
Min  1Q  Median  3Q Max
-2.1347 -0.2156 -0.0067  0.3615  2.0840

Random effects:
 Groups   NameVariance  Std.Dev.  Corr
 lot  (Intercept) 2.870e-04 0.0169424
  months  4.135e-07 0.0006431 -1.00
 Residual 1.950e-03 0.0441644
Number of obs: 39, groups:  lot, 6

Fixed effects:
Estimate Std. Errordf t value Pr(>|t|)
(Intercept) 0.258312   0.018661  4.82  13.842 4.59e-05 ***
type 0.135021   0.028880  6.802000   4.675  0.00245 **
months  0.008048   0.001259 11.943000   6.390 3.53e-05 ***
type:months -0.006937   0.002991 28.91  -2.319  0.02767 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) typPPQ months
type -0.646
months  -0.825  0.533
type:month  0.347 -0.768 -0.421

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] [R-sig-ME] lme function to obtain pvalue for fixed effect

[li li]

Thanks so much for replying.
Yes LimerTest package could be used to get pvalues when using lmer
function. But still the summary and anova function give different
pvalues.
   Hanna

2015-05-26 15:19 GMT-04:00, byron vinueza :
> You can use the lmerTest package .
>
>
>
>
>
> Enviado desde mi iPhone
>
>> El 26/5/2015, a las 13:18, li li  escribió:
>>
>> Hi all,
>>  I am using the lme function to run a random coefficient model. Please
>> see
>> output (mod1) as below.
>>  I need to obtain the pvalue for the fixed effect. As you can see,
>> the pvalues given using the summary function is different from the
>> resutls given in anova function.
>> Why should they be different and which one is the correct one to use?
>>   Thanks!
>>  Hanna
>>
>>
>>> summary(mod1)
>> Linear mixed-effects model fit by REML
>> Data: minus20C1
>>AIC   BIC   logLik
>>  -82.60042 -70.15763 49.30021
>>
>> Random effects:
>> Formula: ~1 + months | lot
>> Structure: General positive-definite, Log-Cholesky parametrization
>>StdDev   Corr
>> (Intercept) 8.907584e-03 (Intr)
>> months  6.039781e-05 -0.096
>> Residual4.471243e-02
>>
>> Fixed effects: ti ~ type * months
>> Value   Std.Error DF   t-value p-value
>> (Intercept) 0.25831245 0.016891587 31 15.292373  0.
>> type 0.13502089 0.026676101  4  5.061493  0.0072
>> months  0.00804790 0.001218941 31  6.602368  0.
>> type:months -0.00693679 0.002981859 31 -2.326329  0.0267
>> Correlation:
>>   (Intr) typ months
>> type-0.633
>> months -0.785  0.497
>> type:months  0.321 -0.762 -0.409
>>
>> Standardized Within-Group Residuals:
>>  MinQ1   MedQ3   Max
>> -2.162856e+00 -1.962972e-01 -2.771184e-05  3.749035e-01  2.088392e+00
>>
>> Number of Observations: 39
>> Number of Groups: 6
>>> anova(mod1)
>>numDF denDF   F-value p-value
>> (Intercept) 131 2084.0265  <.0001
>> type1 4   10.8957  0.0299
>> months  131   38.3462  <.0001
>> type:months 1315.4118  0.0267
>>
>> ___
>> r-sig-mixed-mod...@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] lme function to obtain pvalue for fixed effect

[li li]

Hi all,
  I am using the lme function to run a random coefficient model. Please see
output (mod1) as below.
  I need to obtain the pvalue for the fixed effect. As you can see,
the pvalues given using the summary function is different from the
resutls given in anova function.
Why should they be different and which one is the correct one to use?
   Thanks!
  Hanna


> summary(mod1)
Linear mixed-effects model fit by REML
 Data: minus20C1
AIC   BIC   logLik
  -82.60042 -70.15763 49.30021

Random effects:
 Formula: ~1 + months | lot
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 8.907584e-03 (Intr)
months  6.039781e-05 -0.096
Residual4.471243e-02

Fixed effects: ti ~ type * months
 Value   Std.Error DF   t-value p-value
(Intercept) 0.25831245 0.016891587 31 15.292373  0.
type 0.13502089 0.026676101  4  5.061493  0.0072
months  0.00804790 0.001218941 31  6.602368  0.
type:months -0.00693679 0.002981859 31 -2.326329  0.0267
 Correlation:
   (Intr) typ months
type-0.633
months -0.785  0.497
type:months  0.321 -0.762 -0.409

Standardized Within-Group Residuals:
  MinQ1   MedQ3   Max
-2.162856e+00 -1.962972e-01 -2.771184e-05  3.749035e-01  2.088392e+00

Number of Observations: 39
Number of Groups: 6
> anova(mod1)
numDF denDF   F-value p-value
(Intercept) 131 2084.0265  <.0001
type1 4   10.8957  0.0299
months  131   38.3462  <.0001
type:months 1315.4118  0.0267

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] dotplot

[li li]

Hi all,
  I wrote the following code and have two questions:
  (1) As you can see, I would like different colors for different
types. It does not come out that way in the graph from this code.
Anyone know how to fix this?
   (2) How do I made the lots number on x axis eligible to read?
   (3) How do I remove the verticle line in the graph.
Thanks for your help.
  Hanna


lot <- as.character(c("D130" ,"D135" ,"D1300010", "D1300015",
"D1300020",
 "D1300025" ,"D1300030" ,"D1300035", "D1300040",  "D1300045",
 "D130" ,"D135" ,"D1300010", "D1300015",  "D1300020",
 "D1300025" ,"D1300030" ,"D1300035", "D1300040",  "D1300045",
 "D130" ,"D135" ,"D1300010", "D1300015",  "D1300020",
 "D1300025" ,"D1300030" ,"D1300035", "D1300040",  "D1300045"))

val <- rnorm(30)
type <- rep(1:3,each=10)
tmp <- as.data.frame(cbind(val, type, lot))
tmp$sybm <- rep(c(1, 12, 16), each=10)
tmp$color <- rep(c("blue", "red", "green"), each=10)
tmp$val <- as.numeric(tmp$val)
dotplot(val ~ lot, tmp,
ylab = "values", xlab="lot", aspect=1,
pch=tmp$sybm, color=tmp$color)

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] predict function in regression analysis

[li li]

Hi all,
  I have the following data in which there is one factor lot with six
levels and one continuous convariate time.
I want to fit an Ancova model with common slope and different intercept. So
the six lots will have seperate paralell
regression lines.I wanted to find the upper 95% confidence limit for the
mean of the each of
the regression lines. It doesnot seem straightforward to achieve this using
predict function. Can anyone give some suggestions?

  Here is my data. I only show the first 3 lots. Also I show the model I
used in the end. Thanks very much!
 Hanna

  y lot time
 [1,] 4.5   10
 [2,] 4.5   13
 [3,] 4.7   16
 [4,] 6.7   19
 [5,] 6.0   1   12
 [6,] 4.4   1   15
 [7,] 4.1   1   18
 [8,] 5.3   1   24
 [9,] 4.0   20
[10,] 4.2   23
[11,] 4.1   26
[12,] 6.4   29
[13,] 5.5   2   12
[14,] 3.5   2   15
[15,] 4.6   2   18
[16,] 4.1   2   24
[17,] 4.6   30
[18,] 5.0   33
[19,] 6.2   36
[20,] 5.9   39
[21,] 3.9   3   12
[22,] 5.3   3   15
[23,] 6.9   3   18
[24,] 5.7   3   24


> mod <- lm(y ~ lot+time)
> summary(mod)
Call:
lm(formula = y ~ lot + time)
Residuals:
Min  1Q  Median  3Q Max
-1.5666 -0.3344 -0.1343  0.4479  1.8985
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.743730.36617  12.955 2.84e-14 ***
lot2-0.475000.41129  -1.155   0.2567
lot3 0.412500.41129   1.003   0.3234
lot4 0.961090.47943   2.005   0.0535 .
lot5 0.981090.47943   2.046   0.0490 *
lot6-0.098910.47943  -0.206   0.8379
time 0.025860.02046   1.264   0.2153
---

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Random effect in ancova model

[li li]

Hi all,
  I have the following data and would like to consider the following model:
  value = type + batch (type) + beta1* month +beta2*type*month+ error
Here type is fixed effect and batch is a random effect neste within type,
and month is continuous.
Since there is a random effect in this model. I am not too sure about the
right syntax in R.
Can anyone give some suggestions?
  Thank you.
Hanna




 > mydata
  value type  batch month
 [1,]   0.21 178380 0
 [2,]   0.31 178380 3
 [3,]   0.31 178380 6
 [4,]   0.31 178380 9
 [5,]   0.31 17838012
 [6,]   0.31 17838015
 [7,]   0.41 17838018
 [8,]   0.51 17838024
 [9,]   0.21 189617 0
[10,]   0.31 189617 3
[11,]   0.41 189617 6
[12,]   0.41 189617 9
[13,]   0.41 18961712
[14,]   0.41 18961715
[15,]   0.41 18961718
[16,]   0.41 18961724
[17,]   0.21 197367 0
[18,]   0.31 197367 3
[19,]   0.31 197367 6
[20,]   0.41 197367 9
[21,]   0.41 19736712
[22,]   0.41 19736715
[23,]   0.41 19736718
[24,]   0.41 19736724
[25,]   0.42 286571 0
[26,]   0.42 286571 3
[27,]   0.32 286571 6
[28,]   0.42 286571 9
[29,]   0.42 28657112
[30,]NA2 28657115
[31,]   0.42 297299 0
[32,]   0.42 297299 3
[33,]   0.42 297299 6
[34,]   0.42 297299 9
[35,]   0.42 29729912
[36,]NA2 29729915

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Question with uniroot function

[li li]

Thank you.



2015-04-16 12:33 GMT-04:00 William Dunlap :

> Use optimize() to find the minimum and feed that value into uniroot().
>
>  Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
>  On Thu, Apr 16, 2015 at 7:47 AM, li li  wrote:
>
>>  Hi Jeff,
>>   Thanks for the reply. I am aware that the sign needs to be different at
>> the ends of the starting interval.
>>
>>Another question:
>>
>> Is there a way to set the right end point ( (the "upper" argument in the
>> uniroot function below) as the point where the function takes on its
>> minimun, for example my function f1 below?
>>
>> Thanks very much!
>>
>>
>>
>> u1 <- -3
>> u2 <- 4
>> pi0 <- 0.8
>>
>> f1 <- function(lambda,z,p1){
>> lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*pi0}
>>
>>  x <- seq(-20,20, by=0.1)
>> y <- numeric(length(x))
>> for (i in 1:length(x)){y[i] <- f1(x[i],p1=0.15,lambda=0.998)}
>> plot(y ~ x, ylim=c(-1,1))
>> abline(h=0)
>>
>>
>> a <- uniroot(f1, lower =-10, upper = 0,
>>tol = 1e-20,p1=0.15,lambda=0.998)$root
>>
>>
>>
>>
>>
>> 2015-04-15 22:57 GMT-04:00 Jeff Newmiller :
>>
>> > You really need to read the help page for uniroot. The sign needs to be
>> > different at the ends of the starting interval. This is a typical
>> > limitation of numerical root finders.
>> >
>> ---
>> > Jeff NewmillerThe .   .  Go
>> Live...
>> > DCN:Basics: ##.#.   ##.#.  Live
>> > Go...
>> >   Live:   OO#.. Dead: OO#..  Playing
>> > Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
>> > /Software/Embedded Controllers)   .OO#.   .OO#.
>> rocks...1k
>> >
>> ---
>> > Sent from my phone. Please excuse my brevity.
>> >
>> > On April 15, 2015 7:20:04 PM PDT, li li  wrote:
>> > >Hi all,
>> > >In the following code, I am trying to use uniroot function to solve for
>> > >the root (a and b in code below) for function f1.
>> > >I am not sure why uniroot function does not give the answer since when
>> > >we
>> > >look the graph, the function does cross 0 twice.
>> > >Any suggestion?
>> > >   Thanks.
>> > >   Hanna
>> > >
>> > >u1 <- -3
>> > >u2 <- 4
>> > >pi0 <- 0.8
>> > >
>> > >f1 <- function(lambda,z,p1){
>> > >lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*pi0}
>> > >
>> > >a <- uniroot(f1, lower =-10, upper = 0,
>> > >   tol = 1e-20,p1=0.15,lambda=0.998)$root
>> > >
>> > >b <- uniroot(f1, lower =0, upper = 10,
>> > >   tol = 1e-20,p1=0.15,lambda=0.998)$root
>> > >
>> > >x <- seq(-20,20, by=0.1)
>> > >y <- numeric(length(x))
>> > >for (i in 1:length(x)){y[i] <- f1(x[i],p1=0.15,lambda=0.998)}
>> > >plot(y ~ x, ylim=c(-1,1))
>> > >abline(h=0)
>> > >
>> > >   [[alternative HTML version deleted]]
>> > >
>> > >__
>> > >R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> > >https://stat.ethz.ch/mailman/listinfo/r-help
>> > >PLEASE do read the posting guide
>> > >http://www.R-project.org/posting-guide.html
>> <http://www.r-project.org/posting-guide.html>
>> > <http://www.r-project.org/posting-guide.html>
>> > >and provide commented, minimal, self-contained, reproducible code.
>> >
>> >
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> <http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Question with uniroot function

[li li]

Hi Jeff,
  Thanks for the reply. I am aware that the sign needs to be different at
the ends of the starting interval.

   Another question:

Is there a way to set the right end point ( (the "upper" argument in the
uniroot function below) as the point where the function takes on its
minimun, for example my function f1 below?

Thanks very much!



u1 <- -3
u2 <- 4
pi0 <- 0.8

f1 <- function(lambda,z,p1){
lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*pi0}

 x <- seq(-20,20, by=0.1)
y <- numeric(length(x))
for (i in 1:length(x)){y[i] <- f1(x[i],p1=0.15,lambda=0.998)}
plot(y ~ x, ylim=c(-1,1))
abline(h=0)


a <- uniroot(f1, lower =-10, upper = 0,
   tol = 1e-20,p1=0.15,lambda=0.998)$root





2015-04-15 22:57 GMT-04:00 Jeff Newmiller :

> You really need to read the help page for uniroot. The sign needs to be
> different at the ends of the starting interval. This is a typical
> limitation of numerical root finders.
> ---
> Jeff NewmillerThe .   .  Go Live...
> DCN:Basics: ##.#.   ##.#.  Live
> Go...
>   Live:   OO#.. Dead: OO#..  Playing
> Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
> /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
> ---
> Sent from my phone. Please excuse my brevity.
>
> On April 15, 2015 7:20:04 PM PDT, li li  wrote:
> >Hi all,
> >In the following code, I am trying to use uniroot function to solve for
> >the root (a and b in code below) for function f1.
> >I am not sure why uniroot function does not give the answer since when
> >we
> >look the graph, the function does cross 0 twice.
> >Any suggestion?
> >   Thanks.
> >   Hanna
> >
> >u1 <- -3
> >u2 <- 4
> >pi0 <- 0.8
> >
> >f1 <- function(lambda,z,p1){
> >lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*pi0}
> >
> >a <- uniroot(f1, lower =-10, upper = 0,
> >   tol = 1e-20,p1=0.15,lambda=0.998)$root
> >
> >b <- uniroot(f1, lower =0, upper = 10,
> >   tol = 1e-20,p1=0.15,lambda=0.998)$root
> >
> >x <- seq(-20,20, by=0.1)
> >y <- numeric(length(x))
> >for (i in 1:length(x)){y[i] <- f1(x[i],p1=0.15,lambda=0.998)}
> >plot(y ~ x, ylim=c(-1,1))
> >abline(h=0)
> >
> >   [[alternative HTML version deleted]]
> >
> >__
> >R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >https://stat.ethz.ch/mailman/listinfo/r-help
> >PLEASE do read the posting guide
> >http://www.R-project.org/posting-guide.html
> <http://www.r-project.org/posting-guide.html>
> >and provide commented, minimal, self-contained, reproducible code.
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Question with uniroot function

[li li]

Hi all,
   In the following code, I am trying to use uniroot function to solve for
the root (a and b in code below) for function f1.
I am not sure why uniroot function does not give the answer since when we
look the graph, the function does cross 0 twice.
Any suggestion?
   Thanks.
   Hanna

u1 <- -3
u2 <- 4
pi0 <- 0.8

f1 <- function(lambda,z,p1){
lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*pi0}

a <- uniroot(f1, lower =-10, upper = 0,
   tol = 1e-20,p1=0.15,lambda=0.998)$root

b <- uniroot(f1, lower =0, upper = 10,
   tol = 1e-20,p1=0.15,lambda=0.998)$root

x <- seq(-20,20, by=0.1)
y <- numeric(length(x))
for (i in 1:length(x)){y[i] <- f1(x[i],p1=0.15,lambda=0.998)}
plot(y ~ x, ylim=c(-1,1))
abline(h=0)

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Wrong results from anova

[li li]

0.003725088 is the mean square error.





2015-04-14 13:44 GMT-04:00 Michael Dewey :

> See in-line
>
>
> On 14/04/2015 18:23, li li wrote:
>
>> Hi all,
>>I have following data. When I perform an anova, the residual sum of
>> square returns to be zero.
>> But this is wrong, since we can hand calculate the RSS to be around
>> 0.0308.
>>Did anyone come across the same problem before? Any suggestions?
>>Thanks.
>>  Hanna
>>
>> ydata
>>>
>>  Y sample.Y
>> 1   0.477  0.5
>> 2   0.477  0.5
>> 3   0.478  0.5
>> 4  27.320   27
>> 5  27.420   27
>> 6  27.300   27
>> 7  29.440   29
>> 8  29.620   29
>> 9  29.610   29
>> 10 35.840   35
>> 11 35.900   35
>> 12 35.850   35
>>
>>> fit.Y <- aov(Y~sample.Y, data=ydata)
>>> summary(fit.Y)
>>>
>>  Df Sum Sq Mean Sq F value Pr(>F)
>> sample.Y 3   2203   734.2  190706 <2e-16 ***
>> Residuals8  0 0.0
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>>
>  734.2 / 197096
> [1] 0.003725088
>
> Did you mean 0.0038 above?
>
> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/
>> posting-guide.html <http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
> --
> Michael
> http://www.dewey.myzen.co.uk/home.html
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Wrong results from anova

[li li]

Hi all,
  I have following data. When I perform an anova, the residual sum of
square returns to be zero.
But this is wrong, since we can hand calculate the RSS to be around 0.0308.
  Did anyone come across the same problem before? Any suggestions?
  Thanks.
Hanna

> ydata
Y sample.Y
1   0.477  0.5
2   0.477  0.5
3   0.478  0.5
4  27.320   27
5  27.420   27
6  27.300   27
7  29.440   29
8  29.620   29
9  29.610   29
10 35.840   35
11 35.900   35
12 35.850   35
> fit.Y <- aov(Y~sample.Y, data=ydata)
> summary(fit.Y)
Df Sum Sq Mean Sq F value Pr(>F)
sample.Y 3   2203   734.2  190706 <2e-16 ***
Residuals8  0 0.0
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] xlims of box() function

[li li]

Hi all,
  It looks like the limits for x axis for box function in R is from -1 to
1. How can I change this limit?
  Thanks for your help.
 Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] title of r plots

[li li]

Thanks very much.
Also How do add an empty space when using expression()?
When I do the following, it returns an error message.


plot(1,1,main=expression(-70*degree*C%+-%10*degree*C/Ambient Condition))



   Hanna


2015-02-27 23:03 GMT-05:00 William Dunlap :

>  plot(1,1,main=expression(-70*degree*C%+-%10*degree*C/Ambient))
>
>  Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
>  On Fri, Feb 27, 2015 at 7:27 PM, li li  wrote:
>
>>   Hi all,
>>   I would like to add "-70°C ± 10°C/Ambient"  as the title of my plot.
>> Could anyone give some help on this?
>>   Thanks.
>>Hanna
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> <http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] title of r plots

[li li]

 Hi all,
  I would like to add "-70°C ± 10°C/Ambient"  as the title of my plot.
Could anyone give some help on this?
  Thanks.
   Hanna

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Finding MLE

[li li]

I think I made an error in my funciton. Now it works.


library(stats4)
n <- 8
ll<- function(a,b,x){
-sum(log(gamma((n-1)/2+a-1)/(gamma((n-1)/2)*gamma(a))*1/(2*b^a)*(x/2)^((n-1)/2-1)*(1/b+x/2)^(-((n-1)/2+a-1}
fit <- mle(ll, start=list(a=3, b=1), fixed=list(x=c(2,3)))





2014-11-02 22:04 GMT-05:00 li li :

> Thanks Bert for the reply. I still get a message when adding the start
> argument.
>
>
>
> > n <- 8
> > x0 <- c(2,3)
> >
> > ll<- function(a,b,x=x0,size=n){
> +
> -sum(log(gamma((n-1)/2+a-1)/(gamma((n-1)/2)*gamma(a))*1/(2*b^a)*(x/2)^((n-1)/2-1)*(1/b+x/2)^(-((n-1)/2+a-1}
> >
> > fit <- mle(ll, start=list(a=3, b=1), nobs = length(x0))
> Error in validObject(.Object) :
>   invalid class “mle” object: invalid object for slot "fullcoef" in class
> "mle": got class "list", should be or extend class "numeric"
> >
> >
>
>
>
>
>
>
> 2014-11-02 14:57 GMT-05:00 Bert Gunter :
>
> You do not appear to provide initial values for a and b , i.e. the
>> "start" argument for mle.
>>
>> Cheers,
>> Bert
>>
>> Bert Gunter
>> Genentech Nonclinical Biostatistics
>> (650) 467-7374
>>
>> "Data is not information. Information is not knowledge. And knowledge
>> is certainly not wisdom."
>> Clifford Stoll
>>
>>
>>
>>
>> On Sun, Nov 2, 2014 at 10:36 AM, li li  wrote:
>> > Hi all,
>> >   I am trying to use the mle function in R to find the maximum
>> likelihood
>> > estimator. The ll function below is the negative of the log likelihood.
>> > Suppose x0 is the observed values, I want to find the maximum likelihood
>> > for a and b. After running the code below, I get the error message
>> "Error
>> > in eval(expr, envir, enclos) : argument is missing, with no default".
>> >   Could anyone familiar with this function give some suggetion? Thanks
>> very
>> > much!
>> >  Hanna
>> >
>> >> n <- 8
>> >> x0 <- c(2,3)
>> >>
>> >> ll<- function(a,b,x=x0,size=n){
>> > +
>> >
>> -sum(log(gamma((n-1)/2+a-1)/(gamma((n-1)/2)*gamma(a))*1/(2*b^a)*(x/2)^((n-1)/2-1)*(1/b+x/2)^(-((n-1)/2+a-1}
>> >>
>> >> fit <- mle(ll, nobs = length(x0))
>> >
>> > [[alternative HTML version deleted]]
>> >
>> > __
>> > R-help@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
>> <http://www.r-project.org/posting-guide.html>
>> > and provide commented, minimal, self-contained, reproducible code.
>>
>
>

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Finding MLE

[li li]

Thanks Bert for the reply. I still get a message when adding the start
argument.



> n <- 8
> x0 <- c(2,3)
>
> ll<- function(a,b,x=x0,size=n){
+
-sum(log(gamma((n-1)/2+a-1)/(gamma((n-1)/2)*gamma(a))*1/(2*b^a)*(x/2)^((n-1)/2-1)*(1/b+x/2)^(-((n-1)/2+a-1}
>
> fit <- mle(ll, start=list(a=3, b=1), nobs = length(x0))
Error in validObject(.Object) :
  invalid class “mle” object: invalid object for slot "fullcoef" in class
"mle": got class "list", should be or extend class "numeric"
>
>






2014-11-02 14:57 GMT-05:00 Bert Gunter :

> You do not appear to provide initial values for a and b , i.e. the
> "start" argument for mle.
>
> Cheers,
> Bert
>
> Bert Gunter
> Genentech Nonclinical Biostatistics
> (650) 467-7374
>
> "Data is not information. Information is not knowledge. And knowledge
> is certainly not wisdom."
> Clifford Stoll
>
>
>
>
> On Sun, Nov 2, 2014 at 10:36 AM, li li  wrote:
> > Hi all,
> >   I am trying to use the mle function in R to find the maximum likelihood
> > estimator. The ll function below is the negative of the log likelihood.
> > Suppose x0 is the observed values, I want to find the maximum likelihood
> > for a and b. After running the code below, I get the error message "Error
> > in eval(expr, envir, enclos) : argument is missing, with no default".
> >   Could anyone familiar with this function give some suggetion? Thanks
> very
> > much!
> >  Hanna
> >
> >> n <- 8
> >> x0 <- c(2,3)
> >>
> >> ll<- function(a,b,x=x0,size=n){
> > +
> >
> -sum(log(gamma((n-1)/2+a-1)/(gamma((n-1)/2)*gamma(a))*1/(2*b^a)*(x/2)^((n-1)/2-1)*(1/b+x/2)^(-((n-1)/2+a-1}
> >>
> >> fit <- mle(ll, nobs = length(x0))
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@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
> <http://www.r-project.org/posting-guide.html>
> > and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@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.

[R] Finding MLE

[li li]

Hi all,
  I am trying to use the mle function in R to find the maximum likelihood
estimator. The ll function below is the negative of the log likelihood.
Suppose x0 is the observed values, I want to find the maximum likelihood
for a and b. After running the code below, I get the error message "Error
in eval(expr, envir, enclos) : argument is missing, with no default".
  Could anyone familiar with this function give some suggetion? Thanks very
much!
 Hanna

> n <- 8
> x0 <- c(2,3)
>
> ll<- function(a,b,x=x0,size=n){
+
-sum(log(gamma((n-1)/2+a-1)/(gamma((n-1)/2)*gamma(a))*1/(2*b^a)*(x/2)^((n-1)/2-1)*(1/b+x/2)^(-((n-1)/2+a-1}
>
> fit <- mle(ll, nobs = length(x0))

[[alternative HTML version deleted]]

__
R-help@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.

[R] Making several plots using a loop function

[li li]

Hi all,
   When making a bunch of plots using a loop function, how to add title to
reflect different plots.
Specifically, for the code below, I generated 9 plots.  I would like to add
a title to each plot.
For example, the titles will be respectively, plot1, plot 2, … plot 9.
Thank you very much!
Hanna



par(mfrow=c(3,3), pty="s", pch=16, col="blue")
for ( i in 1:9) {
lm1 <-  lm(log(run_res[2:11, i]) ~ log(level))
plot(log(level), log(run_res[2:11, i]), xlab="dllution levels",
ylab="%PS-80")
abline(lm1)}

[[alternative HTML version deleted]]

__
R-help@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.

[R] rsm package

[li li]

Hi all,
  I fit a complete second order model for my response with three
predictors. However, the lack of fit of the model is
still significant. I wish to add the third order terms also. Is there a way
to do that using rsm function?
  Thanks.
 Hanna

[[alternative HTML version deleted]]

__
R-help@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.

[R] Proportional Odds Model

[li li]

Hi all,
  I am trying to perform an ordinal regression using proportional odds
model.
First I would like to test whether the assumption of proportional odds is
supported by the data.
I got the following error message when I try to fit one of the relevant
models. Can anyone give
some suggestions on how to fix this. Thank you.
   Hanna


 fit2 <- vglm(PLANHORIZON ~ factor(TotalCompanion) +
factor(NumTrips)+factor(Numnight)+factor(FMWEBSITE)+factor(OTHERWEB)+factor(BROCHU)
+ +factor(TRAAGENT) + factor(WOM) + factor(AAA) + factor(GUIDEBOOK) +
factor(ADVERTISE) + factor(NEWSPAPER) + factor(INFOOTHER) +
+ ExpTotal + factor(AGE) + factor(CHILD) + factor(INCOME), +
family=cumulative(parallel=F),  data=two)

Error in if (any(smallmu <- (mu * (1 - mu) < double.eps))) { :
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
residuals,  :
  fitted values close to 0 or 1
2: In Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
residuals,  :
  fitted values close to 0 or 1

[[alternative HTML version deleted]]

__
R-help@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.

[R] Change color of the boxplot outliers

[li li]

Hi all,
  Is there a way to change the color of the boxplot plots outliers?
  Thanks.
Hanna

[[alternative HTML version deleted]]

__
R-help@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.

[R] Bayesian Variance Components Estimation

[li li]

Hi all,
   Does anyone know whether there is an R function available to find the
Bayesian variance components estimators
for random effects or mixed models?
  Thank you very much.
Hanna

[[alternative HTML version deleted]]

__
R-help@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.

[R] How can C++ read the R object written into socket with saveRDS or save

[Rong lI Li]

Hi, all,

Recently, I met one issue when using socket between R & C++ to transmit R
object. Would you pls help give me some suggestions? Many thanks!

[Background]:
I create a socket connection between R & C++ binary first, and then, want
to use saveRDS() or save() in R to save the object into connection
directly. So that the C++ binary can read the object, and send it to
another remote R.

[What I did so far]:
1. I used socketConnection in R and listen/accept in C++, to establish one
blocking socket.
2. I used saveRDS to save the R object into socket directly
3. I want to use "recv()" in C++ to receive the R object.

[Issues I met]:
I found actually, the saveRDS writes the R object with XDR format. I could
not know how many bytes are sent into socket, when calling saveRDS to save
R object. So in the C++ binary, I could not know exactly how many bytes I
should receive from the socket. It is not safe for me, to always use a
pre-defined buffer size to read from the socket.

Any suggestions for this? Are there safe way for me to read the R object
from the socket?
I do not do any conversion with the received data, and only need to
transfer them into a remote R to do the function execution.

=

Rong "Jessica", Li (ÀîÈÙ)
Platform Symphony TET, CSTL, IBM Systems &Technology Group, Development
Tel:86-10-82451010  Email:rong...@cn.ibm.com
[[alternative HTML version deleted]]

__
R-help@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.