title: Segregation Analysis for predicting age-specific cancer risk

text: I am relatively new to the worlds of bioinformatics and genetics research. I have been tasked with presenting to my lab the potential value of a paper that uses Complex Segregation Analysis for a risk analysis. The hope is that we could adopt this methodology to study our own set of cancer genes; unfortunately I have never heard of segregation analysis before two weeks ago and am hoping someone one here can help. My basic understanding of segregation analysis is thus:   1. Start with pedigrees of the phenotype you're interested in ( _in this paper, age of cancer diagnosis is the phenotype_ )   2. Come up potential explanations of the phenotype's inheritance pattern or lack thereof; these are called "models" ( _and the source of most of my confusion about this paper_ )   3. Run computer programs specifically designed for testing which of your models best fits the pedigree data based on likelihood estimates Within the context of this paper, I am confused about the following: **What models did they actually use?** On page 499, they mention using "two main models of genetic susceptibility" and throughout the paper and supplemental materials there are mentions of models that imply a greater number and variety than the "two main models" described. **Was the risk assessment part of the segregation analysis?** My reading of the paper is that the segregation analysis is providing support for the claim that PALB2 mutations are associated with a given age-specific risk factor for breast cancer; and that the various models used are mainly variations of a single equation with different values and coefficients plugged in. But I am not sure I am reading that right.


title: What is a "(in)dependent p-value"?

text: Does someone know what "independent p-values" or "p-values independently distributed" means? And the "dependence structure among test statistics"? Thanks in advance


title: Backing out the standard deviation from information on baseline mean/s.d., and coefficient mean/s.d

text: I am trying to run a power calculation for a randomized control trial. For this I need a mean and standard deviation for our 'baseline'. There are papers out there which would have a mean and standard deviation, but what we would want would be to use their final outcomes as our baseline. In other words, these papers applied treatment $A$. I would like to apply treatment $B$ on top of treatment $A$, so treatment $A$ outcomes will be my new baseline. In looking at their results tables, I can find the mean and SD of their baseline, plus the coefficient of the effect and that coefficient's SD. To find the new mean of my baseline, I would add the coefficient to the mean. However, is there a way to back out the SD of the new mean from the information I have from these papers? To give an example: one paper has the mean in control as 5847 (SD=16784). The change of the mean under treatment $A$ is 267 (SD=527). So my new baseline would be 5847 + 267 = 6114. But would there be any way of getting this SD without the data? (Almost all of it is not available publicly.) Thanks so much!


title: Dealing with missing data in an exponential smoothing model

text: There does not seem to be a standard way to deal with missing data in the context of the exponential smoothing family of models. In particular, the R implementation called _ets_ in the _forecast_ package seems to just take the longest subsequence without missing data, and the book "Forecasting with Exponential Smoothing" by Hyndman et al. does not seem to talk about missing data at all. I would like to do a little more, if my users explicitly ask me to (and if the missing data don't occur too close together or in too many periods that are exactly one season apart). In particular, what I have in mind is the following. During simulation, whenever I would encounter a missing value $y_t$, I would substitute the current point forecast $\tilde y_t$ for $y_t$, so that $\varepsilon_t = 0$. This would, for example, make the data point not be considered for the parameter optimization process. Once I have a reasonable fit for the parameters, I can estimate the standard deviation of the errors (assumed to be normal with mean $0$) and verify that using values for $\epsilon_t$ generated from that distribution do not decrease the likelihood by a big factor. I would use such values for forecasting (using simulation) as well. Are there known pitfalls with this method?


title: Is the Dynamic Time Warping constraint Itakura Parallelogram based on series length?

text: I'm looking into the Itakura Parallelogram global constraint for Dynamic Time Warping. I am confused about the maximum width of the parallelogram, is it solely based on the length of the two time series being compared, or can it be varied to achieve a bigger or smaller warping window?


title: Estimation of white noise parameters in Gaussian random walk model

text: I want to estimate the parameters (mean , variance ) $e(t)$ for the random walk model   $X (t) = X (t-1) + e(t)$. (where $e(t)$ is the white noise with a Normal distribution). By using the fact that the $E(X(t)) = E (X(t-1))$ mean can be easily proven to be zero. Can someone help me to come up with a formal proof for the estimation of std. deviation/variance? Intuitively my guess is that variance should be defined as : $\sum\frac{(X(t) - X(t-1))^2}{N-1}$ but I don't know how to formally prove it.


title: Marginal distribution of a function of order statistics

text: From the joint distribution of any two order statistics, say $Y_j$ and $Y_k$, $j<k$ I would like to derive the distribution of $Z=F(Y_k)-F(Y_j)$. The initial pdf is: $$f_{Y_j,Y_k} (y_j,y_k) =\frac{n!}{(j-1)!(k-1-j)!(n-k)!} \times [F(Y_j)]^{j-1} [F(y_k)-F(y_j)]^{k-1-j} \times[1-F(y_k)]^{n-k} f(y_j) f(y_k)$$ Now I make the substitutions $Z=F(Y_k)-F(Y_j)$ and $ Y=F(Y_j)$ and the strategy is to get joint distribution of $Z$ and $Y$ and inegrate out $y$ to get the marginal distribution of $Z$. Now the joint distribution is $$f_{Z,Y} (z,y)=\frac{n!}{(j-1)!(k-1-j)!(n-k)!} z^{k-1-j} y^{j-1}(1-z-y)^{n-k} $$ and unless I am mistaken the support is $\{0<z,y<1 \}$, yes? If that is the case, now I need to evaluate the integral $$f_Z (z)=\int_{0}^{1} \frac{n!}{(j-1)!(k-1-j)!(n-k)!} z^{k-1-j} y^{j-1}(1-z-y)^{n-k}dy$$ This is the part where I would appreciate some help. The integral looks like a beta pdf but that $z$ in the integrand confuses me. Are there any tricks here? Thank you.


title: Proportional odds assumption violated but not sure what to do

text: I run an ordinal regression model and I wanted to check the proportional odds assumption. In order to do that I used the `VGAM` package and I run olr twice, the first under the assumption and the second without the assumption. Below is the code and the results               > fit1 <- vglm(stage ~Ki67+Cyclin_E,family=cumulative(parallel=T))     > summary(fit1)          Call:     vglm(formula = stage ~ Ki67 + Cyclin_E, family = cumulative(parallel = T))          Pearson Residuals:                          Min       1Q  Median      3Q    Max     logit(P[Y< = 1]) -3.1177 -0.43593 0.37246 0.53111 1.4927     logit(P[Y< = 2]) -3.8479  0.14119 0.18785 0.28679 1.9903          Coefficients:                     Estimate Std. Error  z value     (Intercept):1  2.2414705   1.091225  2.05409     (Intercept):2  3.2164214   1.178916  2.72829     Ki67          -0.1157273   0.039889 -2.90124     Cyclin_E       0.0085266   0.028626  0.29786          Number of linear predictors:  2           Names of linear predictors: logit(P[Y< = 1]), logit(P[Y< = 2])          Dispersion Parameter for cumulative family:   1          Residual deviance: 50.82946 on 62 degrees of freedom          Log-likelihood: -25.41473 on 62 degrees of freedom          Number of iterations: 5                > fit2 <- vglm(grade ~Ki67+Cyclin_E,family=cumulative(parallel=F),maxit=50)     > summary(fit2)          Call:     vglm(formula = grade ~ Ki67 + Cyclin_E, family = cumulative(parallel = F),          maxit = 50)          Pearson Residuals:                          Min       1Q   Median      3Q    Max     logit(P[Y< = 1]) -1.1870 -0.65271 -0.23199 0.44910 3.4798     logit(P[Y< = 2]) -2.6235 -0.70599  0.27305 0.72691 2.8544          Coefficients:                    Estimate Std. Error   z value     (Intercept):1 -0.059702   0.928078 -0.064328     (Intercept):2  2.687277   1.050909  2.557097     Ki67:1        -0.100832   0.047754 -2.111483     Ki67:2        -0.101817   0.036567 -2.784377     Cyclin_E:1     0.018768   0.022708  0.826474     Cyclin_E:2    -0.015416   0.022927 -0.672390          Number of linear predictors:  2           Names of linear predictors: logit(P[Y< = 1]), logit(P[Y< = 2])          Dispersion Parameter for cumulative family:   1          Residual deviance: 78.93318 on 78 degrees of freedom          Log-likelihood: -39.46659 on 78 degrees of freedom          Number of iterations: 34       In order to check if the difference of the two models is significant I run the next command               pchisq(deviance(fit2)-deviance(fit1),df=df.residual(fit2)-df.residual(fit1),lower.tail=FALSE)     [1] 0.03072927      As you see the result is that the 2 models differ and so the proportional odds assumption isn't true. But if you see the coefficients about the Ki67 (Cyclin is not significantly important so i guess i can skip it) they are almost the same. In that case what should I do? I believe that I could stick with the model under the po assumption but I'd like to know what others think


title: How does the Bayes' theorem equation generalize all sorts of regression/classification models?

text: I have been reading “Pattern Recognition & Machine Learning” written by Christopher M. Bishop for some time, but I am still a beginner. I wish to get a bigger view that summarizes regression and classification models. The book begins with Bayes’ theorem stating that: $$ p(\textbf{w} \mid D) = \frac{p(D \mid \textbf{w})p(\textbf{w})}{p(D)} $$ Where $D$ is the observed data, $w$ is model parameters. From this equation, how does it generalize all sorts of regression or classification models? Specifically, I wish to know how it relates to following model types:   * Clustering, e.g. k-means clustering   * Density estimation, e.g. k-nearest neighbours   * Classification, e.g. Gaussian classifier   * Regression, e.g. Bayesian linear regression I know each of the technique, but how does the likelihood and prior probability functions in the Bayesian relate to them?


title: Convert ARIMA Model in R to Java Code (only the forecasting)

text: I would like to convert an ARIMA model developed in R using the `forecast` library to Java code. Note that I need to implement only the forecasting part. The fitting can be done in R itself. I am going to look at the `predict` function and translate it to Java code. I was just wondering if anyone else had been in a similar situation before and managed to successfully use a Java library for the same. Along similar lines, and perhaps this is a more general question without a concrete answer; What is the best way to deal with situations where in model building can be done in Matlab/R but the prediction/forecasting needs to be done in Java/C++? Increasingly, I have been encountering such a situation over and over again. I guess you have to bite the bullet and write the code yourself and this is not generally as hard as writing the fitting/estimation yourself. Any advice on the topic would be helpful.


title: Equivalence of with and without replacement sampling

text: After some heavy reordering and canceling of factorials, I discovered that the following experiment is approximately equivalent for $m \ll n < N$ if conducted with or without replacement: > In $n$ turns, draw marbles (with/without replacement) from an urn containing > $m$ white and $N-m$ black marbles. Count the white marbles. Now, the "with replacement" part is a Bernoulli trial: > Throw $m$ times an unfair coin with success probability $n/N$. Count the > number of successes. Are there any pitfalls in this approximation? I am pretty sure that this is textbook knowledge. Can you recommend a good book where this approximation is derived? I did the exercise in order to Derive househould weights from a uniformly distributed person sample.


title: What should I call these growth rates?

text: I have a table of sales growth rates by month calculated thus: $$\text{growth rate}_i=\dfrac{\text{sales in month } i-\text{sales in January}}{\text{sales in January}}$$ What should I call this table? The best I've come up with it "sales growth rates from January by month". Is this ambiguous? Is there a technical term for this? (I've had almost nothing to do in my life with anything that counts money or statistics for that matter).


title: Apply an exponential weighted moving average to data

text: I have series of returns on a stock. I want to apply exponential smoothing to my data using a rolling window with R. I don't how to do that. I define a function this way:               weight <- function(data) {       result1 <- (rep(0.94,36)^rev(c(0:35))       result2 <- result/sum(result1)      and I do `rollapply` for:               data.rdt <- rnorm(300,0,1)          roll <- roolapply(data.rdt,width=36,FUN=weight,by.column = FALSE,by=1)      but it doesn't work. Please help me.


title: Growth mixture models, "structural parameters", and using a parameter with variance fixed at zero as mediator in regression model?

text: I am examining a paper that uses growth mixture modeling. They have estimated intercept, linear slope, and quadratic growth parameters. They have also estimated random effects for the intercept and linear slope. However, they have set the variance of the quadratic slope to zero (which I understand to be using fixed effects for that parameter). They have stated that because they fixed the variance of the quadratic slope to zero, this parameter could not be used in the 'structural part of the model'. I am just wanting to learn a bit about why this may be. I presume that 'structural parameters' are just normal parameters of the model? Is there a better definition for this? What are non-structural parameters? Also, is it correct that this parameter cannot be used in regression analyses because its variance is fixed at zero? I would assume that you could used fixed effects parameters in regression? Any guidance would be much appreciated. Cheers! -George


title: Minimum sample size to observe certain outcomes

text: Given a categorical random variable with large number of possible outcomes and a sufficiently large sample, I observe that 90% of the sample falls into some relatively small number of categories, say, 30 out of 1000. How can I estimate minimum sample size such that these 30 categories are observed? A sort of a solution is, I guess, to produce smaller datasets resampling with replacement from the original one.


title: Power of chi-square test for goodness-of-fit as function of sample size

text: I would like to calculate the power of a Chi-Square test for goodness-of-fit as a function of sample size for a specified alpha-value (say 0.01). Specifically, I am referring to power as the probability that the test will correctly reject the null hypothesis. In most places I look, I can find only vague references that this can be done and/or am referred to software. It seems to me that this should not be difficult. The specific question I would like to be able to answer (without using some sort of opaque software) follows: Given an $\alpha$ value, for what sample size, $n$, can I expect a particular power, $\beta$? Even better, how might I be able to calculate an ROC?


title: How to interpret the regression slope of standardised variables

text: my y variable is crop production (in kg/ha) and x variable is rainfall (mm). In a linear regression, if the slope is 0.5. I would say, 'if rainfall increases by 1 mm yield increases by 0.5 kg/ha.' Now if I do the regression again but this time I have standardised by x and y variable by centering and dividing by 2 std. dev., how do I say the previous statement 'if rainfall increases by -------- yield increases by ------.' Thanks


title: Goodness of fit for nonlinear model

text: We have fitted a nonlinear function to observed data. The next step should be the assessment of the goodness of fit of this function (like $R^2$ for linear models). What are the usual ways to measure this? **Edit 1:** The fitting was performed as follows:   1. Perform a linear regression with independent variables _A_ and _B_.   2. Calculate distribution's parameters from regression parameters. (The distribution is nonlinear and has variable _C_ as an input.)   3. _Assess goodness of fit of nonlinear distribution by comparing estimated to observed data._ **Edit 2:** Examples for the steps mentioned above:   1. Regression model: $log(y) = \beta_0 + \beta_1 \centerdot log(a) + \beta_2 \centerdot log(b)$   2. $\rho = -\frac{\beta_0}{\beta_1}$ and $\theta = \beta_2$ for the following nonlinear distribution: $f(a) = \rho \centerdot a^{-\theta}$   3. _Assess the goodness of fit of $f(a)$ with a given set of $(a, f(a))$ observations._


title: Predict future student outcomes (binary and continuous) with historic cross-sectional data?

text: Using Stata 11.2, I would like to develop 2 analytic models that could be implemented by school administrators to flag students for intervention. I'm wondering if it would be possible to develop these models based on historical cross-sectional data consisting of 650,000 unique observations (11th grade students); it was collected from 2009 through 2011. First, I want to predict risk of dropping out of school (with the goal of intervening in 5% of students with the highest risk of dropout). Second, I want to predict absence (with the goal of intervening among students in the top 5% of predicted absentee time). Outcome variables: dropout (binary; 0=no, 1=yes) and absence (continuous; cumulative hours absent from school) Predictor variables: sex (binary), track (3-categories), GPA (continuous) So far, I've just worked on predicting dropout, but am not sure if what I'm doing is theoretically correct or statistically sound. (It's been a long time since I took stats!) Here what I've done, using the Stata command:               logistic dropout i.sex i.track gpa      The output shows all independent variables are significantly associated with dropout, based on p<0.05 and 95% CIs for the ORs that do not contain 0. Sex (OR=0.95), Track (OR=0.76), GPA (OR= 1.52). Now I'm not sure about how to proceed in calculating the predicted probabilities of dropout, and then figuring out which students are at greatest risk of dropout. Should the command be `predict phat` ? I think this gives the predicted prob of dropout for each level of each variable, holding all other variables at their mean. Then I would just categorize the predicted values into 20 categories, and the students in the top category would be those who should be targeted for intervention? I would greatly appreciate recommendations on how to conduct both analyses (binary and continuous outcomes).


title: Idea of errors as 'shocks' to a regression?

text: So I always understood the error term to be the difference between the observed value from the true, yet unobservable, function value. However, I often here, especially related to the economics fields, articles and people talking about error terms as 'shocks'. Can anyone offer an intuitive, non-math intensive, explanation?


title: R package fitdistrplus: function fitdist shows unusual error

text: I used the package R fitdistrplus in order to estimate parameters of different distribution function. In one case it shows a strange error message I can't interpret:               > fitdist(data=rf,distr="gev",method="mle",start=list(loc=0,scale=1,shape=0))     Fitting of the distribution ' gev ' by maximum likelihood      Parameters:            estimate Std. Error     loc   -1800.9939106        NaN     scale  2371.4799521        NaN     shape    -0.4909284 0.03581115     Warning messages:     1: In sqrt(diag(varcovar)) : NaNs produced     2: In sqrt(1/diag(V)) : NaNs produced     3: In cov2cor(varcovar) :     diag(.) had 0 or NA entries; non-finite result is doubtful      For another dataset, the command works. Is it maybe due to the length of the series (187)?


title: BRT predictions on zero-inflated gaussian fish abundances include negative results

text: hopefully someone can point me in the right direction here. I'm using boosted regression trees (BRT) to assess the relative importance of a number of environmental factors (sea bottom temperature, sea bottom salinity, substrate grain size, depth, distance from shore, maximum water speed @ seabed) on fish abundances (4 different species of rays, plus all 4 summed) from 1440 sample sites in the Irish Sea, then creating a predicted surface of abundance estimates for a grid of all of those environmental variables across the whole Irish Sea. Concentrating on all rays summed, the range of abundances goes from zero (loads of these: zero inflated) to 126. My problem is that the range of predictions produced by the model goes from -7.06 to 36.38, which assumedly has to be wrong? For the BRT I simply used the count data, with distribution="gaussian"... but I'm yet to find any other examples of people doing this approach. Elith et al 2008, who built the code expanding on Ridgeway's original BRT work, reduce their data to presence/absence and use binomial (Bernoulli). In his thesis comparing BRTs to GAms & GLMs, Abeare uses gaussian, but only after removing the zeroes (i.e. zero truncated). I'm wondering if anyone might know why the predictive results would be negative? Or if anyone could recommend a 'best' way of proceeding? My supervisor suggested lognormally transforming the data, as Abeare did, but this simply produces smaller negatives... Thanks in advance!


title: Outlier detection using regression

text: Can regression be used for out lier detection. I understand that there are ways to improve a regression model by removing the outliers. But the primary aim here is not to fit a regression model but find out out liers using regression


title: GLM for count data

text: I ran an experiment with an eye tracker and my data frame has this look:                              Condition   DWellsAOI1 DwellsAOI2  TotalDwells      Participant1       1             12         13            25     Participant2       2            100         11           111     Participant3       1             50         50           100      and so on. `DWellsAOI1` counts the number of dwells on AOI1. Each participant belongs to one condition only, and the duration of the experiment is not fixed, so different `TotalDwells`. I was trying to check if there is a significant difference among the between conditions in terms of `DwellsAOI1`, and my first approach was to compute the percentages (`DwellsAOI1/TotalDwells`) and run anova. However, data violates both the assumptions. I searched on Google and I found that I can use a generalized linear model (GLM) with binomial family. Currently I'm running it like this:               mod <- glm(cbind(DwellsAOI1,TotalDwells-DwellsAOI1) ~ Condition,                    data=df, family=binomial("logit"))     summary(aov(mod))      Is it the correct way? Thanks!


title: Neural net model - error during training

text: I'm getting started with R, I really like it but recently I found myself in a corner. I'd like to build neural network model that predicts heat consumption. I have historical data that contains outside air temperature (model input) and heat demand values (model output) in megawatts (hourly data from past 4 years). I would like to use my model to predict heat demand for 24h ahead based on air temperature weather forecast (also 24 hours ahead). Here is my code:               > data <- read.delim("C:/.../data.csv", dec=",")     > require(neuralnet)     > trainset<-data[1:26208,]     > testset<-data[26209:26232,]     > net<-neuralnet(heat~temp,trainset,hidden=5,threshold=0.01)      .. and I get error that 'algorithm did not converge' This is my first attempt to build model. That's why I want to use only one input parameter (air temperature), In the future, I would like to use few more inputs like wind speed, info about working days and holidays and so on. Do you have any idea what I did wrong? Is it the problem with number of hidden neurons or layers? I also tried to use other values of parameter 'hidden' and I still got errors. Here is my data: click Thank you in advance for your help.


title: Given two independent normal random variables $X$ and $Y$, what is $P(X\leq x\mid X>Y)$?

text: As the title says, I'm looking for the distribution of $X$ given that $X>Y$.


title: Correlation for grouped data

text: I have a dataset composed by 200.000 projects. Each project is defined according to its size and the presence of active users (1 if the project has active users, 0 otherwise). I've ordered the projects according to their sizes (ascending order), then I've grouped the projects in 10 groups (20.000 projects per group), and for each group, I've counted the number of projects having active users. What I have now this: `group|avg project size| #projects with active_users 1 58.7 1685 2 76.5 1804 3 86.1 1808 4 98.4 1893 5 114.8 1886 6 137.6 1949 7 172.1 1922 8 231.2 1999 9 373.2 1949 10 2698.4 1524` I would like to check if a correlation exists between the average size of projects in a group and the number of projects having active users. Which test it's better to use in this case? Finally, if I found a correlation, how powerful would be this correlation with respect to a correlation found on the raw data? I would appreciate any kind of help/suggestion thanks in advance


title: T test assumptions

text: I'm performing a t-test on a time series with a sliding window (i.e. every N samples, perform a t test). I know that overall, the samples are roughly normally distributed, however adjacent samples are correlated, so the N samples (taken for the t test) may not be normally distributed. Are the t test assumptions violated in this case?


title: How to Build a Foresight System?

text: For a research project, I'm asked to find ways to build an economic foresight system. For example, for the production of cheese. We will have data about the market indicators, like price, demand etc. And we want to build a system that estimates the future state of these variables. Machine Learning and Bayesian Networks might be an approach, but I think they mostly assume linearity. I see that complex systems is a relevant approach, but I don't know where to start. All the papers about complex systems I found are verbal, instead of mathematical / implementation etc. Can you suggest some readings to get an idea for how to construct a system with given parameters, estimates the future ? Thanks for any help!


title: Significance factors

text: How can I determine a significance factor? For example, If 1/1 samples returns positive, the ratio is one. However, if 9/10 samples return positive, the ratio is .9, but is better because of the number of samples tested. Is there any way to statistically rank these test scenarios. I know it may be up to subjectivity, especially when deciding harder ones, such as 9/10 and 85/100, but any approach would be appreciated.


title: How exactly do neural network and support vector machine methods reduce dimensionality?

text: Please, how artificial neural network and support vector machine methods reduce dimensionality of variables? And based on what do we select the extraction features (retained variables)?


title: Why break down the denominator in Bayes' Theorem?

text: (I'm a newbie at stats. I'm a mathematician and a programmer and I'm trying to build something like a naive Bayesian spam filter.) I've noticed in many places that people tend to break down the denominator in the equation from Bayes' Theorem. So instead of this: $\frac{P(A|B)\cdot P(B)}{P(A)}$ We are presented with this: $\frac{P(A|B)\cdot P(B)}{P(A|B)\cdot P(B)+P(A|\neg B)\cdot P(\neg B)}$ You can see that this convention is used in this Wikipedia article and in this insightful post by Tim Peters. I am baffled by this. Why is the denominator broken down like this? How does that help things at all? What's so complicated about calculating $P(A)$, which in the case of spam filters would be `The probability that the word "cheese" appears in an email, regardless of whether it's spam or not`?


title: Doing regressions on samples from a very large file: are the means and SEs of the sample coefficients consistent estimators?

text: I have a fairly larege file 100M rows and 30 columns or so on which I would like to run multiple regressions. I have specialized code to run the regressions on the entire file, but what I would like to do is draw random samples from the file and run them in R. The strategy is: randomly sample N rows from the file without replacement run a regression and save the coefficients of interest repeat this process M times with different samples for each coefficient calculate the means and standard errors of the coefficeints over M runs. I would like to interpret the mean computed over M runs as an estimate of the values of the coefficients computed on the whole data set, and the stadard errors of the means as estimates of the standard errors of the coefficients computed on the entire data set. Experiments show this to be a promising strategy, but I am not sure about the underlying theory. Are my estimators consistent efficient and unbiased? If they are consistent how quickly should they converge? What tradeoffs of M and N are best? I would very much appreciate it if someone could point me to the papers, books etc. with the relevanth theory. Best regards and many thanks, Joe Rickert


title: Interpreting Log-Transformed Percentages in OLS

text: In a log-log model, such as $\log(y) = b_0 + b_1 \log(x)$, I know that with OLS the standard interpretation is a "1% increase in x is associated with a $b_1$% increase in y." I have three related questions:   1. If x is a percentage variable that was logged to correct skew (ex: x ranges from .01 to .99), what is the correct interpretation of the resulting regression? A 1% increase in x no longer seems easily interpretable in this instance.   2. How does the interpretation change (if at all), if the model is differenced? For example: $\log(y_t) - \log(y_{t-1}) = b_1 \left( \log(x_t) - \log(x_{t-1}) \right)$. This is frequently observed in panel data. To me, this seems to be also modeling a percentage relationship between y and x, and it's unclear if the interpretation would differ.   3. Do native predict functions in either R or Stata handle logarithmic transformations accurately, or do they need to be exponentiated/corrected for bias? To sum up, I would like to know how to accurately generate the predicted percent change in y given a (percentage point) value of x -- in contexts where x may be logged and/or differenced.


title: Fair Ranking of Ratings Based on Number of Votes

text: I'm attempting to take a list of products and rank them according to rating and votes. This is tricky! Lets suppose product **A** is rated at average **5.6** at **310** votes. Product **B** is average **7.3** at **10** votes. Product **C** is **3.4** at 1 vote. Product **B** has hardly earned a winning place with 10 votes, but product **A** is lesser rated with a much bigger voting base, showing it's persistence. Product **C** shouldn't be ranked poorly just because one (possibly ignorant) user is pissed at it. Maybe it's new to the market - much newer than product **A**. I was hoping to find the most reasonable "modifier" to get them to rank according to earned position. Perhaps just a small modifier? `Votes/Rating` to distinguish them?


title: Computing probability distribution function for uniform random variables and Y=1-X

text: X is uniform random variable in [0,1] and Y=1-X. How do I calculate the distribution function F(X,Y)? I can see that Y is also uniformly distributed and can draw the intervals. But I am unable to compute the distribution function when x+y>1 & x,y in [0,1]. How do I compute it for this specific case?


title: How to make Marsaglia polar method return values [0, 1)?

text: I've implemented the Marsaglia polar method to generate random numbers that are normally distributed. Unfortunately the code as shown on this website returns numbers that are not within the range [0, 1). Code:                    double x1, y1, w;               do {                  x1 = 2.0 * ranf() - 1.0;                  y1 = 2.0 * ranf() - 1.0;                  w = x1 * x1 + y1 * y1;          } while ( w >= 1.0 || w == 0 );               w = sqrt( (-2.0 * ln( w ) ) / w );          return x1 * w;      where ranf() returns values [0, 1). For example, if (x1, y1) was (-0.43458512585358, -0.07521858582050478), this returns a value of -1.7830255550765148. Obviously not a value within my expected range. I've seen some implementations that multiplies this by the standard deviation and adds the mean. But if I want to get numbers that range from [0, 1), what should I use as my standard deviation and mean? currently, I'm using a standard deviation of 1/8.125 ~= 0.121212... and 0.5 as the mean, but I've only stumbled on this by experimentation. What is the official way for me to get this appropriately ranged? To answer the question in the comments below: I want to generate normal random variables within the range 0 and 1. Obviously I'd like the mean to be at 0.5 and values no be distributed around this mean. More for my enlightenment, will the polar method really return numbers over the entire real line? If so, then it should be easy for me to scale the result down to [0, 1), assuming I don't hit floating point underflow due to the division to scale things down. New edit: What I'm trying to do is take Aniko's advise in his answer to my other Power Factor Problem question. So my first thought was that simulate my distribution of velocities, was to implement my on GaussianRandom class by deriving from Random and overriding Random.Sample() and the 3 other methods as recommended by MSDN. Since Random.Sample() is suppose to return a double in the range [0, 1), then I've hit this stumbling block. I've come to the same conclusion as mbq's answer after staring at the ln(x) graph where x is in (0, 1] and sqrt(x) where x is in (0, 2]. From a practical/programming standpoint, I can't override Random.Sample(), because later when I actually call GaussianRandom.NextDouble(), there is no way for me to pass in the actual standard deviations for my velocities. It looks like I'll have to use the Marsaglia polar method in the raw instead of wrapped within the Random class. My "cartridge" generator will have to `return x1 * w * stdDev + mean`.


title: Relationship between poisson and exponential distribution

text: The waiting times for poisson distribution is an exponential distribution with parameter lambda. But I don't understand it. Poisson models the number of arrivals per unit of time for example. How is this related to exponential distribution? Lets say probability of k arrivals in a unit of time is P(k) (modeled by poisson) and probability of k+1 is P(k+1), how does exponential distribution model the waiting time between them?


title: Survival analysis in R with left-truncated data

text: I am doing a survival analysis in R with the `survival` package. I think I am working with left-truncated data, but I'm not entirely sure how to handle it. I have a cohort of patients diagnosed between 1990 and 2012. All the patients have a well-defined time of diagnosis (entry time). However, the outcome of interest (specific worsening of disease) has only been documented from the year 2000 and onwards. For patients diagnosed before 2000, it is therefore not known whether the outcome has occurred before that time - only after. My first thought was that I needed to restrict the analysis to the time period from 2000, only including patients diagnosed after that point in time. After doing some reading, it appears to be unnecessary to exclude patients diagnosed before 2000. This seems to be left-truncation and that can be dealt with in `coxph` using `Surv(time1, time2, event)`, where time1 is left-truncation time (time from diagnosis to the start of documentation of the outcome) and time 2 is the time-to-event (from time of diagnosis). Here are two examples of patients in my dataset: Patient #1: Diagnosed in 1999. Outcome observed in 2001. Left-truncation time: 1 year (to 2000). Time-to-event: 2 years. Patient #2: Diagnosed in 2001. Outcome observed in 2005. Left-truncation time: 0 years. Time-to-event: 4 years. For these patients, I suppose their survival times (in years) in the survival object would be (respectively):               Surv(time1 = c(1,0), time2 = c(2,4), event = c(1,1))      Is this an example of left-truncated data? If so, is this the correct way to handle it?


title: Free treemapping software

text: Is there any free desktop software for treemapping? Every one I've found appears to have a commercial license. Many Eyes is free and has treemapping features but what I'm really looking for is downloadable desktop software. Thanks in advance.


title: Variance Covariance for logit with elastic net

text: How do you calculate an estimate for the variance covariance matrix of a logistic regression with elastic net regularization? Starting from the variance-covariance matrix of a plain vanilla logistic regression, how does the formula need to be augmented: $\hat{\Sigma}=-\left[\hat{p}(1-\hat{p})XX^{T}\right]^{-1}$   where $\hat{p}=\frac{1}{1+exp(X^{T}\hat{\Theta})}$   If you don't know the answer for elastic net, how would this be implemented for a logistic regression with ridge regularization? Update: This link suggests OLS ridge is:   M = $(X'X+λI)^{−1}X'$   var(β)= ${\sigma}^{2}MM'$


title: How to compare ordering of two groups with widely different scales

text: I'm trying to compare if the order in col 1 is comparable to the order in col 2. The simplest thing I can do is just to compare both ranks using something like Friedman's test. However, Col1 and Col2 is on a continuous basis. And for instance, the .94 and .93 in column 2 are closer together than 13.06 is to .94. Thus it's not quite ordinal data. What would be a good way to compare the order but while taking into account the magnitude? Thanks!                          Col1         Col2     7          1610635       1.61     8          3214571       2.04     9          2188633       1.30     4          1585069      13.06     5          2009116      14.74     1          2091101      12.50     2          2422071      14.64     3          1633422      13.05     10         3428149       0.93     11         2001198       1.20     12         2191020       0.94     13         2294395       8.83     14         2921134       9.18     15         2659132       6.60        * Background information: As requested, here's the background. Column 1 represents the overall counts of a particular type of T-cell using one type of instrument which reads raw counts. Column 2 represents the number of t-cells found using a different type of instrument which provides things on a percentage basis. We are looking to compare that the two instruments roughly act the same -- they should pick up roughly the same number (rather on a percentage or on a raw count). I don't have any other data -- just these two columns. Thanks!


title: Hierarchical Bayesian Model (?)

text: Please apologize my butchering of statistical lingo :) I have found a couple of questions on here that are related to advertising and click through rates. But none of them helped me very much with my understanding of my hierarchical situation. There's a related question Are these equivalent representations of the same hierarchical Bayesian model?, but I'm not sure if they actually have a similar problem. Another question Priors for hierarchial Bayesian binomial model goes into detail about hyperpriors, but I'm not able to map their solution to my problem I have a couple of ads online for a new product. I let the ads run for a couple of days. At that point enough people have clicked on the ads to see which one gets the most clicks. After kicking out all but the one that has the most clicks I let that one run for another couple of days to see how much people actually purchase after clicking on the ad. At that point I know if it was a good idea to run the ads in the first place. My statistics are very noisy because I don't have a lot of data since I'm only selling a couple of items every day. Therefore it is really hard to estimate how many people buy something after seeing an ad. Only about one in every 150 clicks results in a purchase. Generally speaking I need to know if I'm losing money _on each ad_ as soon as possible by somehow smoothing the per-ad group statistics with global statistics over all ads.   * If I wait until every ad has seen enough purchases I'll go broke because it takes too long: testing 10 ads I need to spend 10 times more money so that the statistics for each ad gets reliable enough. By that time I might have lost money.   * If I average purchases over all of the ads I won't be able to kick out ads that just aren't working as well. Could I use the global purchase rate ($ per click) and use it as a prior for $N$ sub-distributions? That would mean that the more data I have for each ad, the more independent the statistics for that ad get. If nobody has clicked on an ad yet, I assume that the global average is appropriate. Which distribution would I choose for that? If I have had 20 clicks on A and 4 clicks on B, how can I model that? So far I have figured out that a binomial or Poisson distribution might make sense here:   * `purchase_rate ~ poisson` (?)   * `(purchase_rate | group A) ~ poisson` (estimate the purchase rate only for the group A?) But what do I do next to actually calculate the `purchase_rate | group A`. How do I plug two distributions together to make sense for group A (or any other group). Do I have to fit a model first? I have data that I could use to "train" a model:   * Ad A: 352 clicks, 5 purchases   * Ad B: 15 clicks, 0 purchases   * Ad C: 3519 clicks, 130 purchases I'm looking for a way to estimate the probability of any one of the groups. If a group has only a couple of datapoints I essentially want to fall back to the global average. I know a bit about Bayesian statistics and have read a lot of PDFs of people describing how they model using Bayesian inference and conjugate priors and so on. I think there is a way to do this properly but I can't figure out how to model it correctly. I would be super happy about hints that help me formulate my problem in a Bayesian way. That would help a lot with finding examples online that I could use to actually implement this. **Update:** Thanks so much for responding. I'm beginning to understand more and more little bits about my problem. Thank you! Let me ask a few questions to see if I understand the problem a bit better now: So I assume the conversions are are distributed as Beta-distributions, and a Beta distribution has two parameters, $a$ and $b$. The $\frac{1}{2}$ $\frac{1}{2}$ parameters are hyperparameters, so they are parameters to the prior? So in the end I set the number of conversions and number of clicks as the parameter of my Beta distribution? At some point when I want to compare different ads, so I would compute $P(\mathrm{conversion} | \mathrm{ad}=X) = \frac{P(\mathrm{ad}=X | \mathrm{conversion}) * P(\mathrm{conversion})}{P(\mathrm{ad}=X)}$. How do I compute each part of that formula?   * I think $P(\mathrm{ad}=X | \mathrm{conversion})$ is called likelihood, or "mode" of the Beta distribution. So that's $\frac{\alpha - 1}{\alpha + \beta - 2}$, with $\alpha$ and $\beta$ being the parameters of my distribution. But the specific $\alpha$ and $\beta$ here are the parameters for the distribution just for the ad $X$, right? In that case, is it just the number of clicks and conversions this ad has seen? Or is how many clicks/conversions _all_ ads have seen?   * Then I multiply with the prior, which is P(conversion), which is in my case just the Jeffreys prior, which is non-informative. Will the prior stay the same as I get more data?   * I divide by $P(\mathrm{ad})$, which is the marginal likelihood, so I count how often this ad has been clicked? In using Jeffreys' prior, I'm assuming that I'm starting at zero and don't know anything about my data. That prior is called "non-informative". As I continue learning about my data, do I update the prior? As clicks and conversions come in, I have read that I have to "update" my distribution. Does this mean, that the parameters of my distribution changes, or that the prior changes? When I get a click for ad X, do I update more than one distribution? More than one prior?


title: When should a weighted KNN be used (or not)?

text: By default, machine learning packages turn inverse distance weighting off for KNN. To me, it seems that inverse distance weighting is always a good option. Why would we not want to use IDW with KNN? [And why would we want to?] Thanks


title: Does multi-response linear regression introduce artificial correlation among response residuals?

text: This is a multi-response linear regression experiment. I generate a random matrix of gaussian noise $\mathcal{N}(0,1)$ sized 10x6 called $Y_N$. I then generate an independent response $X$ via either uniform or gaussian noise, 10x1. I then regress $X$ on $Y_N$, generating $\hat{Y}$. I obtain my residuals, $Y_{res}$, by subtracting $\hat{Y}$ from $Y_N$. Standard linear regression. Since $Y_N$ is 10x6 we have 6 response vectors, therefore $Y_{res}$ also has 6 vectors each containing the residuals corresponding with each response in $Y_N$. I find that the correlation among the columns of $Y_{res}$ is greater than the correlations among the columns of $Y_N$. This is true for >95% of the random experiments performed. Indeed, it looks like the regression is introducing artificial correlation, but why? Is there any way to overcome this limitation? Extending this experiment to multivariate regression, (e.g. X 10x2), further exacerbates the issue. To be more specific, I compute the correlation among the columns of $Y_N$ ($Y_{res}$) as follows: Average Square Correlations:               C = covariance(YN)    % This is a square symmetric covariance matrix.          D = C.^2 % Square each cell so that we only have positive numbers          SC(i,j) = D(i,j)/var(i)*var(j)           % i,j refer to respective cells of C, i,j selected so that only the upper triangle of matrix D (omitting diagonal) is included in computation.       (Then take the average of all `SC(i,j)`'s to obtain average squared correlation.)


title: The distribution of an estimator

text: I was reading the Wikipedia article on Robust Statistics: http://en.wikipedia.org/wiki/Robust_statistics#Properties_of_M-estimators What is meant by the distribution of an estimator, as used here, > It can be shown that M-estimators are asymptotically normally distributed, > so that as long as their standard errors can be computed, an approximate > approach to inference is available. Since M-estimators are normal only > asymptotically, for small sample sizes it might be appropriate to use an > alternative approach to inference, such as the bootstrap.


title: Comparing clusterings: Rand Index vs Variation of Information

text: I was wondering if anybody had any insight or intuition behind the difference between the Variation of Information and the Rand Index for comparing clusterings. I have read the paper "Comparing Clusterings - An Information Based Distance" by Marina Melia (Journal of Multivariate Analysis, 2007), but, other than noticing the difference in the definitions, I don't understand what is it that the variation of information captures that the rand index does not capture.


title: Categorical variables in multinomial logistic regression end up converted into binary variables

text: When I run multinomial logistic regression with some of the explanatory variables as categorical, my algo (glm) turns them in binary variables, automatically. For examples if one categorical variable X has three values a, b anc c, then my output shows cofficient and t-values for x.a, x.b and x.c. But in fact I want t-value at the level of x itself so that I can see if variable X is significant or not in determination of dependent variable. Can you please suggest some way so that I can see output directly at the level of X?


title: Comprehending output from mediation analysis in R

text: I'm trying to get my head around the mediation package in R, using the vignette for the package. I'm struggling to understand the output of the `mediate()` function.               require("mediation")     require("sandwich")     data("framing")     med.fit <- lm(emo ~ treat + age + educ + gender + income, data = framing)     out.fit <- glm(cong_mesg ~ emo + treat + age + educ + gender + income,                     data = framing, family = binomial("probit"))     summary(out.fit)     # OR for sending a message to a senator for treated compared to untreated.      exp(as.numeric(coef(out.fit)[3]))           # mediation     med.out <- mediate(med.fit, out.fit, treat = "treat", mediator = "emo",                        robustSE = TRUE, sims = 100)     summary(med.out)     ...                                   Estimate 95% CI Lower 95% CI Upper p-value     ACME (control)             0.0802       0.0335       0.1300    0.00     ACME (treated)             0.0808       0.0365       0.1311    0.00     ADE (control)              0.0142      -0.1030       0.1325    0.78     ADE (treated)              0.0147      -0.1137       0.1403    0.78     Total Effect               0.0949      -0.0316       0.2129    0.14     Prop. Mediated (control)   0.7621      -2.0926       4.9490    0.14     Prop. Mediated (treated)   0.7842      -1.9272       4.6916    0.14     ACME (average)             0.0805       0.0350       0.1304    0.00     ADE (average)              0.0145      -0.1087       0.1364    0.78     Prop. Mediated (average)   0.7731      -2.0099       4.8203    0.14     ...      Does this mean that 8.08 % of the effect of treatment is mediated through emotional state, among those that are treated? Or is this a change in the coefficient of `treat`? If someone could explain the output it would be much appreciated.


title: Unstable models, repeated crossvalidation, feature selection

text: I'm still trying to classify few (about 200) samples in a high dimensional feature space (dim=19) into 3 (very unbalanced) classes. I use an implementation of Least Squares SVM with one vs one coding with class weighting. I want to:   1. estimate the best model parameters   2. estimate the overall generalization performance My procedure is as followed: lets call this the inner loop:   1. Estimate the hyperplane-parameters via a grid search one each binary model (e.g. 3 models for 3 classes). This is done on the given data via a 5 fold crossvalidation. To be sure I'm doing it right: Each model is "optimized" separatly. The cost-criterion is the balanced accuracy (see http://en.wikipedia.org/wiki/Accuracy_and_precision#In_binary_classification) to deal with the unbalanced classes. I perform a 5 fold cv on each parameter from the grid search. So I have to make e.g. N 5-fold cv, where N is the "gridsize" for the gridsearch. I'm doing this because using no cv (though it would be computational much better) might lead to an biased estimator, right?   2. Now assume I have estimated the "right" hyperplane parameters for each model. This still does not say anything about the accuracy for the 3 classes case. But I think, "local" optimization is still better than no parameter selection? Anyway, to get a good accuracy I'm doing a second optimization step with the complete model (3 submodels for each class) using a simplex algorithm. Again I use 5-fold-cv to estimate the balanced accuracy.   3. after the last step I have a model that has "optimal" parameters. end inner loop begin outer loop   1. now to estimate the overall performance I do another crossvalidation in an "outer loop". So I split my whole dataset into 5 folds, pass only the training set to the inner loop, where the "inner cvs" and parameter selection is made. To get the generalization performance I use the test set to calculate the balanced accuracy for each fold. Each fold I save the confusionmatrix. At the end I simply add all the confusionmatrices to get the overall performance. end outer loop Now my problems are:   1. Is the procedure correct? It is already computationally heavily.    2. My model's parameters and the performance in each outer cv loop is quite different. First I did not use an outer cv loop, only a hold out partition. I experienced very different performances each time I run the script. So I assume my models are unstable. Should I therefore repeat the outer loop? Say I've made 10 iterations/repetitions of the outer loop, what informations do I get?   3. I additionally would like to include a feature selection in the inner loop. I could use `sequentialfs` from matlab. But I'm not sure where in the whole procedure the feature selection takes place. Since `sequentialfs` uses internal 10 fold cv, there would be another cv loop. That irritates me. So what is the correct order to perform all the steps? Any advice is appreciated, thank you very much and sorry for my englisch.


title: Correlation or regression? Or MANOVA?

text: I have made a summative scale based on a likert scale. That have created 7 multi item scales all based on theory and validated by a reliability analysis using cronbach's alpha showing all values above 0,7 for all the 7 multi item scales, which consist of 5-6 questions that have been added and divided by the number of questions to get the average. I am interested in testing these multi-item scales (if scale A has an effect on scale Z) and have done so by doing 12 correlations, linear regression and multiple regression. Is that correct?


title: From joint to marginal likelihood

text: I am trying to calculate the marginal maximum likelihood estimate of a parameter ρ given a vector of data _x_ for person _i_ across _n_ items. This is in the context of Item Response Theory. The first (usual IRT) parameter, θ, is the ability (theoretically unbounded but usually between -3 and 3). The second parameter, ρ, is the probability that the respondent answers in a non- optimal way. I want to integrate θ out of my joint likelihood equation: ![likelihood](http://lehrfeld.me/images/lik.jpg) where the probability function in the integrand is ![probability](http://lehrfeld.me/images/prob.jpg) The second term in the integrand, _f_ (θ), will just be a numerical approximation to the probability density function of θ. The three dots in the first set of curly brackets above can take three different forms depending on certain conditions: 0 (in the "normal" condition, so person _i_ will definitely not respond in a non-optimal way), the formula in the second set of curly brackets, or one minus the formula in the second set of curly brackets. My problem is just in wrapping my head around how I can accomplish these calculations. I am using R, but I suppose any good explanation of how I can do this computationally would be sufficient for me to understand enough to try to implement it. It is my understanding that in order to be able to carry out these calculations, I need to express the first term in the integrand (i.e., the joint probability function) as a function of 0. Then, I multiply this function of 0 by many numerical approximations to the second term (pdf of 0). Finally, I integrate each of these products over 0. Does this mean that the integrand is entirely a function of 0, and once I integrate 0 out I essentially go from a formula to a number? At which point I would sum of the results over all _n_ items, and then implement an optimization routine to search for the most likely value of ρ within the bounds allowed?


title: How do I detect state change in multivariate time series?

text: I have a multivariate time series . For each row in the data we have the values of inputs and a label for stability (0 or 1 ) . What are the algorithms that can detect the stability for an unlabelled time series, using this historical data. The sample data looks like this ![enter image description here](http://i.stack.imgur.com/lSiwX.png) Is this a machine learning problem or a time series classification problem.What are the methods that can be used and how?


title: How to fit two or more datasets with different occurence for regression

text: I want to run a regression in R with different datasets. The question is whether stock performance (daily log return) is influenced by factors like interest rates (the one set by fed or ECB), size of the board of directors (e.g. has been long 7 and then switches to 8) and whether there is an audit committee or not (binary, `1` or `0`). I guess in theory I could just add the interest rate, size of board, etc. to the daily log returns – meaning the interest rate stays fixed for 80 days and then suddenly switches to a new value. However, because this is a large dataset of 10y daily log returns, and I just cannot figure out how to do it automatically in Excel, I am wondering whether it is possible to do that in R or to run a regression on the different data sets. Note: every dataset has a corresponding date. Since 2 days ago, I've been searching the web but cannot find anything which solves my problem, so your help is really appreciated. If you have any links to papers or the like, it would be totally sufficient :) Although I did find a few papers on that topic, they were cutting-edge theory, which is way too complicated for me. Thank you!


title: How to draw a Hasse diagram of a large number of pairwise comparisons?

text: I have 100 pairwise comparisons where my results show that 33% of the pairwise comparison tests are non-significant.   * How can I draw draw a Hasse diagram to represent this?   * Are there any tips to follow when drawing a Hasse diagram?


title: Predict feature combination with highest probability

text: I trained a Support Vector Machine with the **caret** package in **R**. My **dataset** looks the following:               +----+----+----+----+----+----+----+---+     | x1 | x2 | x3 | x4 | x5 | x6 | x7 | y |     +----+----+----+----+----+----+----+---+      I have **7 X values** and one output class **Y** , which is a 2 level factor variable with "yes" and "no". So an example of this data would look like:               +----+----+----+----+----+----+----+-----+     | x1 | x2 | x3 | x4 | x5 | x6 | x7 | y   |     +----+----+----+----+----+----+----+-----+     | 3  | 3  | 2  | 3  | 6  | 7  | 2  | yes |     +----+----+----+----+----+----+----+-----+     | 6  | 3  | 2  | 6  | 2  | 6  | 7  | no  |     +----+----+----+----+----+----+----+-----+      **My Question:** What method can I use to predict the best combination of values for a new sample? So I want two find out which combination of X3 to X7 would be the best if x1 and x2 stay the same. I could do it manually and create a dataframe with every single possible combination and apply the predict function to every row and choose the one with the highest probability for "yes". An example would be:               +----+----+----+----+----+----+----+     | x1 | x2 | x3 | x4 | x5 | x6 | x7 |     +----+----+----+----+----+----+----+     | 3  | 4  | 1  | 1  | 1  | 1  | 1  |     +----+----+----+----+----+----+----+     | 3  | 4  | 2  | 1  | 1  | 1  | 1  |     +----+----+----+----+----+----+----+     | 3  | 4  | 3  | 1  | 1  | 1  | 1  |     +----+----+----+----+----+----+----+     | 3  | 4  | 4  | 1  | 1  | 1  | 1  |     +----+----+----+----+----+----+----+     | 3  | 4  | 5  | 1  | 1  | 1  | 1  |     +----+----+----+----+----+----+----+      Here I leave x1 and x2 always the same, and just change x3 to x7 to create all the possible combinations. But this is very slow. What method can I use to find the best combination? EDIT: I can use the predict function to also return the probability for each class. This would look like:               +---+--------------+--------------+     |   | no           | yes          |     +---+--------------+--------------+     | 1 | 0.6819019473 | 0.3180980527 |     +---+--------------+--------------+     | 2 | 0.6098350889 | 0.3901649111 |     +---+--------------+--------------+     | 3 | 0.4135923061 | 0.5864076939 |     +---+--------------+--------------+     | 4 | 0.6513802695 | 0.3486197305 |     +---+--------------+--------------+      Regards


title: Feature selection for SVM and Maximum Entropy

text: In text classification problems where the number of features >> number of documents, is it useful to perform feature selection with filters (e.g. Information Gain) when using Naive Bayes. However, what about SVM and Maximum Entropy classifiers ?


title: Finding the uncertainty in a Binomial probability estimate

text: In order to make some predictions for my work I have modeled a process using a binomial distribution, but in my case every single experiment must be a success and I am just changing the probability that each individual trial is a success. Effectively what I am looking at is: $$ P({\rm success}) = p^m $$ Where $m$ is the number of trials I need to succeed at. If $m$ and $p$ are known with certainty is there any way to estimate the variance of $P({\rm success})$, or am I fundamentally misunderstanding something?


title: kalman filter initialization parameter

text: Hello I don't have any idea of how to start implementing Kalman Filter in python! I have a DataFrame ( table) with in one column my forecast values and in another column my actual datas (real). The data are the same active power By using pykalman, I want to improve my forecast results and I don't know how to get this parameters: .transition_matrix, .observation_matrix, .initial_transition_covariance, .initial_observation_covariance, .transition_offsets, .observation_offset, .initial_state_mean, .initial_state_covariance


title: First Differences and Fixed Effects

text: I'm aware of the fact that first differences and fixed effects are both designed for the same solution -- removing unobserved unit-level effects. However, I'm unclear on what happens when you include a unit-level dummy in a first differences model (I've seen this done for error-correction models as well as elsewhere). Does including the dummy essentially demean the differences (sort of a soft unit-level time trend?). Or does it simply introduce bias? Thanks


title: What is the difference between a "learner" and "classifier" in supervised learning?

text: This question stems from Pedro Domingos' excellent paper "A Few Useful Things to Know About Machine Learning." The paper is extremely clear and well- written, but I still have a clarification question. Namely, what is the difference between what Domingos describes on page 1 as a "learner" and a "classifier"? I have largely taken these to by synonymous. The sentence that threw me for a loop in the paper is (again, on page 1): "The test of the learner is whether this classifier produces the correct output _yt_ for future examples **xt**." I thought this was simply the test of the classifier. Any clarification or further reading suggestions would be greatly appreciated.


title: Understanding repeated covariance types in SPSS?

text: I am working in SPSS on a repeated measures linear mixed model and I am having a really hard time wrapping my head around how to select a "repeated covariance type". The options are:   * Ante-Dependence:First Order   * AR1   * AR1 Heterogenous   * ARMA   * Compound Symmetry   * Compound Symmetry: Correlation Metric   * Compound Symmetry: Heterogenous   * Diagonal   * Factor Analytic: First Order   * Factor Analytic: First Order, Heterogeneous   * Huynh-Feldt   * Scaled Identity   * Toeplitz   * Toeplitz: Heterogenous   * Unstructured   * Unstructured:Correlations For my study, I am looking at bird calls and how they change under certain manipulated conditions.   1. Can anyone recommend any resources that clearly explain the different covariance types? I am struggling to understand their differences and in what instances you should use each. I've looked at the SPSS manual and that hasn't been helpful.   2. With so many options it seems there might be more than one "right answer" - how do you determine which one to use then? Obviously once I have better resources I will be able to eliminate some options and narrow it down, but what parameter should I be looking at to assess if I have selected the right covariance type or not? Is there a test I can do along with the LMM that will give me a value related to the covariance that I can look at to evaluate its appropriateness?


title: Intuitive explanation of Bayesian logistic regression?

text: I'm looking for an intuitive explanation of Bayesian Logistic Regression (I'm using it for texts if that's relevant). It seems that this article presents it, but it's, uh, _way_ too mathy. Thanks!


title: Correlation of Ratios substantially lower compared to the numerator alone

text: I have four variables. The first two (V1, V2) are correlated .98. The third and fourth variables are correlated (.99). When I correlate the ratios V1/V3 and V2/V4 the correlation drops to .59. Any thoughts why? V1-V3 and V2-V4 are correlated about .98. Coudl that have anything to do with it? And how? Thank you, G


title: Chi Square test to measure degree of association

text: I have a categorical variable with binary outcomes 0,1 (say Y), I want to study the association with other variables X1,X2,X3,X4,X5..... most of Xi's are categorical and few of them are continuous. I used Chisquare test to study the same, which tell wheather I can consider two variables independent or not. Further to rank the Xi's in order of association with Y . Now my question is can I use the pvalue from the ChiSquare test or standardized residual of the chi Square Test to measure the degree of association between the variables ? As far as I could understand chi sqaure test can be used just to tell whether the variable are independent or not but not the degree of association. Is there any other test to do so or chi square can do that?


title: How to perform statistical test for enrichment of gene lists?

text: In a typical ChIP-Seq experiment, we found a transcription factor `ABC` peaks on 590 genes. 242 genes out of the 590 genes are classified as `ORF`s. If the number of `ORF`s contained in the genome is 5420, and the total number of genes in the genome is 6226, are the total `ABC` bound genes enriched in `ORF`s? Can we test whether `ABC` association with ORFs is by chance? Which statistical test should be used? Can it be done in R, and if so, how? Is there any other way to do this? I got to know from this forum that tests for enrichment of gene lists include a hypergeometric test or a one-sided Fisher's exact test. Though I am not very familiar with R, based on other examples, I tried to use R for Fisher's exact test for count data, and my output is this:               fisher.test(matrix(c(242,5178,348,458),nrow=2,ncol=2),alternative="greater")           Fisher's Exact Test for Count Data          data:  matrix(c(242, 5178, 348, 458), nrow = 2, ncol = 2)          p-value = 1          alternative hypothesis: true odds ratio is greater than 1          95 percent confidence interval:      0.05224986        Inf          sample estimates:          odds ratio     0.06156519      My analysis, if correct, suggests that factor `ABC` is related to `ORF`s only by chance. Is this conclusion right? If not, what should be the conclusion?


title: High p-values for logistic regression variable that perfectly separates?

text: I'm using R to run some logistic regression. My variables were continuous, but I used cut to bucket the data. Some particular buckets for these variables always result in dependent variable being equal to 1. As expcted, the coefficient estimate for this bucket is very high, but the p-value is also high. There are about ~90 observations in either these buckets, and around 800 total observations, so I don't think it's a problem of sample size. Also, this variable should not be related to other variables, which would naturally reduce their p-values. Are there any other plausible explanations for the high p-value? Example:               myData <- read.csv("application.csv", header = TRUE)     myData$FICO <- cut(myData$FICO, c(0, 660, 680, 700, 720, 740, 780, Inf), right = FALSE)     myData$CLTV <- cut(myData$CLTV, c(0, 70, 80, 90, 95, 100, 125, Inf), right = FALSE)     fit <- glm(Denied ~ CLTV + FICO, data = myData, family=binomial())      Results are something like this:               Deviance Residuals:           Min        1Q    Median        3Q       Max       -1.53831  -0.77944  -0.62487   0.00027   2.09771            Coefficients:                     Estimate Std. Error z value Pr(>|z|)         (Intercept)     -1.33630    0.23250  -5.747 9.06e-09 ***     CLTV(70,80]     -0.54961    0.34864  -1.576 0.114930         CLTV(80,90]     -0.51413    0.31230  -1.646 0.099715 .       CLTV(90,95]     -0.74648    0.37221  -2.006 0.044904 *       CLTV(95,100]     0.38370    0.37709   1.018 0.308906         CLTV(100,125]   -0.01554    0.25187  -0.062 0.950792         CLTV(125,Inf]   18.49557  443.55550   0.042 0.966739         FICO[0,660)     19.64884 3956.18034   0.005 0.996037         FICO[660,680)    1.77008    0.47653   3.715 0.000204 ***     FICO[680,700)    0.98575    0.30859   3.194 0.001402 **      FICO[700,720)    1.31767    0.27166   4.850 1.23e-06 ***     FICO[720,740)    0.62720    0.29819   2.103 0.035434 *       FICO[740,780)    0.31605    0.23369   1.352 0.176236         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1037.43  on 810  degrees of freedom     Residual deviance:  803.88  on 798  degrees of freedom     AIC: 829.88          Number of Fisher Scoring iterations: 16      FICO in the range [0, 660) and CLTV in the range (125, Inf] indeed always results in Denial = 1, so their coefficients are very large, but why are they also "insignificant"?


title: Statistical significance of coherence values

text: I have two signals sampled at 100 Hz during 3 seconds. So I have 300 sample points at each signal series. I want to calculate coherence values between these two time series and would like to know whether this coherence value is statistically significant. I found two resources online:   1. http://journals.ametsoc.org/doi/pdf/10.1175/1520-0469%281979%29036%3C2020%3ACSL%3E2.0.CO%3B2   2. http://pmc.ucsc.edu/~dmk/notes/cohere_signif/ In the first resource, there is a parameter _n_ (it is the sample size as I could interpret). A significant coherence value for that _n_ is: $$c = \sqrt{1-\alpha^{1/(n-1)}}$$ and this formula does not include any frequency component, unless _n_ is related to frequency instead of sample size. In the second resource, the function returns degree of freedom and significance level at each frequency. However I am not sure whether this function returns the correct computation or not. Besides I couldn't interpret the results of this function. **My question:** how can I find the significance of my computed coherence value at a certain frequency (or frequency band) preferably using MATLAB in the light of the above resources?


title: Logistic regression and ordinal independent variables

text: I have found this post: > Yes. The coefficient reflects the change in log odds for each increment of > change in the ordinal predictor. This (very common) model specification > assumes the the predictor has a linear impact across its increments. To test > the assumption, you can compare a model in which you use the ordinal > variable as a single predictor to one in which you discretize the responses > and treat them as multiple predictors (as you would if the variable were > nominal); if the latter model doesn't result in a significantly better fit, > then treating each increment as having a linear effect is reasonable. > > – @dmk38 Dec 12 '10 at 5:21 Could you please tell me where can find something published that supports this claim? I am working with data and I would like to use ordinal independent variables in logistic regression.


title: Choosing a test based on by data and how to test when I have several groups?

text: I'm trying to analyze the effect of the laying order on a certain concentration of immune factor. I have a large set of data (1000+ eggs). I want to know if egg 1 (first-laid egg) will have higher concentrations than the last-laid egg (egg 7, 8, 9, etc) or even if the concentrations tend to increase/decrease with laying order. I have no idea where to start and I think I'm confusing myself because I'm looking at several nests to see if the effect happens in most cases. If anyone could help me out, or even give me a push in the right direction, I would really appreciate it. Thank you! Some of my data looks like this: ![enter image description here](http://i.stack.imgur.com/E5HwC.jpg)


title: Does GARCH(p,0) make sense at all?

text: I wonder if a GARCH model with only "autoregressive" terms and no lagged innovations makes sense. I have never seen examples of GARCH(p,0) in the literature. Should the model be discarded altogether?      E.g. GARCH(1,0):   $$ \sigma^2_t = \omega + \delta \sigma^2_{t-1} $$   From the above expression one can derive (by repeated substitution) that   $$ \sigma^2_t \rightarrow \frac{ \omega }{1-\delta } $$ for all t, if an infinite past of the process is assumed. In other words, GARCH(1,0) implies homoskedasticity and thus the "autoregressive" term, and indeed the whole model, becomes redundant.      Of course, when estimating models in practice, we do not have infinite past; but for long enough time series this approximation should be reasonably representative. Is this right? Should we never use GARCH(p,0)?


title: how to handle the imbalanced data in regression analysis

text: The problem here is very similar to the problem asked by someben in 2012 (link:Sampling for Imbalanced Data in Regression). It involves the linear regression analysis using an unbalanced dataset. Say, you have known that the variables (e.g., X and Y) are somewhat linearly related, and some data for X and Y are available. However, 1) the numbers of X and Y are not equal; and 2) X and Y are disorded, i.e., we don't know the exact matching (i.e., X1 matches Y1, Y2, or else? we just don't know). In this context, how to do the regression analysis? (Note that: if we have known the exact matching between X and Y , it is no longer a big problem because we can solve it using EM algorithm).


title: Cluster Matching/comparison

text: **Description** I have several datasets(from different subjects) with the same type of data. For each dataset I cluster the data using affinity propagation. Clustering is based on similarity distance between data points and I'm using a type of histogram intersection for the distance (dataset are histograms). For each dataset I obtain a certain number of clusters. Let's assume that stability of the cluster solution is not a problem. I want to compare/match lets say cluster A from dataset 1 with clusterD from dataset 2 and be able to determine if statistically different. **My approach** What I have tried with some mixed results is based on a Kolmogorov-Smirnov test to compare two distributions. If I am comparing ClusterA with ClusterD for example, I build distribution 1 by calculating the distance of every point of clusterA to its exemplar, and distribution2 as the distance of each point of ClusterD to the exemplar of clusterA. And then run the KS test on these distributions. **Question** I don't know if this makes any sense statistically. I have been trying to find a formal prof of this approach or something similar. But I'm willing to try any other method more suitable to the problem. If any of you could point me to the right way or offer some suggestions I would be very grateful.   Thanks so much for your time.


title: Should I run Hierarchical regression to test moderation? Why?

text: I am an undergraduate student who currently preparing for my thesis paper. In my design, there is total 4 variable (1 predictor, 1 moderator and 1 outcome variable, all of them is in continuous scale. I also have 1 control variable which is in categorical scale). I am examining the relationship between the Better-than-average bias and performance appraisal system satisfaction while controlling for gender, and testing the Equity Sensitivity as the moderator in this relationship. Is it appropriate for me to run Hierarchical regression? If yes, may I know why? Is there any references that i can use to cite the reasons?


Sequential pattern mining finds patterns in unlabeled data. So I assume its an unsupervised learning algorithm. But most references to pattern mining never mention the learning part. Why? Can sequential pattern mining be called an unsupervised learning method? This Wikipedia articles lists the Apriori algorithm as an unsupervised method, but no citations.


title: How can I calculate the probability of model $g_i$ given a set of $n$ models and AIC values?

text: ## Background I have a data set $Y$:               set.seed(0)     predictor <- c(rep(5,10), rep(10,10), rep(15,10), rep(20,10)) + rnorm(40)     response  <- c(rnorm(10,1), rnorm(10,4), rnorm(10,2), rnorm(10,1))     plot(predictor, response)      and a set of models g$_i$:               fits <- list()     fits[['null']]      <- lm(response ~ 1)     fits[['linear']]    <- lm(response ~ poly(predictor, 1))     fits[['quadratic']] <- lm(response ~ poly(predictor, 2))     fits[['cubic']]     <- lm(response ~ poly(predictor, 3))     fits[['Monod']]     <- nls(response ~ a*predictor/(b+predictor),                                start = list(a=1, b=1))     fits[['log']]       <- lm(response ~ poly(log(predictor + 1), 1))      ## Problem I can find the best fit using               library(plyr)         ldply(fits, AIC)      I am revising a manuscript with results originally presented as raw AIC values, but I find the actual AIC values are fairly uninformative and difficult to interpret ## Question Can I calculate the probability of each model given the set of models e.g. g$_i$ $$\frac{P(Y|\textrm{g}_i)}{\sum_{j=1:n}P(Y|\textrm{g}_j)}$$ An approach that most similar to using AIC wins, because I would like to minimize the amount of revisions to the methods and results that will be required. ## Other considerations:   1. I have considered using ANOVA although this does not appear to have as straightforward of an interpretation as the above, or does it?:              anova(fits[[1]],fits[[2]],fits[[3]],fits[[4]],fits[[5]],fits[[6]])        2. A Bayesian approach would require more work on both implementing statistics and because it would fundamentally change the interpretation from the likelihood interpretation of $P(Y|g)$ to the Bayesian $P(g|Y)$.


**What exactly is a higher-order correlation?** Is the term just a catch-all for all dependencies that are not captured by second-order correlations? Or do people mean something more precise than that? Can all higher-order correlations be captured by second-order correlations on some nonlinear function of the data? Is a higher-order correlation necessarily a function of more than two variables? Or can one talk about a "higher-order correlation" in relation to simply two variables? Thanks in advance!


title: Distinction between linear and nonlinear model

text: I have read some explanations about the properties of linear vs nonlinear models, but still I am sometimes not sure if a model on hand is a linear or a nonlinear one. For example, is the following model linear or nonlinear? $$y_t=\beta_0 + \beta_1B(L;\theta)X_t+\epsilon_t$$ With: $$B(L;\theta)=\sum_{k=1}^{K}b(k;\theta)L^k$$ $$L^kX_t=X_{t-k}$$ Where $b(k;\theta)$ represents (a decaying) Exponential Almon Polynomial function of the form: $$b(k;\theta)=\frac{exp(\theta_1 k+\theta_2k^2)}{\sum_{k=1}^{K}exp(\theta_1k+\theta_2k^2)}$$ In my view, my main equation (the first one) is linear with respect to $X_t$, because this term is just multiplied with a weight. But I would say the weighting function (the last equation) is nonlinear with respect to the parameters $\theta_1$ ans $\theta_2$. Can someone explain to me if my main function is a linear or a nonlinear one and what does it mean for the estimation procedure - do I have to apply linear or nonlinear least s squares method?. Furthermore, what is the discernible feature by means of which I can definitely identify if a function is a nonlinear or linear one? Thanks!


title: Obtaining prediction equation from principal components analysis using 1st orthogonal component

text: Suppose you've collected data on factors $X_1, X_2, X_3, \ldots X_p$ and you're interested in creating a latent factor score $L = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \ldots \beta_p X_p$ . For instance, these variables might have been part of a regimen of pain testing and you wish to "score" individuals for their overall pain sensitivity. How does one extract the $\beta$s as from the question above and calculate their standard errors? For instance, in R, I can run:               x <- matrix(rnorm(1000), 100, 10)     p <- princomp(x)      What's the linear equation that gives me the first orthonormal component of the $x$ matrix?


title: How does sampling scheme affect some statistics of interest?

text: I have a sample of data sets A = {a1, a2, a3, ...aN} We will change each data value following the same pattern, say a1’ =a1*10, a2’ = a2*10. Which statistics evaluated from this data set will be changed? Assume that I increase the size of this data set simply by make copies for each data, A’ = {a1, a1, a1, a2, a2, a2, a3, a3, a3,...} How will this changing affect the different statistics, like mean, standard deviation, etc.?


title: How to use cross validation to estimate autocorrelation?

text: I have a stationary time series and I want to calculate the autocorrelation coefficient of order 1. For that I use OLS. I know the autocorrelation parameters change over the time. Thus it is not trivial to fix a window length for the regression. Moreover I can easily do over-fitting by finding the windows which optimize a criterion such as mean squared error, for example. Can I use cross validation to find the optimal window length?


title: How do I handle very different weights in a least squares fit?

text: I'm performing a weighted linear least squares fit, where the weights correspond to the number of counts of a specific observation. Due to the nature of the data, it is possible that a small handful of observations get weighted much more than the rest, so that the regression is dominated by, say, two data points only, which obviously leads to bad results. I've run a few tests and I noticed that by e.g. capping the weights, I can improve the results, though this feels rather hack-ish to me. Are there better ways to avoid a small number of data points to out-weigh the rest of the data, or, if capping the weights is a good approach, is there non- _ad hoc_ way to determine the optimal value of the cap? **EDIT** Here's what some sample data and fits look like with (a) no weights, and (b) weights. The data were generated from fairly realistic simulations, so I know the ground truth (red line). ![enter image description here](http://i.stack.imgur.com/IKCE5.png) ![enter image description here](http://i.stack.imgur.com/kaHPa.png) My problem is that the weight of some data points (#2 at x=10 for example) can be sufficiently large to dominate the fit. However, I also don't want the very-low-count data to weigh in too much, otherwise I get a really crappy fit as well.


title: How to visualize an interaction effect from a regression?

text: I have a regression with a harmonic effect of day of the year, which interacts with other variables. I am not sure how to interpret the coefficients. My model is:               m1 <- lme(lcount ~ AirT + sin(2*pi/360*DOY) + cos(2*pi/360*DOY) +                AirT*sin(2*pi/360*DOY) + AirT*cos(2*pi/360*DOY) + RainAmt + RainAmt*AirT,               random = ~1|plot))      I get significant interaction effects of air temperature with the linearized harmonic day of the year (DOY) function. My response variable is the log of animal counts on each day. I want to describe how the effect of air temperature on animal observations changes depending on the day of the year. Does anyone have suggestions on how I can interpret my beta values and/or how I can visualize the effect? I am using R but am not that skilled. The package I used for analyzing my data is nlme. EDIT: **My primary goals are (1) to describe the response of animals to environmental variables and (2) to predict future activity periods (i.e. when and under what conditions should a research bother trying to catch these animals).** So if there is a better way to model this data, I would be interested in hearing it (such as cubic splines - see comments below).


title: Should I use multi-label classification?

text: I have a classification problem with 2 classes (positive and negative). Usually, in such classification problems, all the samples will be labelled either 'positive' or 'negative'. In my dataset, some of the samples possess a combination of both positive and negative characteristics. Formally, if the dataset is $x$, then, > $x=x_1 \cup x_2 \cup x_3$ where $x_1$ is the set of all positive samples, $x_2$ is the set of all negative samples and $x_3$ is the set of samples that contain the characteristics of both the classes. As far as I could think of, this situation could be handled in 2 ways,   1. Ignore $x_3$ (samples that contain characteristics of both classes) and treat the problem as a traditional binary classification.   2. Label the samples in $x_3$ with both labels (positive & negative) and consider this as a multi-label classification problem. I wish to follow the second option, as it is more natural choice. The reason being, ignoring some samples from the dataset gives me a feeling of manipulating the dataset artificially which may affect the performance of the classifier in the real world scenario. In this context, I have the following questions,   1. Is it correct to treat this a multi-label classification problem. If so, is the intuition which is explained in the previous paragraph correct?   2. Is there any other learning paradigm that can handle this scenario? If so, please provide reference to relevant literature.


title: Bayesian networks from a table

text: Could someone help me with question 5.b. I understand that the probability of any of these occuring independently is 0.5 but how do I combine those into a joint distribution function? Is $0.5 \cdot 0.5 \cdot 0.5 \cdot 0.5 = 0.0625$ correct? ![enter image description here](http://i.stack.imgur.com/P1hud.png) Edit: Is it more like this? $S =$ Spam, $D =$ Dear, $F =$ Free, $H =$ Hot $$ P(S, D, F, H) = P(S) \cdot P(D | S) \cdot P (F | S, D) \cdot P(H | S, D, F) $$ Therefore: ![enter image description here](http://i.stack.imgur.com/44ujU.png) or am I missing something?


title: Local log-likelihood for multiclass linear regression model

text: In page 206 of the book 'Elements of statistical learning', the author wrote: > _The local log-likelihood for this $J$ class model can be written_ > > $\sum_{i=1}^NK_\lambda (x_0, x_i)\{\beta_{g_i0}(x_0) + > \beta_{g_i}(x_0)^T(x_i-x_0) - > log[1+\sum_{k=1}^{J-1}exp(\beta_{k0}+\beta_k(x_0)^T(x_i-x_0))]\}$ I can understand that the term $K_\lambda (x_0, x_i)$ is there to weight down the log-likelihood of each individual observation, but really don't know why the term $(x_i-x_0)$ instead of only $x_i$ - which I am expecting. Is this just a typo or am I misunderstanding something?


title: Why do we not look at the covariance matrix when choosing between LDA or QDA

text: I understand the difference between LDA and QDA (linear and quadratic discriminant analysis), being that with LDA assume that your features have the same covariance matrix in each class. I wonder why I have not seen an example yet where they actually calculate the variance covariance matrices in each class and compare them with each other? Or maybe you can use the pairs plot for this (ref below)? Then how do I read it? Or should't I care too much, and just choose the model based on the prediction error? ![enter image description here](http://i.stack.imgur.com/QPvwu.jpg) Thank you!


title: Kalman Filtering, Smoothing and Parameter Estimation for State Space Models in R and C#

text: I am in the process of writing an open source State Space Analysis suite in C# (for fun). I have implemented a number of different Kalman-Based Filters (Kalman Filter, Information Filter and the Square Root Filter) and the State and Disturbance Smoothers that work with these implementations. My tests on these filters (using the Nile data from _Durban and Koopman's (DK)_ book _"Time Series Analysis by State Space Methods"_ and other more complex data) show that the filters and smoothers work and they produce very similar results (as you would expect for a local univariate model). The smoothed output for the basic Kalman Filter for the Nile data looks like (because everyone likes a graphic or two :]) ![Smoothed Nile Data](http://i.stack.imgur.com/MIZWI.png) The dotted line is the 90% confidence interval. Now, my problem; I am now attempting to write the first part of the parameter estimation code and for the parameter estimation of the covariances $\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$ in the linear Gaussian State Space model $$y_t = \mathbf{Z}_{t}\alpha_{t} + \epsilon_{t},$$ $$\alpha_{t + 1} = \mathbf{T}_{t}\alpha_{t} + \mathbf{R}_{t}\eta_{t},$$ where $y_{t}$ is our observation vector, $\alpha_{t}$ our state vector at time step $t$ and $$\epsilon_{t} \sim NID(0, \mathbf{H}_{t}),$$ $$\eta_{t} \sim NID(0, \mathbf{Q}_{t}),$$ $$\alpha_{1} \sim NID(a_{1}, \mathbf{P}_{1}).$$ where $t = 1,\ldots, n$. We have an implementation of the Expectation Maximisation (EM) algorithm and this is estimating $\mathbf{H} = 19550.37$ and $\mathbf{Q} = 7411.44$ with the loglikelihood converging in 32 iterations to -683.1 with a convergence tolerence of $10^-6$. The DK book quotes $\mathbf{H} = 15099$ and $\mathbf{Q} = 1469.1$ and from Tusell's paper and this walkthrough using R he matches this prediction getting $\mathbf{H} = 15102$ and $\mathbf{Q} = 1468$. I have debugged my code and it really seems like my implementation of the EM algorithm is correct. So get to the bottom of what is going on I want to run the same data set using R and KFAS... Using R and the R Package KFAS, I have attempted to reproduce these estimates for $\mathbf{H}$ and $\mathbf{Q}$, but my R knowledge is weak. My current R script is as follows               install.packages('KFAS')     require(KFAS)          # Example of local level model for Nile series     modelNile<-SSModel(Nile~SSMtrend(1,Q=list(matrix(NA))),H=matrix(NA))     modelNile     modelNile<-fitSSM(inits=c(log(var(Nile)),log(var(Nile))),model=modelNile,     method='BFGS',control=list(REPORT=1,trace=1))$model          # Can use different optimisation:      # should be one of “Nelder-Mead”, “BFGS”, “CG”, “L-BFGS-B”, “SANN”, “Brent”     modelNile<-SSModel(Nile~SSMtrend(1,Q=list(matrix(NA))),H=matrix(NA))     modelNile     modelNile<-fitSSM(inits=c(log(var(Nile)),log(var(Nile))),model=modelNile,     method='L-BFGS-B',control=list(REPORT=1,trace=1))$model          # Filtering and state smoothing     out<-KFS(modelNile,filtering='state',smoothing='state')     out      _How can I extend this R script using KFAS in order to estimate the parameters $\mathbf{H}$ and $\mathbf{Q}$?_ and _Should each method of parameter estimation be hitting the same values for $\mathbf{H}$ and $\mathbf{Q}$?_ Thanks for your time.


title: Explanation on One Class SVM

text: I was using One class SVM implemented in Scikit learn, Python for my research work. But I have no good understanding of this. Can anyone please give a simple, good explanation of One Class SVM? Thanks very much in advance.


title: How to test for relationship between cumulative intake and outcome over time in single arm study?

text: I have a one group trial with n = 100. I want to analyze the relationship between the accumulated amount of drug intake (continuous) and the effect (measured by symptom score). For example, for different time points for patient 1:               t1: d_amount 10, S_score 50     t2: d_amount 15, S_score 45     t3: d_amount 20, S_score 40      For example, for different time points for patient 2:               t1: d_amount  5, S_score 50     t2: d_amount 10, S_score 45     t3: d_amount 30, S_score 40      How would I construct a (formally correct) hypothesis proving the statement: "There is a significant negative relation between the amount of drug taken (over time) and the symptom score"?


title: How to analyse individual consistency from three non-normally distributed variables?

text: I would like to know if individuals were consistent across three different days. I measured the latency of individuals for a maximum of ten minutes. Many individuals had a latency of 0 and some individuals never arrived and were therefore given a latency of 601. The data on each of the three days is very non-normally distributed and I do not manage to transform the data (tried log, log10, sqrt etc). Running spearman rank correlations works on two of the three days but is not possible with three variables. Running ICC analyses is possible with three variables but does not seem valid with my non-normally distributed data. How can I analyse the consistency of individuals across the three days? Could I rank the data and then run ICC analysis? Here is the dataset:                    ID d1  d2  d3     1   198 21  5   0     2   192 0   0   0     3   199 1   601 3     4   191 601 601 601     5   197 601 601 313     6   194 0   23  0     7   200 23  263 2     8   143 18  130 2     9   196 19  52  0     10  195 8   94  6     11  194 0   14  0     12  199 0   0   0     13  197 0   19  0     14  129 26  41  1     15  198 13  489 0     16  191 0   0   0     17  170 0   601 1     18  162 0   29  0     19  164 1   93  18     20  168 0   29  0     21  167 0   0   0     22  169 0   601 3     23  166 0   592 0     24  163 23  33  96     25  154 0   601 0     26  155 0   601 0     27  159 0   601 0     28  156 0   90  0     29  158 0   94  0     30  160 0   0   0     31  157 0   60  0     32  141 0   601 0     33  179 0   0   0     34  176 0   43  1     35  161 4   157 0     36  142 1   74  0     37  172 0   2   0     38  178 2   367 3     39  147 1   4   0     40  144 1   601 38     41  143 0   11  0     42  141 0   601 601     43  146 124 601 0     44  142 0   601 0     45  148 0   0   0     46  134 10  73  12     47  180 0   288 0     48  133 0   600 0     49  165 6   577 0     50  131 11  199 57     51  174 42  601 224     52  171 0   346 0     53  140 0   127 0     54  123 2   14  0     55  121 0   10  0     56  127 0   50  0     57  122 7   42  0     58  124 0   131 0     59  125 0   15  0


title: Partial least squares regression for categorical factor in R

text: I adjust the partial least squares regression for one categorical factor (2 levels – `be` or `nottobe`) with with the `pls` package in R. I try to use `round()` function in the predict values for take the decision if the result are the first or second level in my factor. Does this approach sound correct?               require(pls)           #Artificial data            T<-as.factor(sort(rep(c("be", "nottobe"), 100)))           y1 <- c(rnorm(100,1,0.1),rnorm(100,1,0.1))      y2 <- c(rnorm(100,10,0.3),rnorm(100,10,0.6))      y3 <- c(rnorm(100,10,2.3),rnorm(100,11,2.6))      y4 <- c(rnorm(100,5,0.5),rnorm(100,7,0.5))      y5 <- c(rnorm(100,0,0.1),rnorm(100,0,0.1))           #Create the data frame      avaliacao <- as.numeric(T)      espectro <- cbind(y1,y2,y3,y4,y5)      dados <- data.frame(avaliacao = I(as.matrix(avaliacao)), bands = I(as.matrix(espectro)))           #PLS regression     taumato <- plsr(avaliacao ~ bands, ncomp = 5, validation = "LOO", data=dados)      summary(taumato)           #Components analysis      plot(taumato, plottype = "scores", comps = 1:5)                #Cross validation      taumato.cv <- crossval(taumato, segments = 10)      plot(MSEP(taumato.cv), legendpos = "topright")      summary(taumato.cv, what = "validation")      plot(taumato, xlab ="medição", ylab="predição", ncomp = 3, asp = 1, main=" ", line = TRUE)                #Predition for 3 components      T<-as.factor(sort(rep(c("be", "nottobe"), 50)))           y1 <- c(rnorm(100,1,0.1),rnorm(100,1,0.1))      y2 <- c(rnorm(100,10,0.3),rnorm(100,10,0.6))      y3 <- c(rnorm(100,10,2.3),rnorm(100,11,2.6))      y4 <- c(rnorm(100,5,0.5),rnorm(100,7,0.5))      y5 <- c(rnorm(100,0,0.1),rnorm(100,0,0.1))           espectro2 <- cbind(y1,y2,y3,y4,y5)      new.dados <- data.frame(bands = I(as.matrix(espectro2)))      round(predict(taumato, ncomp = 3, newdata = new.dados))##


title: Hypothesis testing if two random variables come from a related underlying function

text: I have observations of two random variables $S_1(t)$ and $S_2(t)$ and I'm trying to build a test to see if they are related. In other words I'm trying to test if unknown functions $F_1(t)$ and $F_2(t)$ are the same function. The parameters $r$ and $b$ are constants with $ r \geq 0$ and $n_1$ and $n_2$ are Gaussian noise. $$ S_1(t) = F_1(t) + n_1 $$ $$ S_2(t) = rF_2(t) + b + n_2 $$ It's been a while since I've taken a statistics class and I'm having issues remembering how to build a test for this. I want to set $H_0: F_1(t) \neq F_2(t)$ and $H_1 : F_1(t) = F_2(t)$ but after that I'm getting lost.


title: Two samples of the same distribution

text: What non-parametric methods exist to assess whether two random samples belong to the same distribution?