title: Correlation of variables accounting for background variation

text: I am investigating the association between fruiting weeds and fox activity. The hypothesis is that fox activity (sightings per week) is positively correlated with availability of fruit (fruits per square meter). However there are also independent landscape-scale variations to fox activity which are unpredictable (due to many factors such as rainfall). I therefore wish to test for an association between fox activity and fruit availability at weedy sites, using data from weed-free sites to discount the effects of independent fox activity variability. I also preferably wish to test this across multiple sites for each treatment (weedy and weed free) to reduce the risk of false-positive correlation results. The data would be:                          Weedy Site 1   Weedy Site 2...(etc)  Weed-free 1      Weed-free 2...(etc)              fruits| foxes     fruits| foxes        fruits| foxes     fruits| foxes      Week 1            Week 2     Week 3...(etc)      Are there any suggestions for testing correlations across replicates, and for removing background variability (trends in activity at weed-free sites) when testing correlations?


title: Radial basis network character recognition

text: I want to develop a simple character recognition program by implementing a given neural network kind; a simple command line-type is enough. The radial basis function neural network was assigned to me and I already studied the weight training, input-to-hidden-to-output procedures but I am still doubtful of in implementing it. My references are (1) and (2). A simple one-dimensional array of a 10 by 10 binary object (that represents a character) is the input. For example, the array below               input = array(          0,0,0,1,1,1,1,0,0,0,         0,0,1,0,0,0,0,1,0,0,         0,1,0,0,0,0,0,0,1,0,         1,0,0,0,0,0,0,0,0,1,         1,1,1,1,1,1,1,1,1,1,         1,0,0,0,0,0,0,0,0,1,         1,0,0,0,0,0,0,0,0,1,         1,0,0,0,0,0,0,0,0,1,         1,0,0,0,0,0,0,0,0,1,         1,0,0,0,0,0,0,0,0,1 )      is the representation of the character "A":               0   0   0   1   1   1   1   0   0   0     0   0   1   0   0   0   0   1   0   0     0   1   0   0   0   0   0   0   1   0     1   0   0   0   0   0   0   0   0   1     1   1   1   1   1   1   1   1   1   1     1   0   0   0   0   0   0   0   0   1     1   0   0   0   0   0   0   0   0   1     1   0   0   0   0   0   0   0   0   1     1   0   0   0   0   0   0   0   0   1     1   0   0   0   0   0   0   0   0   1      I plan to take the total weight of the input and compare it to the training set as in the saved 1-D arrays of the other characters of the alphabet and the one with the closest is the prediction. The problem is I tend to understand algorithms better if presented in a CLRS- manner or similar type as opposed to mathematical formula. I find it hard to understand the explanations in those two papers (which I find the easiest to read among others here in the Google search). Can someone point me to a friendly algorithm for a RBNFF that takes in an array and produces an output of weights? If not, a paper that explains this in Layman's manner would be appreciated.


title: How to characterize the distribution of the intersection of 2 bivariate normals

text: I have two 2-dimensional Gaussian distributions:   $$ D_1 := \mu_1=\pmatrix{x \cr y}, \quad \Sigma = \pmatrix{{\rm var}(x) &{\rm cov}(xy) \cr {\rm cov}(yx) &{\rm var}(y)} \\ D_2 := \mu_2=\pmatrix{x \cr y}, \quad \Sigma = \pmatrix{{\rm var}(x) &{\rm cov}(xy) \cr {\rm cov}(yx) &{\rm var}(y)} $$ and I want to characterize their 'intersection space' either analytically or with computationally generated samples. My ultimate goal is to be able to say something about the collisions of D_1 with D_2. Both are probabilistic representations of moving objects in 2-space. One sampling idea was to sample each space (using a normal random number generator) and find the joint probability for each point. This seems to require me to pick a sensible threshold, which doesn't seem like the right direction. I came across this CV thread, which yields a percentage. Can this be extended to 2 dimensions? Any help would be appreciated by way of solutions, articles or avenues.


title: Randomization test and comparisons of many to 1

text: I have about 80 different groups, each with 3 replicates which I would like to compare against a control having always a value of 1. As I have a low number of replicates and many groups having 2 replicates with 0 values, I would rather do a randomization test than pairwise mann-u test. However, I'm not sure wether I can resample across all groups to then calculate the difference with 1 (the value of the control) or if this would only tell me the probability of having at least one group that is different from one. So far, my strategy was to resample all data points 100000. Then calculate the mean of every third value in all resampled data frames. From each of this I then calculated the difference with respect to one (the control) to get a null distribution of the differences. Finally, I calculated the p-value for each group in my original data set and adjusted by FDR correction. I am afraid this is not valid, and I then will have to resample each group and the control separatel? However that meansI would have to randomise 6 values, 3 from the group and 3 times 1 from the control, for 80 groups as many times as possible?


title: How to "uncluster" a set of financial data?

text: I am attempting to evaluate and compare the profit factor of different "test runs" of a FOREX trading strategy. My problem is that, despite an average time between orders of 2hr+, some of these runs can have 20+ orders in a row, every 5 minutes, in the same direction. I need some way to normalize these clusters that occur without just throwing the data out. I want to treat the cluster as 1 data point by averaging the gain/loss of each trade within the cluster. I was thinking of doing it with a moving time window in the following manner:               For each order:     Weight = 1/(N+1)      N =  count of consecutive orders within 30 minutes of the current order.      But I am not sure if that is correct.


title: Statistical validity of before and after tests

text: I'm working at a company developing web applications and due to the load on our testing team that can execute and analyze A/B or multi-variate tests, we're being asked to run "before and after" tests--also called "launch and learn" tests. With these before and after tests, if we make a change to a web page, we're being asked to measure a key performance indicator (KPI) for a week before the change and then a week after the change. However, in many cases, the traffic to that page is lumpy and not necessarily consistent. As an example, a link on one part of a web page may be served up 1000 times in a 7 day period and clicked on 500 times during that week, but during the 2nd week (8 days, not 7)--once a change has been made and the link has been moved to another part of the page--it might be served up 1200 times and clicked on 550 times. My gut tells me that this is not an "apples to apples" comparison, and so comparing 500 clicks (A, control) to 550 click (B) would not be valid. Can someone please validate or let me know if I'm wrong?


title: The unimodality of a generalised binomial distribution

text: In my research, I'm summing a bunch of non-independent, not identically distributed random (0,1)-variables $X_1,X_2,\ldots,X_n$ with corresponding probabilities $\pi_1,\pi_2,\ldots,\pi_n$. The resulting distribution is multimodal, and I want to show that the multimodality is caused by the dependence. In an attempt to show this, I'm assuming that the probabilities are independent, and want to show that the resulting distribution is unimodal (which I'm not 100% sure is true at this point, although some computational results suggest it is true). So, the distribution is given by \[\mathrm{Pr}[\sum X_i=r] = \sum_{C \subseteq [n]:|C|=r} \big( \prod_{i \in C} \pi_i \big) \big( \prod_{i \not\in C} (1-\pi_i) \big).\] This is a generalised binomial distribution, (we retrieve the binomial distribution in the special case of $\pi_1=\pi_2=\cdots=\pi_n$) I find it hard to believe that this has not be studied previously, which leads me to: > **Question** : Does this distribution have a special name I can search for? > What are some references where I can find this distribution? I found quite a few different generalised binomial distributions via Google, but not one that seems to match the above. Or if you happen to know the answer to the underlying question: > **Question** : Is the distribution unimodal?


title: Why does a paired t-test (when appropriate) result in better variance?

text: I didn't quite get the Wikipedia explanation here: http://en.wikipedia.org/wiki/Paired_difference_test#Use_in_reducing_variance I agree that both the unpaired and paired means are the same...then I see how $\text{var}(\bar{Y_2} - \bar{Y_1})$ will include a covariance term..but what is the alternative variance? Take the paired differences first, and then the variance of the result? And how would that be any different?


title: So how many staticians *does* it take to screw in a lightbulb?

text: Sorry, but the emptyness was a bit overwhelming. And this has been stuck in my head since it got asked at Area51!


title: Non-parametric test for comparing a function of two samples

text: The first thing that comes to mind when comparing some property of two samples is probably the independent samples _t_ -test. Sometimes, people point out that this has the inherent assumption of the _t_ -distribution. Does any one know of a non-parametric alternative to the independent samples _t_ -test?


title: Can statistical prediction be asymmetric?

text: Thus can we have random variables $X_1 , ... , X_n$ , each with zero mean and unit variance, such that for a sizeable representative sample: (1) in the least-squares regression equation for $X_1$ against $X_2 ,..., X_n$ , the variable $X_n$ has a large coefficient, and its contribution is rated highly significant, while the other variables have markedly smaller coefficients and lower significance; (2) in the least-squares regression for $X_n$ against $X_1 ,..., X_{n-1}$ , the variable $X_1$ has a small coefficient and low significance? The reason for this question is that this effect appears in some real data (using SPSS), with five variables. But the result seems paradoxical; so perhaps we are doing something wrong.


title: Error bars when the variables are related

text: I am a stats novice so I'm sorry if my jargon is a bit off. I have several large samples of workers in different demographic categories from two different fields A and B. I want to compare the proportion of people in each sample (not of equal size) and try to figure out if there is a difference in the distributions the samples were drawn from. I computed the average percentage of workers in each demographic category for field A and field B. However, I wanted to construct error bars for my estimates. I realize that because I am dealing with percentages, the values in each category are all related. If one of the values is larger, the others must be small and so on. Does this affect the confidence intervals? Does it make sense to just compute the standard deviation of all "older males in field A" for instance and use that for my error bar?


title: Question on AIC and stepAIC

text: AIC(lm(Fertility ~ ., data=swiss))     [1] 326.0716      ok, since AIC is calculated as               -2*logLik(lm(Fertility ~ ., data=swiss)) + 2*7      Why does stepAIC return a smaller AIC?               stepAIC(lm(Fertility ~ ., data=swiss))     Start:  AIC=190.69


title: Maximum likelihood estimation procedures for state-space linear models

text: State-space models are represented by a state equation and an observation equation (or system of equations to be more precise). These equations are parametarized by components including a transition matrix (FF in some notation) and GG respectively. These component matrices can have large dimensions. Indeed, the log-likelihood function is a non-convex function and its maximization is difficult. Typical optimization methods used include Newton-Raphson and EM algorithm. I am new to using state-space models and am looking for quick and robust optimization procedures used by practitioners of state-space models. Is there any literature or best practices regarding the best ways to estimate these matrices?


title: Regression with repeated data

text: I have a question regarding regression analysis on a dataset were the input values generate different results over time: e.g.               1 2     2 2     3 5     4 1     2 5     3 8      How would I go about doing the regression on such a dataset, since the values change?


title: Why do I get different BIC values when I use regsubsets and lm in R

text: I used `regsubsets` to find a model with lowest BIC; height is our D.V. , the code I typed is below:               male = read.table(file.choose(), header=TRUE)     mreg = regsubsets(height ~ biacromial + pelvic.breadth + bitrochanteric + chest.depth               + chest.diam + elbow.diam + wrist.diam + knee.diam + ankle.diam               + shoulder.girth + chest.girth + waist.girth + navel.girth + hip.girth               + thigh.girth + bicep.girth + forearm.girth + knee.girth + calf.girth               + ankle.girth + wrist.girth + age + weight             , data=male)     plot(mreg)      ![enter image description here](http://i.stack.imgur.com/MjB1g.png) so the best subset for male is: `bitrochanteric,waist.girth+hip.girth+thigh.girth+bicep.girth, calf.girth+weight` I regress the model using `lm`                mreg2 = lm(height ~ bitrochanteric + waist.girth + hip.girth + thigh.girth                   + bicep.girth + calf.girth + weight, data=male)      BIC(reg2)      Then I got a value of 1461.665 ,which is totally different from my graph and so I don't understand at all why it is different.


title: Continuity correction for Pearson and McNemar's chi-square test

text: I think this is a fairly stupid question, but I just realize that I don't quite understand the term continuity correction. I use R and found the same syntax `correct=TRUE` for both `chisq.test` and `mcnemar.test`. Are they referring to different continuity correction method? I have heard of the yate's continuity for pearson chi-square test is not very popular now as it may "over-adjust" the result, then how about that for McNemar's chi-square test? Thanks.


title: Considering residuals as new dependent variables after applying a linear mixed effect model to a set of data

text: I have a dataset including for each subject :   * 22 DV (linears, range [0-1]. These DV are intercorrelated). These measures are volume of different brain regions.   * 1 main effect to test (disease status; 0,1 or 2)   * 2 confounding IV (Age of the subject (linear), sex(binary 1 or 2) ; which have a known impact on my measures)   * 1 "site effect", as both subjects and measures are from different sites (1, 2 or 3) My aim is to evaluate the effect of disease status on each of my 22 measures, after controlling for age, sex and site effects. I have planned to use MANCOVA but I have to deal with non homogenous regression slopes. To overcome this problem, I have been told that I should use a linear mixed effect model, 1 model for each DV :               Measure1 ~ AGE:fixed + SEX:fixed + DISEASE_STATUS:fixed + (SITE:random)      The problem is that I would like to test whether there is a global main effect of the disease on my 22 measures, and in a second part of the analysis check which measure is significantly affected by the disease status. **Questions:**   1. Do you think that I can run my analysis this way : first apply a model for each measure `measure ~ AGE:fixed + SEX:fixed + (SITE:random)` and then perform my analysis on the residuals of this first model (ie, a multivariate ANOVA) to check if there is a disease effect on residuals, and thus on DVs ?   2. If analysis of residuals are not appropriate, do you know if there a way to use some kind of "multivariate linear mixed effect model" ? The reason why I do not want to apply one model to explain each DV is the multiple comparison issue, as I would dramatically inflate my alpha level. Thank you in advance for your advices, and thank you for all the previously published answers that are very useful to non-statisticians researchers ! * * * **Edit** Thank you for your reply. I think I'll follow your advice and I'll use an FDR correction. My first question has been focused on a more general problem concerning the use of residuals as DV. Does the community have an opinion about this strategy?


title: How to structure data for SPSS (percentages for log-linear analysis)?

text: I am a bit desperate because I am writing my Msc thesis and I am not sure how to organize my data for an SPSS usage. My research examines the performance of mobile banner campaigns and how the 3 independent variables: Source(in-app or mobile web promotion), Style (animated or static banner) and Size (full screen or small banner) will affect the dependent variables: Click-Through Rate((total clicks/total impressions)*100) and Conversion Rate((total installs/total clicks)*100). Initially I was considering ANOVA tests like the ones below:               Click-Through Rate = β0 + β1*Source + β2*Style + β3*Size + β4*(Source*Style) +                          β5*(Source*Size) + β6*(Style*Size) + β7*(Source*Style*Size)  + error           Conversion Rate = β0 + β1*Source + β2*Style + β3*Size + β4*(Source*Style) +                        β5*(Source*Size) + β6*(Style*Size) + β7*(Source*Style*Size) + error      But my supervisor advises that I can not use percentages withing ANOVA and I should use loglinear analyses. Can you please advice me if I can somehow still use ANOVA for my tests or if I need to use loglinear analyses and how to structure my data for SPSS use. My supervisor is mentioning I need to organize it into counts, but doesn't give me more details and I am a bit confused: what does that mean?


title: Why does the $\alpha$ parameter grow to $\infty$ in Pareto-distributed random numbers when changing the threshold?

text: I noted a strange fact. Let X be a set of Pareto distributed random number with $\alpha$ and $x_{min}$ defined _a priori_. Now, let $\alpha'$ be the estimated value of the shape and $x_{min}'$ be a new threshold fixed _a posteriori_. If $x_{min}'\rightarrow max(x)$ then $\alpha'\rightarrow +\infty$. where $max(x)$ is the max value in the sample. I will try to explain me better. Look at this code Comparing Pareto fitting methods now imagine to add after the line               hh1<-(matrix(rpareto(100,alpha,0.1),ncol=1))       (where 0.1 is $x_{min}$) this piece of code               hh1<-subset(hh1,hh1 > xmin)       Here xmin=$x_{min}'$ namely the point where i start the fit or if you prefer a cut-off. Of course $x_{min}'$ cannot be greater then $max(x)$. Now imagine to put the code in a loop in order to see the behavior of $\alpha'$ while you increment the cut-off $x_{min}'$. What I would expect is when $x_{min}'\rightarrow max(x)$ then $\alpha'\rightarrow\alpha_0$ because they are independent. but this not happens and $\alpha'$ (the extimated) tend to infinity. Why do we have such a strange behaviour during the estimation? well i understand the behavior in the MLE method, but the others? The same happens with alpha-stable distributions when we try make a measure of the shape parameter in the tails.


title: Plotting gam model output - not component smooth functions

text: I aim to model the temperature variations in two locations in America, example below:               set.seed(10)     RandData <- runif(8760*2)     America <- rep(c('NewYork','Miami'),each=8760)          Date = seq(from=as.POSIXct("1991-01-01 00:00"),                 to=as.POSIXct("1991-12-31 23:00"), length=8760)          DatNew <- data.frame(Loc = America,                         Doy = as.numeric(format(Date,format = "%j")),                         Tod = as.numeric(format(Date,format = "%H")),                         Temp = RandData)     require(mgcv)     mod1 <- gam(Temp ~ Loc + s(Doy) + s(Doy,by = Loc) +       s(Tod) + s(Tod,by = Loc),data = DatNew)          plot(mod1,pages = 1, scale = 0)      Instead of having an output showing the component smooth functions that make up the gam I would like to plot the model output on the original x and y axis i.e. show the temperatures on the y axis. When modelling 1 location I would use something along the lines of:               pred <- data.frame(Doy = DatNew$Doy)     pred <- transform(pred, yhat = predict(mod1, newdata = pred))      However, I do not know how to achieve this if I have several locations i.e. the model depends on the location not solely on the day of year/time of day. How can this be achieved?


title: Significance of transfer entropy calculations

text: I have two time-series U and V (approximately 300 samples) having values between 0 and 5. I wanted to calculate the (first order) transfer entropy $I_{uv}$ (of U with the knowledge of V) and $I_{vu}$ (of V with the knowledge of U). Here are the values I got : $I_{uv}$ = 0.1492 and $I_{vu}$ = 0.1342 Since these values are close to each other, I started questioning their significance. Thus, I calculated the transfer entropy $I_{u-rv}$ of U with the knowledge of a randomized (in time) version of V, and the transfer entropy $I_{v-ru}$ of V with the knowledge of a randomized (in time) version of U. I calculated 500 values of $I_{u-rv}$ and $I_{v-ru}$, which gave me : mean($I_{u-rv}$ ) = 0.1218 - standard deviation = 0.0198 mean($I_{v-ru}$) = 0.1019 - standard deviation = 0.0198 The values of $I_{uv}$ and $I_{vu}$ lie close to (mean+std) so I would say there is an influence of one time serie on the other. From these information what can really be said about the two time-series ? (I will upload a picture of these time-series as soon as I figure a way to do it)


title: Does the 'fit' attribute of a gbm object contain the OOB estimates?

text: I have looked at the help (?gbm) and documentation for gbm.object, it only says that > fit: a vector containing the fitted values on the scale of regression > function (e.g. log-odds scale for bernoulli, log scale for poisson) Is it the OOB estimates like randomForest? That would seem unlikely to me since best_iter is not supplied to get it unless gbm has secretly decided which test I prefer. Is there an efficient way to get the OOB estimates at all since I've already ran 5-fold cv to get best_iter?


title: Remedial-level help needed on repeated measures analysis in R

text: I am very much a novice and really looking for some help. Following advice from a colleague I tried positing this on Stack Overflow but was told I was off topic and should post it here. I'm not trying to get someone to do my work for me - I'm genuinely looking for a bit of help to get started! I'm hoping someone will be kind enough to help me. I have somewhat foolishly offered to help my supervisor with some re-analysis of his research needed to get a paper published (it's been initially rejected because the analysis was too basic) and now realise I have no idea what I'm doing! The experimental design is based on plot treatments:               +/- prior rooting by boar     +/- fencing to exclude boar      Therefore, 4 experimental treatments, replicated 4 times at each of 5 sites giving a total of 80 plots. Data (in this case % cover of bluebells) were then collected in 3 consecutive years for each of these 80 plots (labelled 1-80 in variable "Quad"). Essentially, we want to test the 2 treatment effects (fixed factors rooting and fencing) taking account of random factor Site and the repeated measures over 3 years (so needing to account for any pseudoreplication of Quad). I've made some stabs at this using `nlme` package, e.g.               cover<-lme(fixed=Cover~Root*Fence, random= ~1|Year/Site/Quad, method="ML")     cover2<-lme(fixed=Cover~Root+Fence, random= ~1|Year/Site/Quad, method="ML")     anova(cover,cover2)       but am really not very confident in my approach. Any help would be greatly appreciated, thank you. Data headers: Data Quad: 1 - 80;   Site: 1 - 5;   Root: Rooted or Unrooted;   Fence: Fenced or Unfenced;   Year: 1 - 3;   Cover: 0 - 100 (%)


title: Removing minimally informative bits

text: Each sample in my data set is an $N$-bit bit-vector (say $N$=200). The bits in each sample are not uncorreleated within the sample. I built a matrix $S_{N \times N}$ with each $s_{ij}$ being the proportion of samples where the $i$-th bit and the $j$-th bit were identical. This is a symmetric similarity matrix (obviously, with $1$s on the diagonal). Now I want to throw some of the bits away with as little loss of information of possible.   What is the best way to choose $K$ bits, i.e. $K$ columns+rows of $S$, that when removed, the data-set retains as much of its discrimination as before. In other words, how do I choose the top $K$ bits to remove, that can be best predicted by the remaining bits for each sample.


title: Interaction/Moderation effect in ANOVA - variable is not one of the factors

text: I am running an univariate ANOVA in SPSS to test my hypotheses for my 3x2 between subjects design experiment. The factors for the ANOVA are review valence (3 levels) and presence of brand (2 levels). I want to see how these factors affect purchase probability. I also gathered data for the brand commitment of each participant in the experiment. Brand commitment is measured by giving a score on a scale from 1-9(fully disagree/fully agree) for 3 statements (taken from previous paper). Brand commitment is high/low if the scores are in the upper/lower third of the scale. I want to understand what is the interaction between brand commitment and review valence, and how it affects the DV. I predicted that brand commitment will moderate the effect of the review on the DV. My question is how can I add this brand commitment "interaction or moderation" in the ANOVA or should I use a different method to analyse it. Also, should I compute a new variable that combines the 3 statements in one new variable for the analysis. Any insights that you can provide me with, will be highly appreciated. Stefan


title: Recognize spatio-temporal patterns correlating with events

text: I am trying to recognize a spatio-temporal pattern in my spatio-temporal input sequence X. The occurrence of the pattern is roughly temporally correlated with another event E1. For example I have training data of this kind:               ...0010010000100001001001000010000001000000010010000... x1 \     ...0001001000100000000001000000000001000000100001000... x2  |- X     ...0010000000010000100010000010000010000000000010000... x3 /     ...0001000001000000000000100000000100000000000000010... E1     ------- time ------>       I have training data that contains inputs X and E1. I need a predictor for E1 given X. In the example, the predictor should learn, that E1 gets triggered, if all three components of X are active within temporal vicinity. In the training data, E1 is not super accurately timed, but can happen some time before or after the optimal point in time. Which algorithms would you suggest? PS: X might also be a real signal instead of just binary.


title: Understanding PCA - How to calculate scores

text: I'm looking for advice to whether or not the following method is good and is standard for calculating PCA of the data. So the examples that I will give will be small. Given a matrix of $A = [4, 6, 10; 3, 10, 13; -2, -6, -8]$ and Compute the Covariance matrix given as: $$A = \begin{bmatrix} 10.3& 24.6 & 35.0\\ 24.6& 69.3 & 94.0\\ 35.0& 94.0 &129.0 \end{bmatrix}$$ I then compute the Eigenvalues and Eigen vectors of this matrix. In order to calculate the PCA, I then do the following: 1) Take the square root of the eigen values -> Giving the singular values of the eigenvalues 2) I then standardises the input matrix $A$ with the following: $A - mean(A) / sd(A)$ 3) Finally, to calculate the scores, I simply multiply "A" (after computing the standardization with with Eigenvectors This would then give me the PCA scores? Is this correct or am I missing something?


title: How to set up a model for analyzing three-way contingency table?

text: As an exercise for a statistics course (not homework) I have to test various models for a given contingency table. The contingency table displays data about infant survival rates at two different clinics. There are three variables involved. $X_1$ represents the clinic ("1" or "2"). $X_2$ denotes the kind of care an individual got ("less" or "more", in this exercise). $X_3$ denotes survival ("Yes" or "No"). The exercise asks me to use glm to compute deviance of a couple of models:   * $M_1 = 1\cdot 2\cdot 3$   * $M_2 = 1\cdot 2 + 1\cdot 3 + 2\cdot 3$ and so on. I want to set this up in R but I am not quite sure how to do this. Suppose I want to set up the model $M_1$. I assume multinomial samples for a moment. How would I go about setting this up? Also, is there perhaps good literature on how to properly set up practical tests in R? The course I'm taking focusses a lot on the theoretical aspects but, as a math student, I find I'm lacking a lot of 'practical' knowledge.


title: How to handle outliers in GARCH model?

text: In GARCH model how should I handle outliers? Just remove it from my dataset and skip to next data entry?


title: What is a good modern task load measurement?

text: I will be conducting an experiment during which I need to measure mental load. After some research I decided on subjective measurements, as I can't obtain the equipment needed to attempt objective measurements. I came across the NASA TLX test which seems to be what I need, but it looks a little bit dated. Are there any more modern survey/measurement techniques?


title: How can I transform a column into date format?

text: One of my date column doesn't read in date format, but I need to use the function `WEEKDAY()` on it. How can I transform the data format of that column into date format? **UPDATE:** It is a csv file, following is one of the rows > 01OCT2007,4.92,7.75


title: How can I censor entries at failure or an exit age in Stata?

text: I am attempting to construct a survival analysis in Stata whereby subjects fail if the failure condition is met (Staph) but are also removed from the model once they are over 1100 days old. My current script is - `stset new_end_age, id(CPID) failure(Staph==1) exit(time new_end_age>=1100) origin(time birth_dt) scale(1)` Which results in the following output -               id:  CPID          failure event:  Staph == 1          obs. time interval:  (new_end_age[_n-1], new_end_age]          exit on or before:  time new_end_age>=1100          t for analysis:  (time-origin)          origin:  time birth_dt               54592  total obs.      1898  ignored because never entered     36124  obs. end on or before enter()     16387  obs. begin on or after exit            183  obs. remaining, representing       183  subjects         0  failures in single failure-per-subject data       421168  total analysis time at risk, at risk from t =         0                              earliest observed entry t =         0                                   last observed exit t =     11892      However this is incorrect as I know that there are multiple entries for individual subjects that should remain, and some of these should fail in my time period of interest. I need a command that prioritises the failure condition and if this is not met subjects are removed at `t=1100`.


title: Intuition behind Expected Value of conversion events

text: I'm trying to develop a high level model to value events in a marketing conversion funnel. To take a simple e-commerce example:   1. You start with leads in the form of ad clicks. Some % of these ad clicks result in a conversion to an email address. We'll call this 1st step conversion.   2. Of the users who submitted an email address, some % convert to a paying customer when they make a first transaction. We'll call this second step conversion.   3. Using a probabilistic model, such as Pareto/NBD, we calculate the expected Lifetime Value of our customers, i.e. the expected economic value of only those users who have passed second step conversion by making a transaction. To model the conditional expectation of the economic value of an ad click, I posit the following formula: $$ P(first step conversion=1,|X)*P(second step conversion=1|X,first step conversion=1)*E(LTV|X, second step conversion=1) $$ Where $X$ is a feature vector describing attributes of the source/customer who made the ad click. Now, the tricky part is to train a model in 3 stages to reflect the above formula:   1. Train logistic regression for 1st step conversion, using all leads as instances and feature vector as explanatory variables   2. Train logistic regression for 2nd step conversion, using subset of data, i.e. training on only those instances that passed first step conversion, as well as feature vector as explanatory variables.   3. Train probability model for Expected lifetime value, using only subset of data, i.e. training on only those instances that passed second step conversion. My thought is to integrate these three models by simply taking their product, as in the above formula. Does this intuition make sense, or am I making faulty assumptions in regards to conditional independence, etc.? I would really appreciate a sanity check, and I'd be happy to clarify further if anything is unclear.


title: Can we change the acceptance rate in random walk Metropolis algorithm by changing the parameter of the proposal distribution?

text: Can we change the acceptance rate in the random walk Metropolis algorithm by changing the parameter of the proposal distribution? Let the target distribution be $\pi$. Let $p(x_2 | x_1)$ be the proposal density for a new state $x_2$ at current state $x_1$. The acceptance rate is $$ \alpha = \min(1, \frac{\pi(x_2) p(x_1|x_2)}{\pi(x_1) p(x_2|x_1)}) $$ If I am correct, in the random walk Metropolis algorithm, the proposal density is symmetric in the sense that $p(x_2 | x_1) = p(x_1 | x_2)$, so the acceptance rate doesn't depend on the proposal density, but only on the target distribution $\pi$ to be sampled. So changing the parameter of the proposal distribution will not change the acceptance rate $\alpha$. For example, if the proposal distribution, at current state $x_1$, is a Gaussian distribution centered at the current state with a constant variance, i.e. $N(x_1, \sigma^2)$, which is by the way symmetric in the above sense, will changing the variance $\sigma^2$ of the Gaussian proposal distribution not change the acceptance rate $\alpha$? Thanks!


title: LLR with Positive and Negative Values vs. Dunning method with Entropy-based Calculation

text: Ted Dunning has a blog post about calculating G2 (aka LLR) using Entropy calculations as components. I found this really intriguing. Ted's original post: http://tdunning.blogspot.com/2008/03/surprise-and- coincidence.html And yesterday he was nice enough to help fix my Excel (and subsequent Java) implementation, correcting for +/- signs in the final formula. He even produced a corrected Excel sheet. Discussion with Ted about his formula and his answer to my post: http://math.stackexchange.com/a/693178/35616 But I think there's still a separate issue (or a misunderstanding on my part). Most LLR papers talk about negative and positive values. Things that are more prominent in A vs. prominent in B should change signs, so you know which direction the change is in. A couple examples of this discussion: Can the log likelihood ratio for a simple vs simple hypothesis take a negative value? and http://scg.unibe.ch/archive/papers/Kuhn09aLogLikelihoodRatio.pdf Chi Squared is often compared to G2 (LLR). Although Chi Squared is always positive, since it's squared!, apparently the "G2" test doesn't hold to that, the "2" perhaps being a bit misleading. Looking at Dunning's calculations, I don't see how the sign can ever flip? His contingency tables:                             Corpus A  Corpus B   Row Totals     Target Word     k_11     k_12      totalRow1     Other Words     k_21     k_22      totalRow2     Column totals   col1     col2      grandTotal      You then calculate Entropy for row totals, column totals, and the overall 4 k cells:   * H_rowTotals   * H_colTotals   * H_k These are later combined in the final formula, along with grandTotal. BUT the signs don't change when a word moves from column 1 to column 2:   * For row1, if you transpose k_11 and k_12 (and assuming row 2 is unchanged), then row 1's row's total stays the same, and therefore H_rowTotals doesn't change sign.   * Transposing k_11 and k_12 does change the order of the column totals, BUT the calculation of H_colTotals isn't impacted by order, so the sign doesn't change.   * And then when calculating H_k, it also doesn't care what order the k cells are in, so it also doesn't change sign.   * And of course the grand total is always >= 0. So transposing k_11 and k_12 can't change the sign of H_rows, H_cols, H_k, nor grandTotal, which are the only 4 variable inputs into the final formula, so the final formulate can't change signs. Ted was nice enough to upload a revised Excel sheet, but it actually demonstrates this non-sign change: https://dl.dropboxusercontent.com/u/36863361/entropy-and-LLR-suspect-gist.xlsx But whether you do:                             Corpus A  Corpus B     Target Word      10        0     Other Words       0       10      Or:                             Corpus A  Corpus B     Target Word       0       10     Other Words      10        0      You still get: +27.7 The only theory I can think of is maybe there's some alternative definition of Entropy that is order dependent; or something that approximates it. Another theory I discarded was that maybe you change the signs of one of the columns, but that won't work since it would generate negative numbers in the probability calculation which would then be invalid as inputs to the log function. To be clear, I'm super grateful for the previous help, but this sign change keeps bugging me. Also, posting this on the Stats Stack Exchange site, vs. the general Math site, since I think it's more specific.


title: Terminology problem: "model selection" is the same as traning ?

text: In machine learning we have the following problem: Choosing the optimal model (or training): $$ f^* = \arg\min_{f \in \mathcal{F}} \sum_i l(x_i,y_i) $$ Is the term `model selection" always`exactly" referring to to this? Or something else?


title: How to train in models, with efficient inferences, like belief-propagation ?

text: There are many papers that are devoted to efficient inference in graphical models. Though many of these paper don't explicitly talk about the learning (training, etc) problem. For example: http://videolectures.net/mlss09uk_minka_ai/ I am a little confused on how these models being trained. I thought they are probably doing an EM-like algorithm, i.e.   1. Inference (and calculating all marginals, using VB or EP)   2. Maximizing the likelihood using some blackbox optimization toolboxes using the marginals in the previous case For example, consider different variants of Belief-Propagation. There are HUGE number of variants for BP, but how a graphical model could be trained? Any comments?


title: How to evaluate/validate clusters using multiple clustering methods

text: From some reading I did online, I understand that there are various methods for determining "similarity" used by different clustering algorithms. I am curious if it is good practice to run multiple clustering algorithms/methods (i.e Hiearchical w/ Ward, single linkage, centroid, etc or maybe even K-means) on a dataset and if there is some automated way to to get a "consensus" of clusters. In other words to get some sense of confidence that the right items are clustered together. Items that tend to cluster together using various methods would be considered valid. For example in my example below G and Z tend to cluster together using multiple methods as do S and F. Label = what I am clustering; X & Y are my variables I use to cluster; Cluster1-3 are the results of three clustering algorithms. ![enter image description here](http://i.stack.imgur.com/7m55l.png) _Edit: I removed a side note I had here regarding how large the actual data set I plan to use might be so as not to detract from the main questions._


title: In weighted least squares, how do I weight the residuals to get an accurate "z score"

text: I am regressing spreads in yield curves in certain countries, in the chart below, the Spanish 2-5-10 spread, against the Italian 2-5-10 spread. ![enter image description here](http://i.stack.imgur.com/mteWO.png) I want recent data to count more, so I weight the inputs using a decay weighting scheme with a 1 year halflife. A "simple" regression line and the weighted regression line are shown. I want to calculate the perpendicular distance of the current point (green) from the regression line, in number of standard errors. In the unweighted regression, in R, I will simply say:               > l <- lm(SP ~ IT, data = ss)     > last(l$residuals) / sd(l$residuals)     2013-10-28      -0.1817122       Which gives -.18 standard errors away from the regression line. How do I do this same analysis for the weighted regression though? I am sure the following is incorrect:               > decay     function(len, halflife, sumone = TRUE) {     #function generates an exponentially decaying series         t <- len:1 # generate a series of numbers reverse order so biggest weights last         lambda <- log(2) / halflife #figure out the lambda for the halflife         w <- exp(-lambda * t) #create the weights series           if(sumone) w <- w / sum(w) #normalise sum to 1 if necessary         return(w)      }     > d <- decay(nrow(ss), 260)     > ld <- lm(SP ~ IT, data = ss, weights = d)     > last(ld$residuals) / sd(ld$residuals)     2013-10-28      -0.3667876       I should surely weight the residuals somehow, before doing the above, is that correct? Could I for example take the weighted standard deviation of the residuals that is:               > last(ld$residuals) / wt.sd(ld$residuals, d)     2013-10-28        -0.39717       where my wt.sd function looks like this:               > wt.sd     function (x, wt) {         return(sqrt(wt.var(x, wt)))     }          > wt.var     function (x, wt) {         s = which(is.finite(x + wt))         wt = wt[s]         x = x[s]         xbar = wt.mean(x, wt)         return(sum(wt * (x - xbar)^2) * (sum(wt)/(sum(wt)^2 - sum(wt^2))))     }      Basically, I want to know how to find the distance from the weighted regression line, in standard errors, accounting for the weights.


title: What is a reasonable process to understand a collection of data?

text: Can someone provide their thoughts on a structured process one might go through to understand a collection of data. The scenario is: you've been given a set of data (features and observations - with descriptions) and been told to "tell me what kind of interesting things this data can tell me". I.e., what are interesting questions that this data can answer. The meaning of "interesting" is certainly subjective. This appears to be classical unsupervised learning. **My initial thoughts:**   1. Cluster all pairs of variables to see interesting clusters    2. Run PCA find high-variance groupings Is there a general "how to understand a set of data" process that you've found successful? Thanks


title: Can I prove a curvilinear relationship when the linear independent variable is not significant

text: I am investigating a curvilinear effect between X and Y by using a hierarchic regression analysis. To test for curvilinear effects, the squared term for X was computed (I mean center also variable X). In model 1, the control variables were entered. In Model 2, X (linear) was entered. In Model 3, X (quadratic) was entered. In Model 2, X linear is significant. When the squared term is entered in Model 3, the quadratic term is significant but the linear term is not. Does this prove a curvilinear effect? Or is it essential that in Model 3 both (linear and quadratic) are significant? When I do not mean center the independent variable, Model 3 stated X linear and X quadratic significant. The problem here is multicollinearity issues.


title: Variance of nominal variable?

text: Note that the only measure of location in nominal variable is mode. Is there also any measure of variance? Obviously, some nominal features are more different than other (color of eyes vs. presence of some rare disease).


title: What does it mean to have an F-value of .000 and Sig. of 1.000 (in SPSS)?

text: We conducted this very simple 2x2 Between Subjects Factorial Design Experiment, where the DV were their answers to a classic insight problem. Since there is only one correct answer to the problem, the values for the DV were either 1 if the answer was correct and 0 if it was not. The results I got was an F-value of .000 and Sig. (p value) of 1.000 for one variable and the interactio n. Is there something wrong with my data? Perhaps the DV shouldn't have been dichotomous? What do you suggest I do?


title: Nonparametric alternative to ANOVA for testing the difference between several brand preferences

text: I have collected a data from about 200 customers on their brand preference of 5 international fast food chain shops. Each customer was asked to rank all the 5 brands from highest to lowest (5 for the most preferred and 1 for the least preferred brand). The brand rankings will vary from one person to another. I want to test if there are significant differences between the brand preferences. At a fist glance ANOVA seems to be a good option. We could take the average of the ranks for each brand and then perform ANOVA. But the problem is that, the ranks can not be treated as scale, because the spaces among the ranks may not be equal. (e.g. Brand-2 and Brand-3 may differ a lot to a customer according to his preference than Brand-1 and Brand-2). So, what can be the best alternative to test if the brand preferences are significantly different? I think I should adopt some non-parametric test like test of several medians instead of mean. What would be your suggestion regarding that? Will Kruskal-Wallis be a valid test here? Would you please give me some indications as to how I should perform whatever you suggest in either R or in SPSS?


title: Can we see shape of normal curve somewhere in nature?

text: I do not want to know if some phenomena in nature have normal distribution, but whether we can somewhere see shape of normal curve as we can see it for example in Galton box. See this figure from Wikipedia. ![enter image description here](http://i.stack.imgur.com/9KJbC.jpg) Note that many mathematical shapes or curves are directly seen in nature, for example golden mean and logarithmic spiral can be found in snails. First naive answer is whether nonskewed hills would often "fit" normal distribution :-).


title: is this two dimensional markov chain correct for this queueing system?

text: the problem that i have two single server station with no queuing space a customer goes to station 1 if it is available else it goes to station 2 if it is available or it will be lost output from station 1 goes to station 2 if it is empty or blocking occurs and will stay in station 1 til station 2 become available. if customer visits station 2 first he will be served for two consecutive service times and then leave the system . so to model this i consider that the first station will have 3 states 0(empty),1(customer get service),b(customer finish the service and blocked in station 1) and 4 states for station 2 : 0(empty),1(customer from station 1 get service at station 2),a(customer initially go to station 2 and get his first service time ) , B(customer initially go to station 2 and get his second service time) so is this model correct or it can be simplified more better .![enter image description here](http://i.stack.imgur.com/V3X3y.jpg)


title: One sample median test: Wilcoxon, sign test or chi squared

text: I would like to understand the differences between these tests for one sample median test and when to use each of them. I have been searching information and this is what I got so far:   * Sign test (in R sign.test) does not need a probability distribution symmetric but Wilcoxon (wilcox.test) does.   * Wilcoxon is more powerful. To check the symmetry in the prob distribution, I through an histogram right? What are the assumptions for each of them? When should I use them? I am pretty lost here so any help is welcome.


title: In binary classification, what label imbalance is acceptable to not to care?

text: I have a pretty large set of binary-labeled data (~20K rows, ~100 columns). Approximately 20% of the data is "positive" and 80% is "negative". In my case every classification error has the same meaning. Is this misbalance big enough in order to take special measures such as differential weighting and/or differential sampling. Is there any rule of thumb? Or maybe a publication? Does the answer to this question depend on the classification algorithm?


title: Likelihood Ratio for two-sample Exponential distribution

text: Let $X$ and $Y$ be two independent random variables with respective pdfs: $$f \left(x;\theta_i \right) =\begin{cases} \frac{1}{\theta_i} e^{-x/ {\theta_i}} \quad 0<x<\infty, 0<\theta_i< \infty \\ 0 \quad \text{elsewhere} \end{cases} $$ for $i=1,2$. Two indepedent samples are drawn in order to test $H_0: \theta_1 =\theta_2 $ against $H_1 : \theta_1 \neq \theta_2 $ of sizes $n_1$ and $n_2$ from these distributions. I need to show that the LRT $\Lambda$ can be written as a function of a statistic having $F$ distribution, under $H_0$. * * * Since the mle of this distribution is $\hat{\theta}=\bar{x} $, the LRT statistic becomes (I am skipping a few tedious steps here): $$ \Lambda =\frac{\bar{x}^{n_1} \bar{y}^{n_2} \left( n_1+n_2 \right)}{n_1 \bar{x}+n_2 \bar{y}}$$ I know that the $F$ distribution is defined as the quotient of two independent chi-square random variables, each one over their respective degrees of freedom. Additionally, since $X_i,Y_i \sim \Gamma \left( 1,\theta_1 \right)$ under the null, then $\sum X_i \sim \Gamma \left(n_1 ,\theta_1 \right)$ and$\sum Y_i \sim \Gamma \left(n_2, \theta_1 \right) $. But how can I proceed from here? Any hints? Thank you.


title: Study replication from a Bayesian point of view

text: Consider following situation: Study A compares two groups and finds a mean difference with an effect size of d = .8 (p<.05). Study B is a direct replication, and finds an effect in the same direction of d = .3 (p < .05). From an NHST perspective, there is no single/ simple answer to the question "Does Study B replicate the results of Study A?". Asendorpf et al. (2012) summarized it as follows: > "Given that no single approach to establish replicability is without limits, > however, use of multiple inferential strategies […] is a better approach. In > practice, this means summarizing results by answering four questions: (a) Do > the studies agree about direction of effect? (b) What is the pattern of > statistical significance? (c) Is the effect size from the subsequent studies > within the CI of the first study? (d) Which facets of the design should be > considered fixed factors, and which random factors?". That means, there is no single answer for the question "Has this study been replicated?" (see also Sanjay Srivastava's blog post What counts as a successful or failed replication?). Now I ask myself: Is there probably a simple(r) answer from a Bayesian point of view? Kruschke (2010) describes a "cumulative replication probability", which takes the actual posterior distribution from Study A as the prior for simulated data, which then gives an answer on the probability of replicating a certain decision of Study A (e.g. a Bayesian model comparison). In my scenario, however, we have already two Studies conducted. **Now here's my question: Is there some way to feed Study A's posterior as information into Study B's results and come up with an answer to the question "Is Study B a replication of Study A?"? How could such an answer look like?** (As I am interested but so far unexperienced in Bayesian stats, my description of the situation might be utterly wrong ...) * * * Asendorpf. (2012). Recommendations for increasing replicability in psychology. European Journal of Personality. Kruschke, J. K. (2010). Bayesian data analysis. Wiley Interdisciplinary Reviews: Cognitive Science, 1(5), 658–676. doi:10.1002/wcs.72


title: Are random variables correlated if and only if their ranks are correlated?

text: Assume $X,Y$ are continuous random variables with finite second moments. The population version of Spearman's rank correlation coefficient $ρ_s$ can be defined as the Pearson's product-moment coefficient ρ of the probability integrals transforms $F_X(X)$ and $F_Y(Y)$, where $F_X,F_Y$ are the cdf's of $X$ and $Y$, i.e., $ρ_s(X,Y)=ρ(F(X),F(Y))$. I wonder whether one can generally conclude that $ρ(X,Y)≠0↔ρ(F(X),F(Y))≠0$? I.e., do we have linear correlation if and only if we have linear correlation between the ranks? **Update:** In the comments two examples are given why $\rho(F_X(X),F_Y(Y))=0\rightarrow \rho(X,Y) = 0$ is not true in general, even if $X$ and $Y$ have the same distribution. So the question should be reformulated as $\rho(X,Y) = 0 \rightarrow \rho(F_X(X),F_Y(Y))$? It is also of great interest to me whether this is true/false if $X$ and $Y$ have the same distribution. (Note: If $X$ and $Y$ are positively quadrant dependent, i.e., $δ(x,y)=F_{X,Y}(x,y)−F_X(x)F_Y(y)>0$ then Hoeffding's covariance formula $Cov(X,Y)=∫∫δ(x,y)dxdy$ yields that $ρ(X,Y)>0$ and $ρ(F(X),F(Y))>0$.)


title: How can the number of connections be Gaussian if it cannot be negative?

text: I am analyzing social networks (not virtual) and I am observing the connections between people. If a person would choose another person to connect with randomly, the number of connections within a group of people would be distributed normally - at least according to the book I am currently reading. **How can we know the distribution is Gaussian (normal)?** There are other distributions such as Poisson, Rice, Rayliegh, etc. The problem with the Gaussian distribution in theory is that the values go from $-\infty$ to $+\infty$ (although the probabilities go toward zero) and the number of connections cannot be negative. Does anyone know which distribution can be expected in case each person independently (randomly) picks-up another person to connect with?


title: Should the intersection of 2 event sets be on the same sample space?

text: I'm trying to understand a stochastic process and consecutive events in time. Should the intersection of two events be on the same sample space, on a combination of _simultaneous_ sample spaces, or a combination of _consecutive_ sample spaces? I have a question about $A\cap B$. Should the sets $A$ and $B$ be defined on the same sample space? If not, do these 2 events have to happen at the _same time_? This question arises from this problem. Consider a discrete stochastic process with 2 states, $\{0, 1\}$ and 2 stages $t=1, 2$. We can define 2 sample spaces: for $t=1$, $\Omega_1=\{0,1\}$ and for $t=2$, $\Omega_2=\{0,1\}$.   1. Can we define a _total sample space_ $\Omega_1\times \Omega_2$ even though the 2 events do **not** occur at the same time, but one after the other (consecutively)?   2. Assume that the probability of an event happening in $\Omega_1$ **affects** the probability of events in $\Omega_2$. In other word the **stationary transition** probability has been given, and also assume that the initial probability of states at $t=1$ has been given. Then what is the probability of event $[0,0]$ in $\Omega_1\times \Omega_2$?    I would say, $\tfrac{\mathbb{P}[0 \mid t=2]}{\mathbb{P}[0 \mid t=1]}$.   My answer comes from this equation $\mathbb{P}(A\cap B)=\mathbb{P}(A)*\mathbb{P}(B/A)$.   But if $ P(A\cap B)$ should be defined on the same sample space,or 2 sample spaces both at the same time and not consecutively, my answer is wrong.   The other answer can be $\mathbb{P}[0 \mid t=1]* \mathbb{P}[0 \mid t=2]$ (by having $\mathbb{P}[0 \mid t=1]$ value and transition probability, we can calculate $\mathbb{P}[0 \mid t=2]$ and then multiply by $\mathbb{P}[0 \mid t=1]$).   3. Is there any difference between **transition** probability in a process (like Markov) and **conditional** probability that we define in set theory? I mean differences in the underlying sample spaces?


title: Any out source outlier (anomaly) detection package for Weka?

text: There is KDD'98 cup data waiting for me to run some anomaly detection algorithms but Weka does not include any of them, natively. Is there any 3rd part package that can be integrated to weka?


Can the log likelihood ratio for a simple vs simple hypothesis take a negative value? Since it approximately follows a chi-square distribution for large sample sizes and since a chi-square distribution takes only positive values do we reject or accept the null hypothesis if a negative value is obtained for the log likelihood ratio test under a particular simple vs simple hypothesis? The log likelihood ratio test statistic in this case is $$-2Log\Lambda$$


title: 2x2 ANOVA interactions

text: I'm doing a masters project where I look for effects of font (fluent, disfluent) on two populations' (control, dyslexic) reading comprehension. SPSS only gives me a significant main effect for population, not font, and no significant interaction between population and font. However, when eyeballing the graph, I can see a clear interaction (lines actually cross): as the controls' comprehension gets better on disfluent font texts, the dyslexics' comprehension gets worse (lines goes in the opposite direction) on the disfluent font material. Why is SPSS not reporting significant results for the interaction here? How do I go about reporting my findings?


title: A mixture of mixed-model regression and fixed-effects-only regression: what is it exactly?

text: I have a model with a continuous dependent variable (DV) and 9 independent variables (IVs). My model is a mixed-model, because 7 of the 9 IVs are nested within the two levels of another IV (lets call it X1) which is applied in a repeated-measure design to the same group of subjects (zero at first trial and 1 at the second trial). So I have done the nesting according to the SubjectID. X1 is either 0 or 1 for each subject, but other IVs _except another IV_ are similar for each subject. As I stated right now, a remaining IV (lets call it X2) is **not** nested. As a better explanation, X2 does not follow a repeated- measures pattern; and thus it is not necessarily the same at the first and second trials in each subject (it might be the same for some subjects, and there is some correlation, but it is not intentionally repeated like the other 8 IVs). Now I should build a mixed-model in which the nesting is done according to subjects/trials/X1 (SubjectID). The effect of X1 is being already added as a random factor. Now, all the variables except X2 are properly modeled. However, I am not quite sure about X2. I know I can (and perhaps _should_ ) add X2 as a fixed-effect to the mixed-model. But that way, I will lose some of good information because the model will treat the X2 as a nested variable, while it is not actually nested. I also can run a fixed-effect regression with only X1 and X2 to assess the role of X2 regardless of the nesting. That is somehow good, but that way I cannot add the nested variables (because they will cause atomistic fallacy). So I want something in between. I don't know how to build a mixed model that has the advantages of a mixed-model (has random factor) and also when it comes to X2, does not lose some power (does not consider the nesting according to SubjectID for X2). So I need a model that in most of parts considers SubjectID to nest the data (all the IVs except X2) and at some part (in the case of X2), does not consider SubjectID. Is it possible? Do you have any other suggestions? ...On a second thought, it appears to me that perhaps X2 can be considered repeated-measures as well. It is true that this variable is not like the other ones, exactly repeated within each subject, but there is still some correlation. So perhaps it is good to run a _routine_ mixed-model regression without worrying about losing some power. What do you think? Many thanks.


title: Things to consider before summing binary variables to create a total score

text: I would like to reduce the amount of data I have before conducting an analysis. My data consists of several sets of questions assessing comprehension of various topics (e.g. 4 questions assess comprehension of one concept, 4 questions assess comprehension of another concept and so on). I would like to sum the questions for each concept before running further analyses. But, I am not sure whether it is ok to do this without first conducting a factor analysis or some other reduction technique? If I do go down the road of using a multivariate method like FA to reduce my data, then I run into another problem. My data is categorical and so I don't know whether ordinary FA is suitable to use. I see there is the option of Categorical Principal Components Analysis. But, I am not sure if this is suitable. How do you think I should proceed?


title: Sampling probability

text: Say I start off with N single cells in a sample. These are then allowed to grow untill there are x progeny from each individual cell. How many cells must I take so that I have at least one progeny from each the N original cells in the sample?


title: How to rewrite a sum of probabilities formula as multiplications?

text: I have an equation like that: $p(r|s)= \frac{p(s,r)}{p(s)}=\frac{ \sum_{w,c} p(c,s,r,w)}{\sum_{w,c,r} p(c,s,r,w)} $ I am new to probability and I want to learn that how can I write sum formula as multiples of probabilities? **PS 1:** I want something like: $p(r|s)= \frac{ p(w,r|c,s)p(c,s)p(w,r|~c,s)p(~c,s)...}{ p(s|w,r,c)p(w,r,c)...} $ **PS 2:** Here is my Bayesian network: ![enter image description here](http://i.stack.imgur.com/hvXcg.png) **PS 3:** Because of conditionally independence I have that equations: $p(w|r,s,c) = p(w|r,s)$ $p(r,s|c) = p(r|c)p(s|c)$ $p(r|s,c) = p(r|c)$ $p(s|r,c) = p(s|c)$


title: Apply different transformations to different samples for parametric test?

text: Suppose I have two samples from a population collected from different times, but their distribution looks different (e.g. sample one looks skewed right, while sample two looks bimodal). Can I apply a transformation to sample one, and a different transformation to sample two to make them both normally distributed so I can run parametric tests comparing the two samples?


title: MLE for Naive Bayes in R

text: I am using the `naivebayes` function of the e1071 library. Some example commands are:               model = naiveBayes(Species ~ ., data = iris)     pred = predict(model, iris[,])      My question is: how can I obtain the maximum likelihood estimate for the conditional probability distribution of this model?


title: Testing nonlinearity in logistic regression

text: One of the assumption of logistic regression is the linearity in the logit. So once I got my model up and running I test for nonlinearity using Box-Tidwell test. One of my continuous predictors (X) has tested positive for nonlinearity. What am I suppose to do next? As this is a violation of the assumptions shall I get rid of the predictor (X) or include the nonlinear transformation (X*X). Or transform the variable into a categorical? If you have a reference could you please point me to that too? Many thanks


title: Which method of implementing the Brown's linear exponential smoothing is correct?

text: I am trying figure out what is the difference between Brown's linear model for double exponential smoothing and Holt's model. So the differences can be implemented into a Holt model using if statements. Aside from having alpha for both the trend and smoothing component. I am having a problem where the forecasting occurs when the level equation is added to the trend equation to forecast: IBM SPSS show that when forecasting k steps ahead ((k-1)+(alpha^-1)) is used. But in other books and articles it shows that for the forecast equation for the level and trend is changed to 'at' and 'bt' before being added together. Are the SPSS equation and the other one the same? (I am asking because it is easier to implement one over the other) Also what is recommended for the initial points? And did I miss any differences?


title: Neural networks in multi-class setting: How to train/test with cross-validtion; how to evaluate?

text: I am currently working on a project where neural networks are used for email categorization. There are two things where I still have some problems understanding the general approach: cross-validation and evaluation using micro/macro averaging. But first some background: Let c be the number of categories, e.g. 15; n is the number of manually categorized emails, e.g. 1000. Also, I have a (relatively large) number of different parameter configurations for the trained network and the construction of feature vectors, say p = 200. Feature vectors are binary, based on a simple "bag-of-words" model of the emails. As a result, I would like to get the best parameter configuration for training and the neural networks with the best precision, recall and F1 for this particular configuration. These 3 networks will then be used for the actual classification task of an unseen dataset. So, now for the two questions: 1) Cross-validation It seems like there are two approaches that I could take for training/testing: a) For each parameter configuration p, split my data set into training and test data, e.g. 80:20. Then, train 5 network instances using 5 different 80:20 partitions. Evaluate each trained network instance against the 20% test data and select the best network out of 5 as the best one for this particular configuration. In the end, pick the winner out of all configurations based on its evaluation against the test set in its training run. Or: b) Reserve a proportion of the data as held-out data, e.g. 10%. Split the remainder of the data, e.g. 80% for training, 10% for test, so that I have a 80:10:10 partition. Then, for each config, run training and test on the 80:10 set, i.e. 9 partitions. Pick the best performing network out the 9 as the winner for this config. In the end, evaluate the winner for each config against the held-out data and select the best performing based on this. Which of these approaches is the recommended one and why? I'd also appreciate pointers to any papers that deal with this question. 2) Evaluation in multi-class setting Each email in the pre-categorized training set is assigned to one primary category and might also have one or more secondary categories. The network classifier might assign any number of categories between 0 and c. I would like to test both multi-class and mono- class classifications in this multi-label classification setting. I have read about micro and macro-averaging, but there are two things that aren't totally clear: a) Are the calculations of precision, recall and F1 in the multi-label cases essentially based on c binary classifications, i.e. one- vs-all? Or are the values calculated based on an c*c confusion matrix? b) How are micro- and macro-averaged evaluation meassures calculated differently in the mono-class evaluation (where only the primary category is considered) vs. multi-class evaluation (where both primary and secondary categories are considered)? This has been a long post for my very first question in this forum, so thanks a lot for reading this far and, perhaps, thanks a lot for any help and hints you might be able to provide. Cheers, Martin


title: How to statistically determine if something appears in the beginning of a vector

text: I have multiple vectors of letters, in which an entry can appear multiple times. The vectors are of variable length.               ACBDEABAEE     ACDE     ACBDBBDDE      I want to determine statistically if `A` or `C` occurs mostly at the beginning of the list, whereas `E` occurs towards the end and `D` occurs in the middle. Is there a way to do this?


title: Marginalising probability distribution

text: I am reading a book on probability theory and there is an example to denote conditional independence. A joint distribution is given as follows $$ p(a, b, c) = p(a) p(b) p(c|a,b) $$ Now, assuming none of the variables are observed, the point is to show $a$ and $b$ are conditionally independent. They marginalise over $c$ to get $$ p(a, b) = p(a) p(b) $$ Is this because when you sum over $c$, $p(c|a, b)$ is a valid distribution and sums to 1?


title: Continuous dependent variable that almost looks like a categorical variable

text: I recently conducted an experiment where people had to decide how much money to take from another participant. People could choose to pick any value between 0 and 100 percent. It turns out that around 40% of the participants decided to take nothing, 40% to take everything, and only 20% took a positive amount different from 100. So, in effect, the data that I get looks categorical although it is continuous, (given that people could have chosen any amount between 0 and 100). I was wondering what kind of regressions I could carry out in this situation? The residuals in an OLS regression are highly non-normal, so I don't think that an OLS regression would be a good idea. Moreover, I'm not sure it makes sense to carry out logit regressions, given that the variable is not really categorical. Thank you very much for your help!


title: extract residuals from adonis function in vegan

text: I am using the adonis function in the vegan package to determine effects of different environmental factors in forest plant community composition in different regions. I would like to first use adonis to remove the region effect, this is, to fit a model like adonis_region<- adonis(community ~ region,permutations=999,method="bray") Where community is a presence/absence data matrix of species in forest patches belonging to different regions. Then I would like to use the residuals of this model as the response variable for some other analyses. I know I could use strata to model region as a block variable, but this does not interest me, as afterwards I want to perform another kind of analysis where I would like the region effect to be removed. My problem is that I cannot figure out how to get residual values from the adonis model. Any hint of some other kind of multivariate analysis that could solve this problem is very welcome. Thank you very much!


title: Mutual information between sets of variables

text: I found in this paper (PDF) a mutual information between groups of variables such as these: >   * I(X;Y|F) >   * I(F;Y) >   * I(F;Y|G) > where F is a set of variables, and G is a set of sets of variables. I am neither statistician nor mathematician, but I have a fair knowledge about mutual information, and I do know how to calculate mutual information between two variables I(X;Y) if we have $n$ samples of X and Y. Can anyone help me do the same for the formulas above? Thanks.


title: Is it possible to estimate the convex combination of parameters in the IRLS-Framework?

text: Suppose I want to estimate the parameter $\mathbb{E}(Y)=\mu \ge 0$ with $\mu = a(\alpha)\mu_1(\gamma_1) + \Big(1-a(\alpha)\Big)\mu_2(\gamma_2)$ where $a(\alpha)\in (0,1)$ Using the usual link-functions for $a,\mu_1$ and $\mu_2$, I'll get $a(\alpha) = \frac{\text{exp}(\eta_\alpha)}{1+\text{exp}(\eta_\alpha)}$ and $\mu_i = \text{exp}(\eta_{\gamma_i}), \ \ i=1,2$ where $\eta$ corresponds the the linear-predictor. For the sake of simplicity, I suppose that my mean $\mu$ is a poissonian rate, i.e., $\mathbb{E}(Y)=\mathbb{Var}(Y)=\mu$ In the IRLS-Framework, the weight-matrix $\boldsymbol{W}$ is defined by $\boldsymbol{W} = \boldsymbol{V}^{-1}\Big(\frac{\partial\mu}{\partial\eta}\Big)^2$ Since $\boldsymbol{V}^{-1} = \frac{1}{\mu} \ \ $ and $ \ \ \Big(\frac{\partial\mu}{\partial\eta}\Big)^2 = (\mu_1 - \mu_2)a\cdot (1-a)$ The weight-matrix for alpha, i.e., $\boldsymbol{W}_\alpha$ is $\boldsymbol{W}_\alpha = \frac{(\mu_1 - \mu_2)a\cdot (1-a)}{\mu}$ Now I've to derive the other weight-matrices $\boldsymbol{W}_{\gamma_1}$ and $\boldsymbol{W}_{\gamma_2}$ as well as the linerized forms of the link functions applied to $Y$, i.e., $\boldsymbol{z}_\alpha, \boldsymbol{z}_{\gamma_1}$ and $\boldsymbol{z}_{\gamma_2}$. **My question is, if this approach is feasible in this context?** Before, I used the classical Newton-Raphson-Like algorithms provided by the R-function `optim()`. But I would like to try something different like the IRLS or EM approach.


title: How to test unit root in a timeseries with unknown structural change?

text: Is there a unit root test that works with structural change? I only found Zivot-Andrew test but I would try another test, I do not know if the timeseries has a structural change, so I need a test that works with A break at unknown time.


title: Circular/elliptical tests for datasets with magnitude and direction

text: I am trying to analyse (using R) a set of r(theta) data to see if the magnitude (r) is dependent on direction (theta). I have looked at circular statistics but these only seem to deal with the distribution of the observations theta(i) around the circle. Our null hypothesis is r=1 for all theta, vs its not. If H0 is false then we would like to be able to find the principal directions involved. From theoretical considerations we expect r(180)>1 and r(0)<1 with r(90)=r(270)=1 Some way of graphically representing this would be useful. The data comes from the interaction of water waves and tidal currents with zero degrees representing the waves travelling in the same direction as the current.


title: Comparing order of observed outcomes to predicted order

text: **N.B.** , The questions require a bit of setup, which I've tried to present abstractly; please feel free to comment on anything that needs clarification. I have problem where I am predicting the value of several actions for an ensemble of systems that satisfy a particular set of constraints, then observing the each of action values (via simulation) for sample instances of that system. I'm presenting the results as a series of boxplots where the observed distributions are vertical boxes spreading in the y-direction, and are located horizontally by their predicted values. Each plot in the series compares a different prediction method. Question 1: I think this provides a pretty good visual argument about the relative worth of the prediction methods, but it's not exactly quantitative. I need a statistic that compares the different methods ability to accurately order the outcomes, including accounting for "close" predictions mis-orders when outcomes are also close being not as "bad" as "not close" predictions mis-orders. Suggestions? Question 2: I also have data for a large sample of different constraint sets - i.e., many different ensembles; in this case, each system has a different ordering for actions (not just by prediction method). Any recommendations for visualizing the results comparing the different prediction methods?


title: SAS Proc Mixed model interpretation

text: I'm fitting a linear mixed model by SAS. There are 596 sectors and 8489 subjects. (each sector contains 10~15 subjects). Each subject is measured at most 6 times, so the total number of observation is 50043. SAS code are as follows.               proc mixed data=work.dat2 covtest method=ml maxfunc=1000 ;        class group_k sectorid childuid;        model laz=group_k x1 x2 x4 x6 x1_k x2_k x4_k x6_k                   / solution cl outpm=out;              random sectorid;              repeated / subject=childuid type=cs ;              run;      One of the result tables is as follows. The first table is without the option 'type=cs' (default type=VC), and the second table is with the option 'type=cs'. ![enter image description here](http://i.stack.imgur.com/bUPaz.png) ![enter image description here](http://i.stack.imgur.com/lLdpD.png) I don't understand the relation between the SAS result table and the theory presented here http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect003.htm The covariance matrix for the default 'variance component' is $ \sigma_k^2 1(i=j) $ and that for the type 'compound symmetry' is $ \sigma_1 + \sigma^2 1(i=j) $ In the two tables, I don't understand which corresponds to $ \sigma_k, \sigma_1, \sigma $ respectively.


title: Distinguishing two distributions

text: In my application I'm receiving a continuous stream of 50 changing numbers, mostly these samples are contained with a nice Gaussian distribution, but now and then there are two+ distributions within the samples. The samples are always positive and within the range of 0 to 1. If there is more than one distribution present, the distribution I'm always interested has the lowest mean. So ideally I want to identify the lower bound and roughly the mean of the main distribution. The lower bound is straight-forward, its always the minimum value of the samples, but I'm stumped as to how find roughly the mean of the main distribution when there is more than one distribution present. The really tricky thing is that its a real time system and I have very little time to process the samples before the next set arrives, so other than finding the min, max, mean and at a push median I don't have time to carry out much more processing. Anyone have any suggestions for how I could achieve my goals?! Many thanks and sorry for the stats-newbie nature of my question


title: Interpreting coefficient in a linear regression model with categorical variables

text: I will give my examples with R calls. First a simple example of a linear regression with a dependent variable 'lifespan', and two continuous explanatory variables.               data.frame(height=runif(4000,160,200))->human.life     human.life$weight=runif(4000,50,120)     human.life$lifespan=sample(45:90,4000,replace=TRUE)     summary(lm(lifespan~1+height+weight,data=human.life))          Call:     lm(formula = lifespan ~ 1 + height + weight, data = human.life)          Residuals:     Min       1Q   Median       3Q      Max      -23.0257 -11.9124  -0.0565  11.3755  23.8591           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) 63.635709   3.486426  18.252   <2e-16 ***     height       0.007485   0.018665   0.401   0.6884         weight       0.024544   0.010428   2.354   0.0186 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 13.41 on 3997 degrees of freedom     Multiple R-squared: 0.001425,   Adjusted R-squared: 0.0009257      F-statistic: 2.853 on 2 and 3997 DF,  p-value: 0.05781      In order to find the estimate of 'lifespan' when the value of 'weight' is 1, I add (Intercept)+height=63.64319 Now what if I have a similar data frame, but one where one of the explanatory variables is categorical?               data.frame(animal=rep(c("dog","fox","pig","wolf"),1000))->animal.life     animal.life$weight=runif(4000,8,50)     animal.life$lifespan=sample(1:10,replace=TRUE)     summary(lm(lifespan~1+animal+weight,data=animal.life))          Call:     lm(formula = lifespan ~ 1 + animal + weight, data = animal.life)          Residuals:     Min      1Q  Median      3Q     Max      -4.7677 -2.7796 -0.1025  3.1972  4.3691           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 5.565556   0.145851  38.159  < 2e-16 ***     animalfox   0.806634   0.131198   6.148  8.6e-10 ***     animalpig   0.010635   0.131259   0.081   0.9354         animalwolf  0.806650   0.131198   6.148  8.6e-10 ***     weight      0.007946   0.003815   2.083   0.0373 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 2.933 on 3995 degrees of freedom     Multiple R-squared: 0.01933,    Adjusted R-squared: 0.01835      F-statistic: 19.69 on 4 and 3995 DF,  p-value: 4.625e-16      In this case, to find the estimate of 'lifespan' when the value of 'weight' is 1, should I add each of the coefficients for 'animal' to the intercept: (Intercept)+animalfox+animalpig+animalwolf? Or what is the proper way to do this? Thanks Sverre


title: How do you deal with a "multiple choice" observation in Bayesian inference, when the choices are on a scale?

text: Suppose I have a questionnaire and I ask respondents how often they eat at McDonalds:   1. Never   2. Less than once a month   3. At least once a month but less than once a week   4. 1-3 times a week   5. More than 3 times a week I then correlate these answers with whether the respondents are wearing brown shoes.   1. Brown 65 -- not brown 38   2. Brown 32 -- not brown 62   3. Brown 17 -- not brown 53   4. Brown 10 -- not brown 48   5. Brown 9 -- not brown 6 _**The thing I can't get my head around is this:_** If a respondent picks #5, he (statistically) has a higher probability of wearing brown shoes than not. But, in a sense, his response subsumes responses 2-4, and if you accumulate their statistics (ie, "eats at McDonalds sometimes") he has a higher probability of _not_ wearing brown shoes. Now I realize that there are a bunch of caveats here -- sampling error in the stats, etc. But is there ever a valid argument for "rolling up" the stats (so that the values used in inference for #5 eg, would consist of the sums of the 2-5 values, or some other scheme), or is this concept just a product of my twisted mind? (Note that I'm not talking about "collapsing" the stats into fewer observations, which I assume would be perfectly valid, but rather adjusting the probabilities that are used in inference, based on the knowledge that the various possible mutually-exclusive observations are on a sliding scale.)


title: What does X represent in this linear trendline?

text: I'm trying to extend this graph via the provided trendline formula. I'm inserting the formula into the next row and trying to assign X to be the date value, but it doesn't appear to be working. Do I need to apply a function to the date value in order for it to work? Basically I'm trying to predict what June 17th's downloads will be based on the formula. Anyone know how do to this? ![enter image description here](http://i.stack.imgur.com/OCiv7.png)


title: any alternative to kaggle for data mining & analysis competitions and training?

text: Are there any alternative to www.kaggle.com to test abilities for data mining? Or how can I find my own projects to test myself on data mining projects? Kaggle is a data mining community that constitute some competitions and I am looking for some alternative web site or community to better my data mining and machine learning abilities. If there aren't any alternative, what are the common ways for the data mining experts to advance their skills and practice more. I want to be a member f such a community you advice or use the ways of practising to advance my ML and data mining skills.


title: Exercise in testing a difference of means

text: A random sample of 49 Kwakiutl native Americans reported that they went fishing an average of 68 days each year. Information from a non-native local population showed that they fished an average of 66 days per year (sd = 14). Do the Kwakiutls fish more than the non-Kwakiutls and if so is it a big difference?


title: pairwise t test or Tukey's test?

text: I have conducted an experiment with one group of subjects ($n=14$). I have two within-subjects factors and one dependent variable (called Error). The factor "Type" has 2 levels (Baseline, Operant); the factor "Event" has two levels (Action, Tone). The hypothesis is that there is a significant interaction Type*Event, specifically because:   1. $Error_{(Baseline-Action)}-Error_{(Operant-Action)}\leq 0$   2. $Error_{(Baseline-Tone)}-Error_{(Operant_Tone)}\geq0$ So, I first do a two-way repeated measures ANOVA. In R:               aov1 <- aov(Error~Type*Event+Error(subj/(Type*Event)), data=my_data)     summary(aov)      and the interaction is significant. Now I have to do a multiple comparison. I think it is possible to use both a pairwise t test:               # 'inter' is a column created by interaction(my_data$event, my_data$type)     pairwise.t.test(my_data$Error, my_data$inter, paired=T, p.adj="fdr")      and a Tukey test:               model <- lme(value~inter, random=~1|subj/(inter), data=my_data)     comp <- glht(model, linfct=mcp(te="Tukey"),type="none")     summary(comp)      The pairwise confirms my hypothesis 1 and 2, but the Tukey doesn't (it return a significant p-value only for the hypothesis 2). However, I don't know if I'm mistaken something with the commands for the Tukey. **My questions are:**   * **is this modus operandi correct in my case?**   * **on what basis I can choose which among pairwise t test and Tukey's test is the most appropriate? Are there other alternatives?** **EDIT.** What is the correct procedure for computing effect size from the output of glht (which give z-scores) and from a pairwise comparison t test (which doens't return the t values)?


title: Please help me understand white noise and MA(q)

text: I am reading the section about moving average models in Hyndman & Athanasopoulos _Forecasting: principles and practice_. I am trying to understand the MA(q) model in words. What is white noise? Is this a differenced series which is normally distributed with mean zero? Is it the difference between an observation and the mean of all observations? I do not know what the book is talking about when it says "white noise". I can understand what a differenced series is. I can understand what sum of square error means. But what is this "white noise" and where did it come from? What is an error term? What does it mean? Who made this up? Can I see an actual example that I can work out in Excel? When forecasting with an MA(q) model, do you add the moving average series to the mean to get a forecast? How does it actually work? An Excel document or an example involving actual numbers would really help. I am having a lot of difficulty understanding what is actually going on in the formula (reproduced below). Some examples with actual numbers would be great. $$ y_t=c+e_t+θ_1e_{t−1}+θ_2e_{t−2}+⋯+θ_qe_{t−q} $$


Who are some famous USSR/Russian Statistics/probability academic researchers? Are russian (current Russia and former USSR) strong in this field academically?


title: How to generate ROC Plot for semi-supervised algorithm?

text: By having a data-set 1000 (900 unlabeled, 100 labeled) record data-set for binary classification, I want to apply a semi supervised algorithm. The problem is that I don't know how to get values for the ROC plot? Actually the algorithm tries to find a threshold on a probability value (base on the 100 labeled records) in order to achieve the best F1 score for binary classification. What I'm thinking is to set different values to the threshold but I don't know that after selecting a threshold should I try it on the 100 records in training to get output point for ROC (which means evaluating the performance of the threshold on training set) or I have to test it on 900 records (test) to get the output point for ROC (which is the evaluation of the select threshold's performance on the test set)? Please help me!


title: Difference-in-Difference Regression with Only 4 Data Points

text: Is it possible to construct a difference-in-difference regression with only 4 data points, representing `Group 1`/`2` `pre`/`post` values? For example, Wooldridge describes a setup where there is a random sample of units (e.g. firms) in two periods of time under two different treatments. Say, firms in California in 2012 and 2013 and firms in Wisconsin in the same years. Further, there is a change whereby, say California in 2013 changed their tax policy. So as he describes, it sounds like we are not tracking the same units per state in the two years but a different sample per state per year (i.e. 4 distinct samples). **This leads to my question:** using the example above, can this same setup be used where you have a summary of the outcome variable for each state and year (not samples from each)? So, there are only 4 data points, say, total revenue of fast food restaurants in California in 2012 and 2013 (and the same for Wisconsin). Thus the regression uses 4 data points, but clearly a standard error on the interaction term cannot be estimated.               dat<-data.frame(sales=c(1234,2232,1530,2500),year=c(0,0,1,1),state=c(0,1,0,1))     summary(lm(sales~year*state,data=dat))


title: Feature selection on microarray for supervised classification

text: I'm studying the task of feature selection on biological microarray, thats is high dimensional dataset (thousands of features) with small number of data points (lees than one hundred). This feature selection serves to support supervised classification task. What research or algorithm you can recommend me? For now I learn about permutation-based algorithm, but I'm open to any advice. I need real research paper to understand how it works, I'm not looking for pre-packed programs. **EDIT1:** I have read 2-3 article dealing of feature selection based on correlation between feature. The central idea of this article is to remove the correlated feature because bring redundant information. For me this is not be a good idea because some classifiers could exploit the dependence between freature (Like SVM or Decsion Tree). What you think about? **EDIT2:** Is there a specific term to describe a High Dimensional Dataset with a few number of osservation like microarray? because I need it for improve my search on Google Scholar.


title: Two methods to calculate the confidence interval for number needed to treat yield different results

text: Assume that two drugs were tested. The risk of death for drug 1 is $p_1$ and the risk for drug 2 is $p_2$. We define:   1. Risk difference (RD) $RD=p_1-p_2$   2. The number needed to treat (NNT) $NNT=1/|RD|$ If we know the estimated RD as RD* and its standard error as se(RD*), what is the 95% CI for NNT? I can think of two methods for solving this problem. Which one is correct and why?   1. We first construct the 95% CI for RD, and then obtain 95% CI for NNT by inverting the CI for RD, that is:  step 1: 95% CI for RD: $RD^* \pm 1.96* se(RD^*)$   step 2: 95% CI for NNT: $1/(RD^* \pm 1.96* se(RD^*))$   result: $(7.5, 149)$   2. We first derive se(NNT*) from se(RD*） by the Delta method, and then calculate the 95% CI for NNT by:  step 1: $se(NNT^*)= se(RD^*)~|~d(NNT^*)/d(RD^*))$   step 2: $NNT^* \pm 1.96* se(NNT^*)$   result: $(1.38, 27.07)$ Obviously the two results are quite different. What is the problem in these two methods?


title: Scaled graph laplacian in presence of loops

text: I am interested in spectral clustering so I was looking the scikit-learn code for computing the Lapacian of a graph given its weighted adjacency matrix.               def _laplacian_dense(graph, normed=False, return_diag=False):     n_nodes = graph.shape[0]     lap = -np.asarray(graph)  # minus sign leads to a copy          # set diagonal to zero     lap.flat[::n_nodes + 1] = 0     w = -lap.sum(axis=0)     if normed:         w = np.sqrt(w)         w_zeros = (w == 0)         w[w_zeros] = 1         lap /= w         lap /= w[:, np.newaxis]         lap.flat[::n_nodes + 1] = (1 - w_zeros).astype(lap.dtype)     else:         lap.flat[::n_nodes + 1] = w.astype(lap.dtype)          if return_diag:         return lap, w     return lap      The code appears to zero out the diagonals unconditionally, and computes the degrees based for this modified adjacency matrix with 0 diagonals. This does not make a difference when computing the unscaled Laplacian but it does make a difference when computing the scaled/normed Laplacian. Ulrike von Luxburg's tutorial defines the degree of vertex $i$ given weighted adjacency $(W_{ij})$ simply as $d_i = \sum_{j=1}^{n} w_{ij}$ unlike the code above. Is this a bug? Or there is a reason for computing the Laplacian this way?


title: If the dependent variable is standardized by age and sex, does it still make sense to include these as controls in a multivariate regression?

text: Suppose I have a dependent variable that is a test score, $score^*_{i}$. It has been standardised for both age and gender. I want to measure the effect that a binary variable, say Adoption, has on the persons test score. Will it cause any problems if I include age and gender as controls in the regression of score on adoption (given that these were used to standardise the dependent variable)? That is, I am regressing (using standard OLS): $score^*_{i} = (score_{i} - \alpha_{1}AGE_{i} - \alpha_{2}Gender_{i}) = \beta_{0} + \beta_{1}Age_{i} + \beta_{2}Gender_{i} + \beta_{3}Adopted_{i} + \epsilon_{i}$   where $score_{i}$ is the raw score (i.e. unstandardized score), $\epsilon_{i}$ is the iid error term. I want to know if $\beta_{3}$ is affected by the inclusion of age and gender? My thinking is that this may bias the OLS estimates of the age and gender coefficients, but that it would not affect the the estimate of interest--the Adoption estimate. Is this correct? Or is it wrong altogether to include age and gender as controls in this setting? If so, could you please explain why? Thanks!


title: Imputation of missing parameter value by regressing the other parameters

text: I have a question on imputation of missing values. A support vector machine does not work with missing values, and in practice i thus hence na.omit the rows / cases with missing values on any of the parameters (e.g. X1 to X10). From a theoretical perspective, can one improve the overall classification performance if one imputes the missing values for the values of one missing parameter (lets say X1) through predicting the missing values of this parameter based on e.g. a multiple regression, SVM on the remaining parameters (e.g. X2-X10)? Asking the question, because the performance in practice is often disappointing, and hence wondering whether there is a theoretical reason why my approach per above might be nonsense / not useful. In a similar vein - as a decision tree does work with missing values, and hence should use the full information of the existing dataset, am i right to assume that in this case from a theoretical perspective there is definitely no value-add / benefit in imputing missing values by using the other input parameters?


title: Stat tools and methodologies for principal-agent problem

text: What are the most useful, popular or important statistical tools and/or methodologies for empirical study of 'agent-principal problem'? A simple example will be very useful. References or book recommendations are highly appreciated. Thank you very much!!


title: Rescale vector logarithmically in matlab

text: I'm quite weak in math so I hope you don't mind that I'll try to fully describe how I think when solving this (I think easy) problem. Please tell me if I made some totally invalid assumptions. I have some data - a slope of a linear regression curve of a price of a particular commodity in stock exchange. This is the data plot.               data = importdata('E:\data.txt');     plot(data);      ![Unscaled data plot](http://i.stack.imgur.com/bUgHx.png) These data are used as the training data for a neural network. And since the NN uses tanh as the activation function, I need to rescale these data to the range -1 to 1. I've managed to do that.               dmin = min(data);     dmax = max(data);     dataScaled = ((data-dmin)./(dmax-dmin) - 0.5 ) *2;      And the plot of the scaled data:               plot(dataScaled);      ![Wrongly scaled dataplot](http://i.stack.imgur.com/dYS3l.png) This plot is wrong though, the linear regression slope by nature oscilates around 0, scaling by the maximum and minimum shifted the oscilation center down because the maximum and minimum values were different. So I looked at the maximum and minimum and saw this.               dmax = 3.7348     dmin = -1.8465      I told myself that I cannot rescale by asymetric minimum and maximum...so I thought I'll take the one furthest away from zero and use that with opposite sign for the second value. Since |dmax| > |dmin| then               dmin = -dmax;      Now that both values are the same distance from zero, I'm going to rescale by that and hopefuly the rescaled values will oscillate around zero.               dataScaled = ((data-dmin)./(dmax-dmin) - 0.5 ) *2;     plot(dataScaled);      ![Data scaled the right way](http://i.stack.imgur.com/YkwiL.png) And now the problem. The standard deviation of the scaled data is 0.0553 - that's very wrong (I discovered that in my previous question here). My neural network is not learning anything becuase the standard values are pushed near zero. Well, I guess the NN is learning something but on a scale so small that it makes working with the neural network difficult (at least difficult for me). So I thought I have to ...somehow rescale it, logarithmically, in a way that the logarithm is going from 0 to positive values the same way it goes from 0 to negative values. The standard deviation is 0.0553 ... so most values oscillate between -0.0553 to 0.0553. I'd like to convert the vector in a way that the oscillation from -0.0553 to 0.0553 _stretches_ to the range, say, -0.8 to 0.8 - I'd like to test more ranges to see which one I like the most - and values above or below the standard deviation are _compressed_. How do I do that in matlab? Does my question make sense? Just tell me and I'll be hopefully able to rephrase it. **TL;DR** Just read from "And now the problem".


title: Binary regression with error-in-variables

text: I try to predict the binary, dependent variable $Y$ using a single predictor $X$. $X$ contains some amount of error, due to measurement problems. I would like to come up with a model for $X$, e.g. assuming normality for errors: $X \sim N(X_{true}, \sigma)$, where $X_{true}$ is the true, but unobserved value. My interest lies on $\sigma$ here. Theory tells me that there are subgroup differences in $\sigma$. At best I would like to write an own regression line for $\sigma$ using other covariates. My question is this: Is this at all possible and what are the identification issues? If anyone could point me to relevant literature that could help me solve this, I would very much appreciate this. Thanks.