R help - Oct 2008 - call works with gee and yags, but not geepack

I have included data at the bottom of this email. It can be read in by
highlighting the data and then using this command: dat <-
read.table("clipboard", header = TRUE,sep="\t")
I can obtain solutions with both of these:

library(gee)

fit.gee<-gee(score ~ chem + time, id=id,
family=gaussian,corstr="exchangeable",data=dat)

and

library(yags)
fit.yags <- yags(score ~ chem + time, id=id,
family=gaussian,corstr="exchangeable",data=dat,alphainit=0.05)

However, I am making a mistake with:

library(geepack)
fit.geese <- geese(score ~ chem + time, id=id,
family=gaussian,corstr="exch",data=dat)

I obtain the following error:

Error in geese.fit(x, y, id, offset, soffset, w, waves, zsca, zcor, corp,  :
 nrow(zsca) and length(y) not match

Could someone tell me what I have done incorrectly. Thanks for your
time, Juliet.

Data Below:

id treat time1 time2 time3 time4 chem1 chem2 chem3 chem4 time score chem
1 1 20 18 15 15 1000 1100 1200 1300 0 20 1000
1 1 20 18 15 15 1000 1100 1200 1300 2 18 1100
1 1 20 18 15 15 1000 1100 1200 1300 3 15 1200
1 1 20 18 15 15 1000 1100 1200 1300 6 15 1300
2 1 22 24 18 22 1000 1000 1055 950 0 22 1000
2 1 22 24 18 22 1000 1000 1055 950 2 24 1000
2 1 22 24 18 22 1000 1000 1055 950 3 18 1055
2 1 22 24 18 22 1000 1000 1055 950 6 22 950
3 1 14 10 24 10 1000 1999 800 1700 0 14 1000
3 1 14 10 24 10 1000 1999 800 1700 2 10 1999
3 1 14 10 24 10 1000 1999 800 1700 3 24 800
3 1 14 10 24 10 1000 1999 800 1700 6 10 1700
4 1 38 34 32 24 1000 1050 1400 0 38 1000
4 1 38 34 32 24 1000 1050 1400 2 34 1050
4 1 38 34 32 24 1000 1050 1400 3 32
4 1 38 34 32 24 1000 1050 1400 6 24 1400
5 1 25 29 25 29 1000 1020 1040 1045 0 25 1000
5 1 25 29 25 29 1000 1020 1040 1045 2 29 1020
5 1 25 29 25 29 1000 1020 1040 1045 3 25 1040
5 1 25 29 25 29 1000 1020 1040 1045 6 29 1045
6 1 30 28 26 14 1000 1100 1180 1500 0 30 1000
6 1 30 28 26 14 1000 1100 1180 1500 2 28 1100
6 1 30 28 26 14 1000 1100 1180 1500 3 26 1180
6 1 30 28 26 14 1000 1100 1180 1500 6 14 1500
7 0 20 15 21 20 1000 1200 1200 1000 0 20 1000
7 0 20 15 21 20 1000 1200 1200 1000 2 15 1200
7 0 20 15 21 20 1000 1200 1200 1000 3 21 1200
7 0 20 15 21 20 1000 1200 1200 1000 6 20 1000
8 0 21 27 19 27 1000 900 1075 900 0 21 1000
8 0 21 27 19 27 1000 900 1075 900 2 27 900
8 0 21 27 19 27 1000 900 1075 900 3 19 1075
8 0 21 27 19 27 1000 900 1075 900 6 27 900
9 0 15 22 22 20 1000 1000 700 0 15 1000
9 0 15 22 22 20 1000 1000 700 2 22
9 0 15 22 22 20 1000 1000 700 3 22 1000
9 0 15 22 22 20 1000 1000 700 6 20 700
10 0 39 39 34 1000 950 1033 1025 0 39 1000
10 0 39 39 34 1000 950 1033 1025 2 950
10 0 39 39 34 1000 950 1033 1025 3 39 1033
10 0 39 39 34 1000 950 1033 1025 6 34 1025
11 0 27 27 31 22 1000 950 910 1050 0 27 1000
11 0 27 27 31 22 1000 950 910 1050 2 27 950
11 0 27 27 31 22 1000 950 910 1050 3 31 910
11 0 27 27 31 22 1000 950 910 1050 6 22 1050
12 0 28 24 33 1000 1015 985 0 28 1000
12 0 28 24 33 1000 1015 985 2 24 1015
12 0 28 24 33 1000 1015 985 3 33 985
12 0 28 24 33 1000 1015 985 6

Sorry. I did not output the NAs correctly.

dat <- read.table("clipboard", header = TRUE)

id treat time1 time2 time3 time4 chem1 chem2 chem3 chem4 time score chem
1 1 20 18 15 15 1000 1100 1200 1300 0 20 1000
1 1 20 18 15 15 1000 1100 1200 1300 2 18 1100
1 1 20 18 15 15 1000 1100 1200 1300 3 15 1200
1 1 20 18 15 15 1000 1100 1200 1300 6 15 1300
2 1 22 24 18 22 1000 1000 1055 950 0 22 1000
2 1 22 24 18 22 1000 1000 1055 950 2 24 1000
2 1 22 24 18 22 1000 1000 1055 950 3 18 1055
2 1 22 24 18 22 1000 1000 1055 950 6 22 950
3 1 14 10 24 10 1000 1999 800 1700 0 14 1000
3 1 14 10 24 10 1000 1999 800 1700 2 10 1999
3 1 14 10 24 10 1000 1999 800 1700 3 24 800
3 1 14 10 24 10 1000 1999 800 1700 6 10 1700
4 1 38 34 32 24 1000 1050 NA 1400 0 38 1000
4 1 38 34 32 24 1000 1050 NA 1400 2 34 1050
4 1 38 34 32 24 1000 1050 NA 1400 3 32 NA
4 1 38 34 32 24 1000 1050 NA 1400 6 24 1400
5 1 25 29 25 29 1000 1020 1040 1045 0 25 1000
5 1 25 29 25 29 1000 1020 1040 1045 2 29 1020
5 1 25 29 25 29 1000 1020 1040 1045 3 25 1040
5 1 25 29 25 29 1000 1020 1040 1045 6 29 1045
6 1 30 28 26 14 1000 1100 1180 1500 0 30 1000
6 1 30 28 26 14 1000 1100 1180 1500 2 28 1100
6 1 30 28 26 14 1000 1100 1180 1500 3 26 1180
6 1 30 28 26 14 1000 1100 1180 1500 6 14 1500
7 0 20 15 21 20 1000 1200 1200 1000 0 20 1000
7 0 20 15 21 20 1000 1200 1200 1000 2 15 1200
7 0 20 15 21 20 1000 1200 1200 1000 3 21 1200
7 0 20 15 21 20 1000 1200 1200 1000 6 20 1000
8 0 21 27 19 27 1000 900 1075 900 0 21 1000
8 0 21 27 19 27 1000 900 1075 900 2 27 900
8 0 21 27 19 27 1000 900 1075 900 3 19 1075
8 0 21 27 19 27 1000 900 1075 900 6 27 900
9 0 15 22 22 20 1000 NA 1000 700 0 15 1000
9 0 15 22 22 20 1000 NA 1000 700 2 22 NA
9 0 15 22 22 20 1000 NA 1000 700 3 22 1000
9 0 15 22 22 20 1000 NA 1000 700 6 20 700
10 0 39 NA 39 34 1000 950 1033 1025 0 39 1000
10 0 39 NA 39 34 1000 950 1033 1025 2 NA 950
10 0 39 NA 39 34 1000 950 1033 1025 3 39 1033
10 0 39 NA 39 34 1000 950 1033 1025 6 34 1025
11 0 27 27 31 22 1000 950 910 1050 0 27 1000
11 0 27 27 31 22 1000 950 910 1050 2 27 950
11 0 27 27 31 22 1000 950 910 1050 3 31 910
11 0 27 27 31 22 1000 950 910 1050 6 22 1050
12 0 28 24 33 NA 1000 1015 985 NA 0 28 1000
12 0 28 24 33 NA 1000 1015 985 NA 2 24 1015
12 0 28 24 33 NA 1000 1015 985 NA 3 33 985
12 0 28 24 33 NA 1000 1015 985 NA 6 NA NA

On 30/10/2008, at 9:08 AM, Juliet Hannah wrote:
> I have included data at the bottom of this email. It can be read in by
> highlighting the data and then using this command: dat <-
> read.table("clipboard", header = TRUE,sep="\t")
> I can obtain solutions with both of these:
>
> library(gee)
>
> fit.gee<-gee(score ~ chem + time, id=id,
> family=gaussian,corstr="exchangeable",data=dat)
>
> and
>
> library(yags)
> fit.yags <- yags(score ~ chem + time, id=id,
> family=gaussian,corstr="exchangeable",data=dat,alphainit=0.05)
>
> However, I am making a mistake with:
>
> library(geepack)
> fit.geese <- geese(score ~ chem + time, id=id,
> family=gaussian,corstr="exch",data=dat)
>
> I obtain the following error:
>
> Error in geese.fit(x, y, id, offset, soffset, w, waves, zsca, zcor,  
> corp,  :
>  nrow(zsca) and length(y) not match
>
> Could someone tell me what I have done incorrectly. Thanks for your
> time, Juliet.
I'm pretty sure this is a bug in geese(), which should be reported to  
the
maintainer of geepack.  The problem is with the treatment of missing  
values.

If looks at
dim(na.omit(dat[,c("id","score","chem","time")]))
one
gets 44.
In geese.fit() zsca is set equal to matrix(1,N,1) where N is set  
equal to
length(id).  But id has length 46 whereas the response y has been  
trimmed
down to length 44 by eliminating any rows of the data where any of  
the variables
involved are missing.  Hence a problem.

The solution of the problem requires some code re-writing by the  
maintainer of geepack.

	cheers,

		Rolf Turner

######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}