[R-es] ayuda con dos topicos
Kjetil Halvorsen
kjetil1001 en gmail.com
Sab Jul 18 06:13:39 CEST 2009
Mire abajo.
2009/7/17 Marcia Muñoz <marcarmu en gmail.com>:
> Gracias por la respuestas!
>
> Yo estoy usando library(geoR) y tengo el semivariograma ya, el problema es
> que no me lo da ajustado
>
> Los comandos para hacer el ajuste los reconozco (e.g. corSpher o corRatio),
> no se como escribir la formula...ahi es donde me siento perdida.
¿Puedes dar más detalles de lo que has hecho? específicamente, yo no encuentro
corSpher o corRatio en la ayuda de geoR. En geoR hay varios métodos
para hacer variografía, mire
library(geoR)
?variog
?variog4
?variofIt
?eyefit
los últimos ajusta variograms teóricos a un variograma empírico, en
comparación,
?likfit
ajusta parametros en el variograma como parte de una estimación de un
modelo (gaussiano)
para el superficie, junto con trend parametros o otros. Esto puede
usar mucho más tiempo,
pero es teóricamente preferible. Haz ambos, y comparar resultados!
Unos ejemplos:
library(geoR)
library(help=geoR)
data(parana)
?parana
plot(parana, bor=borders)
> ?variog
> vario.b <- variog(parana)
variog: computing omnidirectional variogram
> plot(vario.b)
> vario.c <- variog(parana, op="cloud")
variog: computing omnidirectional variogram
> plot(vario.c)
> vario.s <- variog(parana, op="sm")
variog: computing omnidirectional variogram
> plot(vario.s, ylim=c(0, 2.5))
> ?variofit
> fit.1 <- variofit(vario.b, cov.model="matern", k=0.5)
variofit: weights used: npairs
variofit: minimisation function used: optim
variofit: searching for best initial value ... selected values:
sigmasq phi tausq kappa
initial.value "6258.48" "285.92" "0" "0.5"
status "est" "est" "est" "fix"
loss value: 7287615843.70907
Warning message:
In variofit(vario.b, cov.model = "matern", k = 0.5) :
initial values not provided - running the default search
> fit.2 <- variofit(vario.b, cov.model="matern", k=1.5)
variofit: weights used: npairs
variofit: minimisation function used: optim
variofit: searching for best initial value ... selected values:
sigmasq phi tausq kappa
initial.value "4693.86" "95.31" "625.85" "1.5"
status "est" "est" "est" "fix"
loss value: 7660811461.0475
Warning message:
In variofit(vario.b, cov.model = "matern", k = 1.5) :
initial values not provided - running the default search
> ?cov.spatial
> fit.3 <- variofit(vario.b, cov.model="gaussian")
variofit: weights used: npairs
variofit: minimisation function used: optim
variofit: searching for best initial value ... selected values:
sigmasq phi tausq kappa
initial.value "6258.48" "285.92" "0" "0.5"
status "est" "est" "est" "fix"
loss value: 4933397738.6335
Warning message:
In variofit(vario.b, cov.model = "gaussian") :
initial values not provided - running the default search
> ?likfit
> likfit(parana, ini.cov.pars=fit.1, cov.model="matern")
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
likfit: estimated model parameters:
beta tausq sigmasq phi
" 247.3" " 348.2" "3643.2" " 560.0"
Practical Range with cor=0.05 for asymptotic range: 1677.720
likfit: maximised log-likelihood = -672.2
>
Kjetil
>
> Gracias por su ayuda
>
> saludos
>
> MCM
>
> El 17 de julio de 2009 18:11, Kjetil Halvorsen <kjetil1001 en gmail.com>
> escribió:
>>
>> Hola!
>>
>> Para dar buen ayuda necesitas en realidad explicar algo más lo que
>> quieres hacer.
>> la cosa es que , en R, o mejor, en CRAN, hay muchas maneras alternativas
>> para
>> estimar (semi)variogramas.
>>
>> Un paqute útil es gstat
>> install.packages("gstat", dep=TRUE)
>> library(gstat)
>> library(help=gstat)
>> vignette("gstat")
>> ?fit.variogram
>>
>>
>>
>> Una alternativa interesante, especialmente si vas a trabajar con GLM's
>> en el contexto
>> espacial, es geoR:
>> install.packages(c("geoR","geoRglm"), dep=TRUE)
>>
>>
>>
>> Si vas a trabajar con conjuntos de datso muy grandes (más que un poco
>> de centenares
>> de puntos), puede ser de auyda fields:
>> install.packages("fields", dep=TRUE)
>>
>> En adición hay otras alternativas!!
>>
>> Kjetil
>>
>> 2009/7/17 <PabloEmilio.Verde en uni-duesseldorf.de>:
>> > Hola,
>> >
>> > El paquete "nlme" tiene implementado algunas estructuras de
>> > covariancia espaciales.
>> >
>> > Por ejemplo si utilizas modelos lineales con effectos aleatorios y
>> > estructura de correlacion espacial podes usar "lme". Si usas modelos
>> > lineales sin effectos aleatorios podes usar "glm".
>> >
>> > Esto es un ejemplo editado desde el help:
>> >
>> > # comienzo del ejemplo ...
>> > fm1 <- lme(weight ~ Time * Diet, BodyWeight,
>> > random = ~ Time)
>> >
>> > plot(Variogram(fm1, form = ~ Time))
>> >
>> > fm2 <- update(fm1, weights = varPower())
>> >
>> > fm3 <- update(fm2,
>> > correlation = corExp(form = ~ Time))
>> >
>> > fm4 <- update(fm3,
>> > correlation = corSpher(form = ~ Time))
>> >
>> > anova(fm1, fm2, fm3, fm4)
>> >
>> > # fin ...
>> >
>> > Si miras los "Task views" de R vas a encontrar mucho mas
>> > informacion sobre estadistica espacial.
>> >
>> > Los "bins" son el numero de clases en un histograma. No
>> > tengo "The R Book" asi que no se en que contexto usan la
>> > palabra "bins".
>> >
>> > Saludos,
>> >
>> > Pablo
>> >
>> >
>> >> Saludos a todos!
>> >>
>> >> Estoy haciendo un semivariograma y me gustaria hacer un ajuste al
>> >> semivariograma, yo quisiera probar varios, uno gausiano, el lineal y
>> >> uno
>> >> exponencial para ver cual es mejor para mis datos. No se cuales son los
>> >> comandos en R para hacer eso, alguien tiene un scrip que me pueda
>> >> servir de
>> >> referencia?
>> >>
>> >> La otra pregunta es sobre lo que significa bins?. Estoy siguiendo el
>> >> libro
>> >> de R (The R Book) para hacer los semivariogramas, pero en la pagina
>> >> 776,
>> >> aparece esto: the trend surface ..., followed by the numbers of
>> >> bins....y en
>> >> el libro aparece un 300 que no se de donde sale. Quizas alguno de
>> >> ustedes
>> >> tenga el libro y pueda darme una ayuda con esto, para entender mejor lo
>> >> que
>> >> se hizo en el ejemplo y ver como lo hago con mis datos
>> >>
>> >> muchas gracias!!
>> >>
>> >> saludos
>> >>
>> >> MCM
>> >>
>> >>
>> >> --
>> >> Marcia C. Muñoz
>> >> Universidad de Puerto Rico
>> >> Departamento de Biologia - FB246
>> >> POBOX 23360
>> >> San Juan, PR 00931-3360
>> >> tel: 787-7640000, 2900
>> >>
>> >> [[alternative HTML version deleted]]
>> >>
>> >>
>> >
>> > _______________________________________________
>> > R-help-es mailing list
>> > R-help-es en r-project.org
>> > https://stat.ethz.ch/mailman/listinfo/r-help-es
>> >
>>
>>
>>
>> --
>> "... an entire human genome would fit on a music CD."
>>
>> --- www.thinkgene.com
>
>
>
> --
> Marcia C. Muñoz
> Universidad de Puerto Rico
> Departamento de Biologia - FB246
> POBOX 23360
> San Juan, PR 00931-3360
> tel: 787-7640000, 2900
>
--
"... an entire human genome would fit on a music CD."
--- www.thinkgene.com
Más información sobre la lista de distribución R-help-es