Thursday, December 07, 2006

multinom

Hi Tessa,
To compare the fit of two models, you can first fit each one separately, and compare them using the anova() function. It will do a likelihood ratio test.
I think you have to construct your response differently, and not use "replica" to group observations.
The response variable has to be a matrix, with a row for each replica, and the counts on each patch in the columns. The explanatory variables will have as many entries as there are lines in the matrix. I now see that your covariates are patch-specific, which is more involved than what multinom() allows. Multinom() only allows an offset with that structure.
Maybe you will have to fit your data using a log-linear poisson model then, or function vglm() from the "VGAM" library can also handle multinomial data. Agresti, Categorical Data Analysis has a lot of info on the analysis of multinomial data.
Scripts for that book are available at http://www-stat.stanford.edu/~owen/courses/306a/Splusdiscrete2.pdf.
Cheers, Tom

Tuesday, December 05, 2006

Statistics Forum Functional Ecology Research School

Dear Tom,

I have tried to use the function multinom on my data. My experiment consists of multiple replicas of batches of 15 mites that distribute themselves over 5 patches. Thus, the response variable is multinomial distributed. I included 3 more explanatory variables as well (idad, nceggs and ideggs, assuming for now that they are independent).

My data looks like this (but then I have some more replica’s)

replica patch ncad idad ideggs nceggs
aaaa a 2 0 0 2
aaaa b 10 5 1 17
aaaa c 0 1 0 2
aaaa d 0 2 12 1
aaaa e 0 7 4 2
bbbb a 1 3 5 1
bbbb b 9 3 1 8
bbbb c 0 0 1 4
bbbb d 2 0 11 6
bbbb e 2 9 4 2

My question are

1. whether I have used the function multinom correctly? I have tried to find documentary which I could understand, but I failed. For example Venables and Ripley spent less then a page on it.

2 how should I define “replica”. Now I have it as a factor, but I wonder if that is correct.

3. how to do modelsimplification on the multinom function? Since the summary is a matrix and the drop1 function doesn’t work for multinom, I don’t know how to get to the minimal model. Maybe the most straightforward thing to do is start with the most simple model and add terms and see if the AIC decreases?

I will sent you my data file by email.

This is what I did:

library(MASS)
library(nnet)
setwd("n:\\Tessa\\experimenten\\SpatialdistrIdNc\\R")
mites=read.table("NcIdtryout.txt", header=T)
attach(mites)
names(mites)
head(mites)
model=multinom(patch~replica*idad*nceggs*ideggs, weight=ncad)
summary(model)

Thanks in advance,

Tessa

Statistics Forum Functional Ecology Research School