[R] lme random effects in additive models with interaction

elifnurdogruoz elifnurdogruoz at mynet.com
Thu Jun 21 21:35:27 CEST 2012


Thanks for your answer, I would like to make clear my question: 

My data is like following and there is a response variable y:

Time	      Size	     Charge	Density	Replication
3             small             +               low		1            
.		.		       .		.		.
.		.           	       .		.		.	
9             small             +               low	    	1      
3             big                 +               low	     	.       
.		.		       .		.		.
.               .           	       .		.       	.           
9             big                 +              	low     	1       
3             small             -               low     	1      
.		.		       .		.		.
.               .           	       .		.       	.                  
9             small             -               low     	1       
3             big                 -               low     	1       
.		.		       .		.		. 
.               .           	       .		.       	.                  
9             big             	-               low 	1
3             small             +               high    	1       
.		.		       .		.		.
.               .           	       .		.		.	
9             small             +               high    	1         
3             big                 +               high    	1         
.		.		       .		.		. 
.               .           	       .		.       	.             
9             big                 +              	high    	1          
3             small             -               high    	1         
.		.		       .		.		. 
.               .           	       .		.       	1                   
9             small             -               high    	1        
3             big                 -               high    	1        
.		.		       .		.		. 
.               .           	       .		.       	.                 
9             big                 -               high		1
3             small             +               low		2            
.		.		       .		.		.
.		.           	       .		.		.	
9             small             +               low	    	2      
3             big                 +               low	     	2     
.		.		       .		.		.
.               .           	       .		.       	.           
9             big                 +              	low     	2       
3             small             -               low     	2      
.		.		       .		.		.
.               .           	       .		.       	.                  
9             small             -               low     	2       
3             big                 -               low     	2       
.		.		       .		.		. 
.               .           	       .		.       	.                  
9             big             	-               low 	2
3             small             +               high    	2       
.		.		       .		.		.
.               .           	       .		.		.	
9             small             +               high    	2         
3             big                 +               high    	2         
.		.		       .		.		. 
.               .           	       .		.       	.             
9             big                 +              	high    	2          
3             small             -               high    	2         
.		.		       .		.		. 
.               .           	       .		.       	.                   
9             small             -               high    	2        
3             big                 -               high    	2       
.		.		       .		.		. 
.               .           	       .		.       	.                 
9             big                 -               high		2


My code with comments:

##this function selects the knots
default.knots <- function(x,num.knots)
{
if (missing(num.knots))
num.knots <- max(5,min(floor(length(unique(x))/4),35))
return(quantile(unique(x),seq(0,1,length=
(num.knots+2))[-c(1,(num.knots+2))]))
}

knots <- default.knots(Time)

z <- outer(Time, knots, "-") 
z <- z * (z > 0)
z<-z^2

i.size50 <- I(Size==small)
i.chargepos <- I(Charge=="+")
i.densitylow <- I(Density==low)

##Create X and Z matrices,  I put interactions because I want intercept to
be zero at time 0.
X <- cbind( I(Time^2),Time*i.size50,Time*i.chargepos,Time*i.densitylow)
Z <- cbind( z, z*i.size50, z*i.chargepos,z*i.densitylow)


K <- length(knots)

## form blocked diagonal matrix Z to specify which columns of Z are used for
each group

block.ind <- list(1:K, (K+1):(2*K),(2*K+1):(3*K),(3*K+1):(4*K))
Z.block <- list()
for (i in 1:length(block.ind))
Z.block[[i]] <-
as.formula(paste("~Z[,c(",paste(block.ind[[i]],collapse=","),")]-1"))

##create dummy grouping variable since groupedData object is required for
lme
group <- rep(1, length(Time))
model.data <- groupedData(y~X|group, data=data.frame(X, y))

fit <- lme(y~-1+X, data=model.data, random=pdBlocked(list(
pdBlocked(Z.block,pdClass="pdIdent"), pdIdent(~-1+ Replication) ))
,control=list(maxIter=1000, msMaxIter=1000, niterEM=1000))

The experiment is repeated twice (Replication 1 and 2) , hence I think that
Replication should be random effect. As you said,  my replications are
randomly chosen from a population and I should make inference about the
population.  I don't have a chance to take more replications. Then, I am
planning to generate new data sets from the fitted model.

--
View this message in context: http://r.789695.n4.nabble.com/lme-random-effects-in-additive-models-with-interaction-tp4634067p4634143.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list