I was told that the response could be a matrix of numbers positives and numbers negatives (at each site) as well as a string of positives and negatives (i.e. 1 row for each individual). This also seemed strange to me and like you, to date, I have only used 1 row per observation. But if I do model the data using 1 obs per row with glmmPQL, i.e.

y<-c(1,rep(0,5),1,1,rep(0,4),rep(1,3),rep(0,3),rep(1,4),rep(0,2),rep(1,5),0)Facility<-sort(rep(1:5,6))Cov<-Facility*10LONG<-c(rep(0,6),rep(1,6),rep(1,6),rep(0,6),rep(1,6))LAT<-c(rep(0,6),rep(0,6),rep(1,6),rep(2,6),rep(2,6))D<-data.frame(y,Facility,Cov,LONG,LAT)
model= glmmPQL(y  ~ Cov,  random = ~ 1 | Facility, correlation=corExp(form=~LONG+LAT|Facility),            data=D,family=binomial(logit))
I get the error 
Error in getCovariate.corSpatial(object, data = data) :   cannot have zero distances in "corSpatial"
It works fine without the spatial correlation... 
Thanks for your help anyway!
HughDate: Thu, 16 May 2013 16:44:40 -0400
Subject: Re: [R-sig-Geo] Spatially autocorrelated random effects in logistic regression.
From: sjmyers3142@gmail.com
To: hughsturrock@hotmail.com

I am not familiar with what you mean by grouping using a matrix of positive and negatives.  I assume this is coding your data in some means other than long form where each observation has its own row.  Have you tried the more traditional data method?  When I run your Y<-cbind(1:5,5:1) in R I get this returned:

> Y<-cbind(1:5,5:1)> Y     [,1] [,2][1,]    1    5[2,]    2    4[3,]    3    3
[4,]    4    2[5,]    5    1> 
Which does not seem like what you have described or anything that I can see entering a logistic regression.



I typically bind all columns into a single data.frame using cbind with one row per observation (with grouping factors repeated among rows as necessary), use 0 and 1 as my response in the binomial case, and then specify that data.frame using data= argument.  This is the method that is shared across all R functions that I have used.  I'd hate for you to be thrown an error over something trivial like that.  

The other thing to look at is that in your toy data set, it seems that you only have one observation per grouping factor that you submit to random= argument.  That might also be throwing an error.  I would try a different example with a few more data points and also not make my output a function of the covariate with no error.  A single observation per group might cause issues because it would result in perfectly fitting models (which binomial models can't handle oddly), also if your Y response values are perfectly predictable (differentiable by rank) by X that will also cause this issue.  It is a weird thing that is found in logistic regression as well.  Data sets where a linear function of the covariates can perfectly divide the Ys into 0 and 1 groups will throw errors.  You need some slop, so to speak.  I can see your Ys (unless the print out above is what it actually is) but I get a feeling that they might fall into this trap.  

As you can tell, I have gotten glmmPQL to work for my data in simpler cases and have used a good bit of other spatial functions, but I am no expert in glmmPQL.  I thought I could help more easily.  Perhaps what I would look at next in terms of trouble shooting (using a slightly more complex toy data set to look at these things) might help.

Seth

On Thu, May 16, 2013 at 4:24 PM, Hugh Sturrock <hughsturrock@hotmail.com> wrote:




Hi Seth,
Thanks for your input on this. Unfortunately it doesn't seem like it will accept 2 grouping factors. In fact, it doesn't let you specify different groups if your response is already grouped, i.e. as a matrix of positives and negatives. For example

Y<-cbind(1:5,5:1)Long<-c(0,1,1,0,1)Lat<-c(0,0,1,2,2)Covar<-seq(10,50,10)Group=factor(rep("a",nrow(Y)))g<-factor(1:5)

All of these produce errors:
# run a model just with random=argument for householdmodel1 = glmmPQL(Y  ~ Covar,  random = ~ 1 | g, correlation=corExp(form=~Long+Lat), family=binomial(logit))

model2 = glmmPQL(Y  ~ Covar,  random = ~ 1 | g, correlation=corExp(form=~Long+Lat|g), family=binomial(logit))

model3 = glmmPQL(Y  ~ Covar,  random = ~ 1 | g, correlation=corExp(form=~Long+Lat|Group), family=binomial(logit)) ## This is the one I really want..!
# The only model that runs ismodel4 = glmmPQL(Y  ~ Covar,  random = ~ 1 | Group, correlation=corExp(form=~Long+Lat), family=binomial(logit)) 

# This seems to be the same as model4, i.e. model4 is borrowing the group level specified as random in the correlation model5 = glmmPQL(Y  ~ Covar,  random = ~ 1 | Group, correlation=corExp(form=~Long+Lat | Group), family=binomial(logit)) 

# It doesn't seem to like the corMatrix command on glmmPQL modelscorMatrix(model5) Error in UseMethod("corMatrix") :   no applicable method for 'corMatrix' applied to an object of class "c('glmmPQL', 'lme')"


I've done a fair amount of reading on the use of glmmPQL but nowhere can I find anything on this specific problem.. 
Date: Thu, 16 May 2013 15:27:56 -0400

Subject: Re: [R-sig-Geo] Spatially autocorrelated random effects in logistic regression.
From: sjmyers3142@gmail.com
To: hughsturrock@hotmail.com


Hi,
I forgot about the bit of "tricking" glmmPQL by adding a dummy variable that is the same for each observation.  Do you have SAS by any chance?  This may be a little easier in Proc Glimmix.


I believe you can have two grouping factors.  You can have one grouping factor where each household has its own grouping, so household one is "1" (as a factor) and household 2 is "2" and there would be multiple observations within each nesting at that grouping level.  This is what you would specify in your random = argument if you want to use a random intercept to account for correlations within a single household.  You can then specify another grouping factor where all observations are a "1" or some other factor.  Then you can use the | symbol within the correlation= argument after the east + north to model between household spatial autocorrelation.  It seems that glmmPQL might be "borrowing" the grouping factor for correlation= from the random=, which is odd if true (I just skimmed an attached chapter, and this seems to be the case, but perhaps you can force two different grouping factors one at the household and one where the entire study area is a group).


Try specifying random as you had it with house as the group and then correlation as 
correlation=corExp(form=~Long+Lat | Group)



I believe the partial sill and nugget effect won't do anything if the coordinates are exactly the same.  I am not 100% certain though.  It is down to how it was coded.  You can post this question to the appropriate mailing list and maybe get the package authors to respond.  I would try it though to see in the following way.


run a model just with random=argument for householdrun a model just with correlation=argument with all data in one big group (no nugget)run a model just with correlation=argument with all data in one big group (WITH nugget)


run a model just with random=argument for household AND with correlation=argument with all data in one big group (no nugget)
run a model just with random=argument for household AND with correlation=argument with all data in one big group (WITH nugget)



If you save each model and look at the correlation matrices for each and compare them, you can likely see what is going on.
Here's a book chapter that I read when I first ran glmmPQL.  The relevant part starts around page 301.


Seth



On Thu, May 16, 2013 at 2:06 PM, Hugh Sturrock <hughsturrock@hotmail.com> wrote:





Yes that's what I thought. I need to specify a different group for each household, however, as your pasted section states, if a grouping factor is specified, the correlation structure is assumed to apply to observations WITHIN the grouping level, not BETWEEN. I would like to do both, i.e. account for correlation within groups and spatial autocorrelation between groups. 


If I run the model with the response as a vector of 1's and 0's I have the problem of cooincident coordinates. If I run the model with the response as a matrix of positives and negatives, when I specify each row of the response as a separate group, it gives me the errors I posted before. The only way seems to be to force each row of the reponse (which is several observations of positives and negatives at each site) to be a single 'observation', create one grouping factor for all sites and then specify the correlation structure which will model the spatial autocorrelation between 'observations' (which in this case are my numbers positive and numbers negative at each site). Does that make sense? I'm just worried that this doesn't account for correlation within sites, but maybe it does, i.e. the correlation at a distance of 0 (i.e. the partial sill of the variogram)..??  


spBayes was suggested to me, however, I really wanted to find a frequentist way round this as my understanding of Bayesian is still quite limited..


Date: Thu, 16 May 2013 13:30:19 -0400


Subject: Re: [R-sig-Geo] Spatially autocorrelated random effects in logistic regression.
From: sjmyers3142@gmail.com
To: hughsturrock@hotmail.com



Hi,
I read the link.  If each household has its own factor indicting a group (your code seems to make them all the same), then the random= argument will take care of within household correlations by assigning a random intercept to each.  The addition of the intercept will allow each household to have higher or lower values, and they move as a block and so this creates the within house correlation (if I add 2 inches to all household heights across the street then that would create a household effect, for instance).  The correlation = argument as you have used it does not know anything about groupings, so it will treat all observations at the finest level (individual observations).  If two observations have the same spatial location (same household) then I believe that the correlation= argument will not do anything to those locations right on top of each other, it will treat them as the same point in space, which they are.  This prevents things like the same data point influencing itself, which it is of course correlated with.  It will, however, incorporate the correlations from other households, which seems what you would like.  glmmPQL makes multiple calls to lme in library(nlme).  If you look at the documentation for correlation structures (such as corExp) you will see this:



"m a one sided formula of the form ~ S1+...+Sp, or ~ S1+...+Sp | g, specifying spatial covariates S1 through Sp and, optionally, a grouping factor g. When agrouping factor is present in form, the correlation structure is assumed to apply


only to observations within the same grouping level; observations with different
grouping levels are assumed to be uncorrelated. Defaults to ~ 1, which corresponds to using the order of the observations in the data as a covariate, and no


groups"
Sorry for the pdf to gmail copy and paste debacle.
As you can see, if you do not impose a group using | then it does not seem to group them.  If it did, then it would have to "borrow" the grouping from the random= argument, which seems nonsensical.  Here's how you can test it.  Run the model and then use corMatrix (see nlme documentation, it is under corMatrix.corStruct).  You can then see the correlation matrices the model used.  You can then see if the model behaved as you'd like.



Seth






On Thu, May 16, 2013 at 12:55 PM, Hugh Sturrock <hughsturrock@hotmail.com> wrote:




Seth, thanks for your help. I was modeling this as follows
Y<-cbind(1:5,5:1) # matrix of numbers positive, numbers negative for each household


House<-1:5Long<-c(0,1,1,0,1)Lat<-c(0,0,1,2,2)Covar<-House*10Group=factor(rep("a",nrow(Y))) # I've read that you have to specify this as the model needs a 'random' input




model1 = glmmPQL(Y  ~ Covar,  random = ~ 1 | Group, correlation=corExp(form=~Long+Lat), family=binomial(logit))summary(model1)
This seems to run fine, however, I'm concerned that it isn't quite doing what I want. How is this model accounting for correlation between observations in the same household? Is it accounting for spatial autocorrelation between households? The response to this post seems to suggest not



http://grokbase.com/t/r/r-help/125x5gnv6m/r-glmmpql-spatial-autocorrelation



Initially I thought the 'random' input should be a factor identifying the house from which the observations arise, i.e. 1:5, but if I specify it like this
H<-factor(House)


model2 = glmmPQL(Y  ~ Covar,  random = ~ 1 | H, correlation=corExp(form=~Long+Lat), family=binomial(logit))
I get the following error:
Error in corFactor.corSpatial(object) : 


  NA/NaN/Inf in foreign function call (arg 1)In addition: Warning messages:1: In min(unlist(attr(object, "covariate"))) :  no non-missing arguments to min; returning Inf


2: In min(unlist(attr(object, "covariate"))) :  no non-missing arguments to min; returning Inf

I'm slightly lost as to whether glmmPQL can be used for what I want and if so whether the model 1 above is the right way to specify it...? Any thoughts?





Date: Thu, 16 May 2013 10:44:28 -0400
Subject: RE: [R-sig-Geo] Spatially autocorrelated random effects in logistic regression.

From: sjmyers3142@gmail.com
To: hughsturrock@hotmail.com


Message: 5Date: Thu, 16 May 2013 01:07:41 +0100



From: Hugh Sturrock <hughsturrock@hotmail.com>



To: "r-sig-geo@r-project.org" <r-sig-geo@r-project.org>



Subject: [R-sig-Geo] Spatially autocorrelated random effects in        logistic        regression.



Message-ID: <DUB116-W3414EE6B601E68C60BB108C8A30@phx.gbl>



Content-Type: text/plain



Hi all,I am trying to model binomial data (infection status) of people living in households. I have run a logistic regression with a household level random effect to account for correlation between individuals within households. When I look at the value of the random effect (intercept) at each household there appears to be residual spatial autocorrelation (using Moran's I and semivariograms) between households. I would like to repeat the analyses accounting for this spatial autocorrelation between households, i.e. between random effects. I'm not a mathematician  but having done some reading it seems as if this can't be done using the command glmmPQL in MASS, but perhaps can be done using the geoRglm package. BUT, having read through the geoRglm vignette, I'm still pretty clueless as to how to code such a model. If anyone could give me an example I would be MASSIVELY appreciative.



Thanks!Hugh



        [[alternative HTML version deleted]]




Dear Hugh,

Does glmmPQL fail because you have multiple instances with the same spatial coordinates?  glmmPQL can handle spatial autocorrelation, but I have not tried it with your sampling design.  Here is an example for glmmPQL in case it helps.




You must have library (nlme) and library (MASS) loaded in R before you do this.  Also, you must have 2 columns representing spatial coordinates.  I'll call them east and north.





Here is a small example with two variables.




spatialregress1<-glmmPQL(outcome ~ var1 + var2, data=yourdata, family=binomial (link="logit"), correlation=corExp (form= ~ east + north))




You can include a nugget effect by correlation=corExp (form= ~ east + north,nugget=TRUE)




I haven't run this sample code, so I hope there is no small typo.  But the general form should be correct.




There are multiple spatial correlation structures from which to choose.  I used corExp in my example because it has often been useful for me.  You already know how to examine residuals, it seems.  If it were me, I would look at the residual variogram from the model without spatial correlation and use library (gstat) to attempt to fit several different spatial structures (exponential, spherical, etc) to the residuals using fit.variogram.  I visually examine them and look for best fit at relatively short lags where autocorrelation is greatest.  This will also allow you to extract start values for the correlation parameters (range, sill, and possibly nugget) and supply those to glmmPQL.  I do not know off the top of my head the order for passing those to corExp, so you might have to do some digging there.




I am not sure with glmmPQL how to extract residuals that are conditioned on the correlation structure, but this might do it.




res_spat<-residuals(spatialregress1, type ='normalized')




Hope that gets you started.

Seth Myers








 		 	   		  

 		 	   		  

 		 	   		  

 		 	   		  
	[[alternative HTML version deleted]]

