आर

Question

में एक वेबुल लिंक फ़ंक्शन के साथ मॉडलिंग डेटा मैं सिग्मोइड वक्र संबंधों का पालन करने वाले कुछ डेटा को मॉडल करने की कोशिश कर रहा हूं। मेरे काम के क्षेत्र (मनोविज्ञान) में, एक वेबुल फ़ंक्शन आमतौर पर प्रोबिट के बजाय ऐसे रिश्तों को मॉडल करने के लिए उपयोग किया जाता है।

मैं आर का उपयोग कर एक मॉडल बनाने की कोशिश कर रहा हूं और वाक्यविन्यास के साथ संघर्ष कर रहा हूं। मुझे पता है कि मुझे फ़ंक्शन को VGAM पैकेज से उपयोग करने की आवश्यकता है, लेकिन मैं एक समझदार मॉडल प्राप्त करने में असमर्थ हूं।

# Data frame example data 
dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L)) 

यहाँ dframe1 में डेटा की एक साजिश है::

library(ggplot2) 

# Plot my original data 
ggplot(dframe1, aes(independent_variable, dependent_variable)) + geom_point() 

यह एक वेइबुल समारोह द्वारा तैयार किया जा पाएंगे, क्योंकि डेटा एक फिट चाहिए यहाँ मेरी डेटा है सिग्मोइड वक्र संबंध।

library(VGAM) 

# Generate model 
my_model <- vglm(formula = dependent_variable ~ independent_variable, family = weibull, data = dframe1) 

# Create a new dataframe based on the model, so that it can be plotted 
model_dframe <- data.frame(dframe1$independent_variable, fitted(my_model)) 

# Plot my model fitted data 
ggplot(model_dframe, aes(dframe1.independent_variable, fitted.my_model.)) + geom_point() 

आप देख सकते हैं, यह मेरा मूल डेटा सब पर प्रतिनिधित्व नहीं करता है: यहाँ डेटा मॉडल और एक प्रतिनिधि साजिश उत्पन्न करने के लिए मेरे प्रयास है। मैं या तो अपना मॉडल गलत तरीके से उत्पन्न कर रहा हूं, या मैं मॉडल की अपनी साजिश गलत तरीके से उत्पन्न कर रहा हूं। मैं क्या गलत कर रहा हूं?

नोट: मैंने इसे और अधिक समझने योग्य बनाने के लिए इस प्रश्न को संपादित किया है; पहले मैं पूरी तरह से गलत फ़ंक्शन का उपयोग कर रहा था (weibreg())। इसलिए, नीचे दी गई कुछ टिप्पणियां समझ में नहीं आ सकती हैं। .....

Answer 1

यहाँ मेरी समाधान है, bbmle साथ।

डाटा:

dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L)) 

का निर्माण एक संचयी वेइबुल कि 0.5 से परिभाषा से 1.0 को जाता है: द्विपद बदलाव के साथ,

wfun <- function(x,shape,scale) { 
    (1+pweibull(x,shape,scale))/2.0 
} 

dframe2 <- transform(dframe1,y=round(40*dependent_variable),x=independent_variable) 

एक वेइबुल (लॉग पैमाने पर प्रासंगिक पैरामीटर) फ़िट:

library(bbmle) 
m1 <- mle2(y~dbinom(prob=wfun(exp(a+b*x),shape=exp(logshape),scale=1),size=40), 
    data=dframe2,start=list(a=0,b=0,logshape=0)) 

भविष्यवाणियां उत्पन्न करें:

pframe <- data.frame(x=seq(-0.2,0.3,length=101)) 
pframe$y <- predict(m1,pframe) 

png("wplot.png") 
with(dframe2,plot(y/40~x)) 
with(pframe,lines(y/40~x,col=2)) 
dev.off() 

Answer 2

आप ड्रैक-पैकेज (खुराक-प्रतिक्रिया-मॉडलिंग) का भी उपयोग कर सकते हैं।

मैं वास्तव में मॉडल के इस प्रकार के लिए एक noob हूँ, लेकिन यह किसी भी तरह मदद करता है perhabs ...

यहाँ मैं फिट एक चार पैरामीटर वेइबुल, asymptotes (के लिए तय मानकों के साथ अन्यथा ऊपरी अनंतस्पर्शी थोड़ा अधिक होगा 1, यह नहीं पता कि यह आपके लिए एक मुद्दा है)। मुझे अभिसरण समस्याओं की वजह से स्वतंत्र परिवर्तनीय (+0.2) को भी बदलना पड़ा ताकि यह = = 0 हो।

require(drc) 
# four-parameter Weibull with fixed parameters for the asymptote, added +0.2 to IV to overcome convergence problems 
mod <- drm(dependent_variable ~ I(independent_variable+0.2), 
      data = dframe1, 
      fct = W1.4(fixed = c(NA, 0.5, 1, NA))) 

# predicts 
df2 <- data.frame(pred = predict(mod, newdata = data.frame(idenpendent_variable = seq(0, 0.5, length.out=100))), 
        x = seq(0, 0.5, length.out=100)) 

ggplot() + 
    geom_point(data = dframe1, aes(x = independent_variable + 0.2, y = dependent_variable)) + 
    geom_line(data = df2, aes(x = x, y = pred)) 

हालांकि मैं बेन बोल्कर से सहमत हूं कि अन्य मॉडल बेहतर अनुकूल हो सकते हैं।

मैं केवल इकोटॉक्सिकोलॉजी (खुराक-प्रतिक्रिया-मॉडल से इस तरह के मॉडल जानता हूं, जहां कोई एकाग्रता में रूचि रखता है जहां हमारे पास 50% मृत्यु दर है [= ईसी 50])।

अद्यतन एक चार पैरामीटर लॉग-रसद मॉडल भी काफी अच्छा (छोटे AIC और RSE तो वेइबुल) फिट बैठता है: फिर मैं यहाँ अनंतस्पर्शी पैरामीटर तय की और चतुर्थ बदल दिया।

# four-parameter log-logistic with fixed parameters for the asymptote, added +0.2 to IV to overcome convergence problems 
mod1 <- drm(dependent_variable ~ I(independent_variable+0.2), 
      data = dframe1, 
      fct = LL2.4(fixed=c(NA, 0.5, 1, NA))) 
summary(mod1) 

# predicts 
df2 <- data.frame(pred = predict(mod1, newdata = data.frame(idenpendent_variable = seq(0, 0.5, length.out=100))), 
        x = seq(0, 0.5, length.out=100)) 

ggplot() + 
    geom_point(data = dframe1, aes(x = independent_variable + 0.2, y = dependent_variable)) + 
    geom_line(data = df2, aes(x = x, y = pred)) 

Answer 3

ठीक है, मैं सिर्फ इस कई महीनों देर भर में आया था, लेकिन आप यह भी GLM साथ psyphy पैकेज से mafc.cloglog लिंक का उपयोग कर सकता है। यदि x क्लोग्लॉग का पालन करता है तो लॉग (x) एक कमजोर मनोचिकित्सक फ़ंक्शन का पालन करेगा। उपर्युक्त प्रतिक्रियाओं के साथ पकड़ है कि आपको अनुपात के लिए परीक्षणों की संख्या की आवश्यकता है। मैंने इसे 100 तक सेट किया है, इसलिए यह परीक्षणों की एक पूर्णांक संख्या प्रदान करेगा, लेकिन आपको इसे वास्तव में उपयोग किए जाने वाले नंबरों के अनुरूप ठीक करना चाहिए। यह करने के लिए कोड यहाँ है।

dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L)) 

library(psyphy) 

plot(dependent_variable ~ independent_variable, dframe1) 
fit <- glm(dependent_variable ~ exp(independent_variable), 
    binomial(mafc.cloglog(2)), 
    data = dframe1, 
    weights = rep(100, nrow(dframe1))) # assuming 100 observations per point 
xx <- seq(-0.2, 0.3, len = 100) 
pred <- predict(fit, newdata = data.frame(independent_variable = xx), type = "response") 
lines(xx, pred)