title: what does an actual vs fitted graph tell us?

text: I am just wondering what we can infer from a graph with x-axis as the actual and y axis as the predicted data? ![signups](http://i.stack.imgur.com/dwlOW.png)


title: Test of significance between gain scores of intervention group and control group

text: I am performing a significance test on the gain score between a Pre-Test (say Y1) and a Post-Test (say Y2) for a Treatment Group and a Control Group. Let me give a brief background on the experiment. My client is working with a large group of schools on improving the learning levels of kids in Math. Out of this large group 653 students were tested in the Pre-Test. There was also a control group which had 213 students who were also tested. After this the intervention group received special coaching (intervention) along with regular classes. And the control group just had regular classes. About 6 months later, a Post-Test was conducted. But this time 739 students (different set of students from the intervention group, may have some overlap) from the intervention group and 210 students (pretty much the same group as the Pre- Test) from the control group were tested. The problem I am facing is how to do a significance test on the difference between the Gain Scores. Here is the method which I am planning to employ. Since I am sampling large enough samples from the Intervention Group and the Control Group for the Pre-Test, and I have point estimates for the Mean and SD, the CLT gives me a normal distribution for the Pre-Test for both Intervention and Control. The same argument applies to the Post-Test. Now I assume that the Pre and Post Test performances are independent. Now I can combine Pre and Post test normals to create a new normal distribution for Y2-Y1. The same argument again applies to control also. Now I want to apply Welch's t-test to check for statistical significance between the groups. But I don't have the degrees of freedom and number of observations for both Intervention and Control. I'll use the Welch-Satterthwaite Equation to get the Degrees of freedom for each group and then add one to get the number of observations. After this I can apply Welch's t-Test. When I am doing this I get extremely high t-statistic values, So I am not really sure whether I am doing it the right way. Could someone please tell me if there is a flaw in my logic.


title: How do I give a larger weight to a smaller population?

text: I have a set of cities ranked 1-10 based on several criteria. I also have the percentage of the total population living in these cities. I want the cities with smaller populations to be able to compete with larger cities, so I want to give the population percentages a weight that will compensate for smaller population size. How do I do this?


title: Data Analysis using R

text: In preparation for my end of semester R exam, I've been working on the following CSV file which was modelled on R, which I'm currently having problems figuring out. The excercise is from the textbook'Applied Statistics: Principles and Examples' by D. David Roxbee Cox. ![enter image description here](http://i.stack.imgur.com/6wvqq.png) _The data is taken from a real study into Light Water Reactor (LWR) power plants constructed in the USA in the 1970s. The aim of the study was to develop a model that could be used to model the capital cost involved in the construction of LWR power plants. The data set includes 32 LWR power plants and contains the following variables:_ The titles PR, NE, CT, BW and PT have outputs that consist of 1 and 0 where 1=Yes and 0=No. The data set includes 32 LWR power plants and contains the following variables: ![enter image description here](http://i.stack.imgur.com/095E0.png) ![enter image description here](http://i.stack.imgur.com/6FatZ.png) Now that I have introduced the domain of the question - here is my problem: I've been asked to 'add suitable terms' to that model in order to investigate whether there is good reason to believe that the impact of a power plant's capacity on construction costs might not be the same for plants that did have a LWR on the same site previously and those that did not. The summary results of my final model in R is shown below: ![enter image description here](http://i.stack.imgur.com/fzSel.png) Therefore, the equation of my final model is:               C = -541.82231+ 26.11545T1+ 5.22321T2+ 0.41890S- 94.33786PR+108.34622NE- 132.58098PT      I tried fitting a model assessing just the S (power capacity) and PR (whether LWR existed on the same site previously) terms against C (construction cost) but the wording of the question states 'add suitable terms' to that model in order to investigate...' causes me to believe what I wanted to do is wrong because I would have been taking out terms. Does anyone know what suitable terms I should add and how I test the correlation?


title: Relationships between the corr(x,y) and corr(x,y) on a subset

text: Suppose two sets of data $\mathbf{X}$ and $\mathbf{Y}$ have a correlation of $\rho_0$. Can anything be said about the correlation between $\mathbf{X}_1$ and $\mathbf{Y}_1$ and $\mathbf{X}_2$ and $\mathbf{Y}_2$ where $\mathbf{X}_1\cup \mathbf{X}_2=\mathbf{X}$ and $\mathbf{Y}_1\cup\mathbf{Y}_2 = \mathbf{Y}$ and $0<|\mathbf{X}_1| = |\mathbf{Y}_1|<|\mathbf{X}|$? In a practical sense, would there be any relationship (or constraint) between $r$ and $r_1$ and $r_2$ in the below example?               df.0 # some data frame containing x and y     df.1 <- df.0[df.0$x>5, ]     df.2 <- df.0[df.0$x<=5, ]     r <- cor(df.0$x, df.0$y)     r1<- cor(df.1$x, df.1$y)     r2<- cor(df.2$x, df.2$y)      I'm thinking something like $r \leq \max\{|r_1|,|r_2|\}$. Quite happy if there is none!


title: Algorithms for regression analysis which can handle large scale datasets

text: I am a CS undergraduate student and for my final project i developed a regression algorithm that is suited for large-scale datasets (i wouldn't say 'Big Data', but still large scale). For the final results, i'm looking for regression algorithms that are known to successfully work on such large datasets, so i can use them as a baseline for my algorithm. I dont have limits to some particular proerties of the dataset, so the algorithm should be general. Does anybody have an idea for regression models i can use? EDIT: To clarify: the best thing for me is to find a paper which shows one (or more) regression models, and their performances on large-scale datasets. If anybody could also point me out to a known public dataset (that is big of course) it will be very helpfull. Thanks for helping, Nir


title: Is ANCOVA proper for my experiment?

text: I've got a serious problem when trying to analyse the data of my psychological experiment. Here are the details of the experiment: In each trial, subjects were presented with lateralized colored circles. They had to press a right button when the circle was red, or a left button when the circle was blue. The compatibility condition determines if the location of the circle and the side of the response correspond (compatible) or not (incompatible). The saturation of the colors were also manipulated (6 saturation levels: 15, 25, 35, 45, 60, 80%). So I used a 2 (compatibility condition: compatible, incompatible) x 6 (color saturation level: 15, 25, 35, 45, 60, 80%) within-subjects design. 13 subjects participated in the experiment. For each subject, there were 224 observations per experimental condition. I would like to analyse the relation between the mean and variance of each experimental condition for each individual. Classicaly, the relation is linear, leading to very high correlations for each subject (r > .85). However, in my experiment, the r values range from 0.15 to 0.95 with a mean of 0.56. The majority of subjects shows an r < .75 (11 out of 13 subjects), which is far from the traditional finding (> .85). So I decided to switch to an other analysis. I think that the data is better described by two regression lines: one for each compatibility condition (but this is just a graphical feeling). My guess is to perform an ANCOVA for each subject with std as the dependant variable, mean as covariate and compatibility as factor. But I don't know if it is the right analysis, and if it could be considered as a within- subjects ANCOVA. Any help would be appreciated!


title: Expected number of ratio of girls vs boys birth

text: I have came across a question in job interview aptitude test for critical thinking. It is goes something like this: > The Zorganian Republic has some very strange customs. Couples only wish to > have female children as only females can inherit the family's wealth, so if > they have a male child they keep having more children until they have a > girl. If they have a girl, they stop having children. What is the ratio of > girls to boys in Zorgania? I don't agree with the model answer given by the question writer, which is about 1:1. The justification was any birth will always have a 50% chance of being male or female. Can you convince me with a more mathematical vigorous answer of $\text{E}[G]:\text{E}[B]$ if $G$ is the number of girls and B is the number of boys in the country?


title: Predicting temperature time series with Holt-Winters

text: I am trying to write a prediction algorithm for a set of temperature data. I settled on Holt-Winters since it seemed to be a simple time series prediction algorithm and I can easily code it up in python to understand what is going on with it. When I am plotting the smoothing function as it learns this is how it looks. As you can see, it follows the original curve pretty well. ![enter image description here](http://i.stack.imgur.com/bZtGz.png) But when I try to plot a future curve for one year (365 days) -- it really falls down and dies. ![enter image description here](http://i.stack.imgur.com/7mwJP.png) And to me intuitively it makes sense why it dies like that. Because if you see the last prediction equation of Holt-Winters it really only makes use of the very very last point in both the curve smoothing and the trend smoothing. And we know that exponential smoothing has a very short memory because of the whole exponential thing. So I am wondering how does one actually go about using Holt-Winters for prediction (specifically for 365 day seasonal data like this). If you know of any other methods which I can look at for this domain (temperature prediction) please let me know. I come from python, so it would be very useful if you can point me to resources such as libraries etc to get the job done.


title: What do these chi-square results mean?

text: I have some items described by 43 categories like this:               Dataset Item Category1 Category2...               Category43     D1      Item1 1        0              0    ...    1 ...     D1      Item3 1        0              0    ...    1 ...     D2      Item4 1        0              0    ...    1 ...     ..      What I did was to create a frequency table like this               Dataset Category1 Category2...               Category43     D1        617       388    ...                 827     D2        1234      7272   ...                 1237      I am testing to see if there is a relationship between the dataset type and the category frequency counts. I have the following data as the output of `dput`:               structure(list(data.OldFrequency = c(617L, 388L, 6L, 9L, 1344L,      857L, 30L, 63L, 60L, 22L, 23L, 107L, 9L, 16L, 9L, 10L, 14L, 28L,      9L, 174L, 245L, 103L, 4096L, 121L, 6L, 48L, 189L, 33L, 1426L,      64L, 16L, 135L, 77L, 26L, 110L, 44L, 75L, 1610L, 1022L, 38L,      1578L, 242L, 67L), data.NewFrequency = c(1220L, 959L, 307L, 29L,      5093L, 771L, 65L, 125L, 120L, 41L, 187L, 203L, 11L, 87L, 20L,      159L, 45L, 68L, 60L, 11L, 644L, 51L, 7053L, 159L, 6L, 162L, 208L,      52L, 3277L, 27L, 594L, 79L, 95L, 119L, 96L, 84L, 180L, 2991L,      2227L, 34L, 2249L, 37L, 29L)), .Names = c("data.OldFrequency",      "data.NewFrequency"), row.names = c(NA, -43L), class = "data.frame")      Running `chisq.test` using this gives me the following:               Pearson's Chi-squared test          data:  d      X-squared = 2551.405, df = 42, p-value < 2.2e-16          Warning message:     In chisq.test(d) : Chi-squared approximation may be incorrect      I am confused on what null hypothesis this is testing and what the implications of this are. Can someone please help me understand how to interpret this? I am not a statistician and would love if someone could explain this in simple words. And how would I fix the warning message?


title: Selecting features using Adaboost

text: How does Adaboost select best features from the sample data (or a unit feature vector)? It would be nice if someone can explain if the above statement is true or not. I've seen the term _features_ and _classifiers_ being used in place of each other in documentation that's available of Adaboost (that's what I thought). For example: www.cs.nyu.edu/~eugenew/publications/viola-facedet04-talk.pdf Please explain!


title: random forest oob error increase as more tree build

text: I am implementing random forest in my code, I find that when I build one tree, the oob error(mean square error since I do regression) is close to zero, while more tree is build, the oob error stablized, this is counter-intuitive, since the textbook teaches me that the oob error should decline as more trees builds, I compare my implementation with R, the my oob is a slightly less than that of R, when mtree=1000, but when mtree=1, my oob is close to zero, while R is quite big. Simply put, my oob increase as more trees build and stabilized. R and text book shows that oob decrease as more trees build and is stabilized. So, is there anything wrong with my implementation? how should I tune my algorithm?


title: SVM importance of predictor variables

text: I am building a model in R using support vector machine (SVM) with KBF kernel. The model seems to work quite well. I would like to assess the relative importance of predictor variables. Can anyone tell me how this can done. Does there exist an R script for that?


title: Specific contrasts in mixed model with interaction

text: I have a dataset consisting of two groups tested across three days. Therefore I run a linear mixed model as follows: `model1 = lme(response ~ day*group, random=~1|ID, data)`, that reveals a significant interaction. I would like to know if there is a difference between the groups in how they differ over D1-D2, D2-D3, and D1-D3, and potentially how groups 1 and 2 differ on each of the days. As I understand that one can only perform (n-1) contrasts, I initially only ran contrasts for the first two comparisons, `contrasts(data$day) = t(ginv(c(-1,1,0),c(0,-1,1))`(group did not need contrasts as there were only two), which showed output for D1D2GR1 vs D1D2GR2 and D2D3GR1 vs D2D3GR2. I then read about the option to create a new variable that combines the two other variables (called 'newvar' with D1G1, D2G1, D3G1, D1G2, D2G2, D3G2 as levels) to be able to create more custom contrasts:               D1D2 = c(-1,1,0,-1,1,0)     D2D3 = c(0,-1,1,0,-1,1)     D1GR = c(-1,0,0,1,0,0)     D2GR = c(0,-1,0,0,1,0)     D3GR = c(0,0,-1,0,0,1)     contrastmatrix = cbind(D1D2,D2D3,D1GR,D2GR,D3GR)     contrasts(data$newvar) = t(ginv(contrastmatrix))     model2 = lme(response ~ newvar, random=~1|ID,data)      **Question1:**   Are both these ways of getting contrasts appropriate and performed in the right way? **Question2:**   Is there a way to get the three contrasts 1) D1D2GR1 vs D1D2GR2; 2) D2D3GR1 vs D2D3GR2; and 3) D1D3GR1 vs D1D3GR2 directly from the new variable? As I have 6 levels, conducting these three contrasts only shouldn't it be possible to do this? **Question3:**   If I would like to run post-hoc tests similar to the contrasts, how could I do that for the first model? And is running it for the model with the 'newvar' as follows the right way? In other words, is it right for the contrasts to do that on the generalised inverse matrix and for the post-hoc tests to do it on the non-generalised inverse matrix?               glht(model2, linfct = mcp(newvar = t(contrastmatrix)))      **Question4:**   Finally, I would like to get effect sizes for the observed effects. Is it right that for the contrasts I can get this by `r = r<-sqrt(t^2/(t^2 + df))` with the t-values and df from the output? How can I calculate this from the z-values from the post hoc output? The z-values look almost exactly the same as the t-values.


title: Sufficient Statistic

text: I have a question: > Does a sufficient statistic have to be one to one? For example, can $T(x) = > x^2$ or $T(x) = |x|$ be sufficient statistics? I know that one to one > functions of sufficient statistics are sufficient.


title: t-test on seasonal effects on water quality

text: I am confused in the use of t-test. I monitored physical parameters of a river water two time (from the same site)- during winter and summer. Now I want to see the seasonal difference on physical properties, say for example water temperature. Should I use Independent Sample t-test; or paired sample t-test??


title: How should very large sample sizes be used?

text: At what point can very large samples, as a proportion of the population, be treated as a census? For example, if your sample contains 90% of the population units, can one dispense with inferential statistical inferences? What about 80%? And so on . . ..


title: How to improve location and scatter estimation conditioning on higher order statistics?

text: Using sample moments, how can the mean and variance estimators be improved if e.g. skewness and kurtosis are known exactly? And what about using estimates for these instead, which should imho be of no help? In general then I would like to know if there are analogous approaches for other kinds of location and scatter estimators, such as L-moments.


title: How to find the function of a skewed distribution curve by calculating k=P for 4 points along the curve

text: Say the minimum of k is 0 and the maximum of k is 4. k is being derived from the probability mass function, though I'm not using the formula because I do not have p of the distribution. If I have k(0) = 0.1852 k(1) = 0.4074 k(3) = 0.0926 k(4) = 0.0092 Is there a function that I can use to calculate the curve? For a little more clarification, what I am really looking for is to calculate the function of a skewed bell curve by using k(0),k(1),k(n), and k(n-1), where n could be any number up to 50. I am calculating these k's manually because the set includes multiple probabilities (ex: a dice roll and a coin flip). As stated earlier, I do not know p for the distribution.


title: Intraclass correlation in the context of linear mixed-effects model

text: Suppose that one investigates the two sources of variability with data $y_{ij}$ acquired from $j$th subject under $i$th session ($i=1,2; j=1,2..., n$). A linear mixed-effects model can be formulated as follows, $$y_{ij} = \alpha_0 + b_i + c_j + \epsilon_{ij}$$ where $\alpha_0$ is a constant, $b_i$ and $c_j$ are the random effects for $i$th session and $j$th subject respectively, and $\epsilon_{ij}$ is a residual term. With assumptions of $b_i ~\sim N(0, \tau_1^2)$, $c_j ~\sim N(0, \tau_2^2)$, and $\epsilon_{ij} ~\sim N(0, \sigma^2)$, the intraclass correlation (ICC) values for sessions and subjects can be defined respectively as $$ICC_{session} = \frac{\tau_1^2}{\tau_1^2+\tau_2^2+\sigma^2}, ICC_{subject} = \frac{\tau_2^2}{\tau_1^2+\tau_2^2+\sigma^2}.$$ And the above ICC values can be obtained through `lmer` in R package `lme4`. I have two questions:   1. Is there a way to test the significance of the above ICC values in the context of LME model, similar to the $F$-stat for ICC defined under random-effects ANOVA?   2. Intuitively the bigger the variability between the two sessions (more different between the two sessions), the higher the session ICC. But how is this intuition consistent with the notion that ICC measures reliability, reproducibility, or consistency? Is the session reliability (or reproducibility, consistency) about any two subjects' responses coming from the same session, not about the difference between the two sessions?


title: What are easy to interpret, goodness of fit measures for linear mixed effects models?

text: I am currently using the R package lme4. I am using a linear mixed effects models with random effects:               library(lme4)     mod1 <- lmer(r1 ~ (1 | site), data = sample_set) #Only random effects     mod2 <- lmer(r1 ~ p1 + (1 | site), data = sample_set) #One fixed effect +                  # random effects     mod3 <- lmer(r1 ~ p1 + p2 + (1 | site), data = sample_set) #Two fixed effects +                  # random effects      To compare models, I am using the `anova` function and looking at differences in AIC relative to the lowest AIC model:               anova(mod1, mod2, mod3)      The above is fine for comparing models. However, I also need some simple way to interpret goodness of fit measures for each model. Does anyone have experience with such measures? I have done some research, and there are journal papers on R squared for the fixed effects of mixed effects models:   * Cheng, J., Edwards, L. J., Maldonado-Molina, M. M., Komro, K. A., & Muller, K. E. (2010). Real longitudinal data analysis for real people: Building a good enough mixed model. Statistics in Medicine, 29(4), 504-520. doi: 10.1002/sim.3775    * Edwards, L. J., Muller, K. E., Wolfinger, R. D., Qaqish, B. F., & Schabenberger, O. (2008). An R2 statistic for fixed effects in the linear mixed model. Statistics in Medicine, 27(29), 6137-6157. doi: 10.1002/sim.3429  It seems however, that there is some criticism surrounding the use of measures such as those proposed in the above papers. Could someone please suggest a few easy to interpret, goodness of fit measures that could apply to my models?


title: Reject null, if and only if q-value is no more than the FDR threshold?

text: In multiple comparisons, suppose we have calculated $q$-values. Is it a $q$-value for each test, or a $q$-value for all the multiple comparisons? Given a false discovery rate (FDR) threshold on FDR, are $q$-values used to determine whether to accept or reject each null? For an individual test, is it rejected **if and only if** its $q$-value is no less than the FDR threshold? Why? Thanks!


title: Point hypothesis and equivalence hypothesis at once

text: Sometimes researchers (especially in collaborations) have two opposite but sound theories: There is a difference between two groups or there is only negligible difference. Now they ask their statistician. How should he approach the situation? If he does a point hypothesis test, he favours the difference party since only their theory can be confirmed in case of rejected point hypothesis whereas not rejecting the point hypothesis teaches nothing. If he does only equivalence testing, he favours the negligible difference party for the analogous reason. Should he do both, of course with multiplicity correction? So a TOST for equivalence _and_ a point hypothesis test? Or a TOST and a respective relevance test? This procedure would have three outcomes:   1. Equivalence of both parameters up to negligible differences.   2. Large enough differences.   3. Nothing to learn since both hypotheses have not been rejected. Is such an "et-et"-approach reasonable? Why do we hardly see such "et- et"-analyses in publications? This is a general question for reasoning. So multiple answers are encouraged and I do not restrict this question to particular models.


title: How to format multi-row time series data for LIBSVM regression

text: I would have expected this to be covered in detail by the LIBSVM tutorial but after hours of wasting time googling for answers I've had to throw in the towel. What I am trying to do is rather trivial really. Right now I have a financial data set in CSV format:               time,indicator1,indicator2,indicator3,roc     20140508190000,0,-2,-1,-0.00361310835712725     20140508200000,1,-2,0,0.0216786501427154     20140508210000,-1,-2,1,-0.0252917584998426     20140508220000,0,-2,0,-0.0108397167220604     20140508230000,-1,-2,0,-0.0505725535527163     201405090,0,-2,0,0.0144560896277541     2014050910000,0,-2,0,-0.00361363061468619     2014050920000,1,-2,0,0.0542122953485871     2014050930000,-1,-2,1,-0.0975645009756423     2014050940000,0,-2,1,-0.0144550448106373     2014050950000,-1,-2,1,-0.198670712324811     ...      So we have the date/time column and then four attributes per row. Later on I am going to want to predict the last column (roc). This data is obviously not normalized and I intend to use phraug to convert my data into a LIBSVM friendly format. But not so fast... Here's the problem. I have not been able to find any examples that explain how to format multi-row input data. In other words, since I plan to use SVR I somehow need to feed my SVM a sliding window of data sets. For instance, let's say I want to give the engine five rows in order to have it predict the output value of the sixth:               20140508190000,0,-2,-1,-0.00361310835712725     20140508200000,1,-2,0,0.0216786501427154     20140508210000,-1,-2,1,-0.0252917584998426     20140508220000,0,-2,0,-0.0108397167220604     20140508230000,-1,-2,0,-0.0505725535527163      The ideal here would be the roc value of the next row: 0.0144560896277541 Then we slide the window down by one row:               20140508200000,1,-2,0,0.0216786501427154     20140508210000,-1,-2,1,-0.0252917584998426     20140508220000,0,-2,0,-0.0108397167220604     20140508230000,-1,-2,0,-0.0505725535527163     201405090,0,-2,0,0.0144560896277541      And the ideal here is: -0.00361363061468619 Etc. - in essence I need to feed LIBSVM _sets of data_ , not just atomic rows of data. Either I can just keep my original format (normalized) and configure LIBSVM's training/prediction routine to do that somehow, _OR I need to convert my five rows into one row_ , and I have no clue how to do that. Actually I really hope that this is done via some kernel hyperparam and not via the CSV. Otherwise it would require a different CSV input file for differing ranges of data sets. And that would be very complicated and error prone. I hope all this made sense. Can someone PLEASE point me in the right direction? Any input would be much appreciated. **UPDATE:** I think I'm starting to see the light! Just to clarify - I have no problem with HOW to format my data _once_ I fully understand the needed format. Let's assume I use regression and a window of 3 days of data with three input values per day. Per my understanding this involves creating a multi-dimensional vector, thus this would be the proper input matrix for a 3x3 vector: response_variable 1:day1-1 2:day1-2 3:day1-3 4:day2-1 5:day2-2 6:day2-3 7:day3-1 8:day3-2 9:day3-3 For SVC my response variable would most likely be an int representing a category, e.g. up/down (-1/1) and for SVR something like a normalized ROC (e.g. 0.014). Can someone please confirm this? Thanks in advance! :-)


title: GEE in SPSS: the case of ever-decreasing QICc values

text: I've been using GEEs in SPSS 22 to analyze my dissertation data and have discovered an interesting problem: when trying to figure out which subset of model factors have the lowest QICc, and therefore best fit my data, the QICc just keeps decreasing as I eliminate factors. Now, I know QICc penalizes for more factors in the model (favoring simplicity), but this is getting a little silly. Just for fun I went through and started eliminating even the significant factors and still the QICc decreased... More background: I'm modeling count data using a negative binomial distribution with a log link and exchangeable covariance matrix (selected using QIC values). I have three factors (let's call them A, B, and C) and two covariates. The factors are all categorical with two levels each. I want to test for the effect of the factors and possible interactions between them, so I started with the full factorial model. The 3-way interaction and two of the three 2-way interactions are not significant, so I decided to remove them from the model and see what happened to the QICc. Please feel free to correct me if this is the incorrect way of doing it... I'm used to AIC and AICc! Here's a short summary of what the pattern in QIC and QICc values I get. The same pattern happens when I use the full log quasi-likelihood function and the kernel log quasi-likelihood function, the raw numbers are just different. So, here's how it goes:               Model Factors = A, B, C, A*B, A*C, B*C, A*B*C   QIC = 0.626   QICc = 16.597       (Where A*B*C, A*B, and B*C are not significant and that doesn't change as I delete factors.)               A, B, C, A*B, A*C, B*C   QIC = 0.624   QICc = 14.597       A, B, C, A*B, A*C        QIC = 0.623   QICc = 12.597       A, B, C, A*B             QIC = 0.624   QICc = 10.607       A, B, C                  QIC = 0.635   QICc = 8.614        ...and so on. QICc keeps dropping by about 2 for each factor or interaction removed. So what is the scoop here? Is it that my model factors are just crap-tastic at explaining my data? But then why the significance? Should I just go with the full model since I really am interested in testing all the possible interactions and GEE really isn't meant for the same type of "model building" I'm used to with other methods? Any help would really be appreciated! Thanks in advance!


title: Predicting user selections based on similar user

text: Lets say you give a set of users a set of polls, or give them a choice of foods to eat, or let them listen to a group of songs (guess like pandora). So looking at the choices that all the users make you can match user 1 with user 5 based on the the 1st 10 choices they make. So you taking user 1 you present him with the next poll (or whatever) and he chooses, then I can take user 5 and present him with the same choices. I think I could then say with a certain probability that user 5 will select the same one that user 1 selected based on prior history. What is this statistical/mathematical process called? How would I code it in general? (I know this isn't the right place for it, but I suppose if I know what the process is called I could find it). Well now I see there is the tag called predictive models, but I am not sure if that is the right term and what specific form I am trying to describe?


title: How to avoid the problem of two-way causality?

text: I am studying the effect of social capital on households' income. I am doing multiple regression to estimate this effect. For this, I have households' income as dependent variable and social capital as explanatory variable. I have also included human capital and some other household characteristics as control variables for multiple regression. Social capital and household income have two-way causality as it is argued that social capital is endogenous variable. This is proven theoretically and empirically in various research studies. To avoid this problem and to get unbiased and consistent estimates. I intend to use two stage linear regression (2SLS) or instrumental variable technique. Can someone suggest which method is preferable in my case and how to perform my analysis exactly in both cases if:   1. I have valid instrument, or    2. I don't have valid instrument for the endogenous variable?


title: A research containing the statistical difference between two samples

text: I have conducted a research in which 120 people participated. For my survey I created two visual images (casual and clown).   * Sample 1 (n=60) conducted the interview with the casual person    * Sample 2 (n=60) conducted the interview with the clown The first part of the survey (which does not include data gathering) the participant had to answer multiple questions based on feedback of either the casual person or the clown. At the end of this narrative the participant had to fill in 8 questions in which the participant had to define the degree of trust it had with either the clown or the casual person. This degree of trust is based on a likert-scale (discrete values) ranging from 0(disagree) to 7 (agree). See example: http://imgur.com/K22QIkh.jpg So the goal of my research is to see if there is a significant difference in trust by:   * Comparing the results of Question 1 between sample 1 with sample 2 (Compare 60/60)   * Comparing the results of Question 2 between sample 1 with sample 2 (Compare 60/60) And so on for all the eight questions. I used a one sample t-test for the data analysis. In SPSS i put all 120 results under each other (0 to 60=sample 1, 61 to 120=sample2) and defined that 0 to 60 is group 1 and 61 to 120 is group 2. Thus, by doing this I could compare the results for all 8 questions between sample 1 and 2. This eventually gave me the following results: See example: http://imgur.com/2gMkrfB.jpg My main question is: With an alpha of 0.05 all the results of the questions seem to be significant. However question 6 (vraag 6) both have the exact same mean, how can this difference be significant? Am I doing something wrong? Any advice is welcome :)


title: Backward Stepwise Selection

text: I need to implement a backward stepwise regression. I read the chapter from "The Elements of Statistical Learning" however the explanation is poor here: > Backward-stepwise selection starts with the full model, and sequentially > deletes the predictor that has the least impact on the fit. The candidate > for dropping is the variable with the smallest Z-score **This is taken from > chapter 3.3.2 pg:59** I need to know what is the "full model" mentioned here? What will I delete at each level acording to the Z-score? I need to know how this algorithm works step by step. **EDIT:** I try to implement the code on matlab. I don't know I'm totally wrong or not :/ Comment on this function please               function [] = BackwardStepWise(X,y,N,p)  % X is Nxp matrice y is Nx1 matrice          X = [ones(N,1) X];          % vector for holding column numbers      v= 1:p;          for k=p-1:-1:1                  %create matrix with the selected columns             T= [ones(N,1)];             for j=1:k                 %add column to the matrix                 T = [T X(:,v(1,j)+1)];             endfor             %evaluate beta             beta = inv(T' * T) * T' * y;              %calculate Z-scores for each column         sigmahat_sq = (y-T*beta)'*(y-T*beta)/(N-p-1);         TT = inv(T'*T);             Zmin = 100000;             Zminindex = -1;         for i=1:size(beta,2)             z(i) = beta(i)/sqrt(sigmahat_sq*TT(i,i));                     if(z(i) < Zmin)                 Zmin = z(i);                 Zminindex = i;             endif         endfor              %drop column which has the smallest Z score (edit v vector)         v2 = [];         for i=1:size(v,2)             if(i == Zminindex)                 continue;             else             v2 = [v2 v(i)];              endif                endfor           v = v2;          endfor          endfunction


title: 2 standard deviation rule

text: I know that 95% of the observations under a normal distribution fall approximately under 2 standard deviations from the mean. Does this change when the distribution has fat tails? What is the name of the theorem that guarantees this? Thanks beforehand.


title: Is p=5.0% significant?

text: Today I was asked, whether a p-value of 0.05 (exactly) is considered significant (given alpha=5%) or not. I did not know the answer and Google turned up both answers: (a) the result is significant if p is less than 5% and (b) if p is less than 5% or equal to 5%. Of course, none of these websites cited anyone. Why should one - it's common knowledge and the 5% are arbitrary, anyway. But that won't help me to tell my students something to remember. So, here are my desperate questions on testing hypotheses: If the p-value is exactly the alpha - do I consider the result significant or not? And what is the authorative citation in this case? Thank you very much


title: Why signs of coefficients change when doing multivariate vs. univariate logit regression?

text: Excuse my dumb question, but I did an univariate logistic regression where the sign of the coefficient of my variable was negative (and it was significant). Once I have input it into a multivariate regression, the sign of the coefficient of the same variable became positive (still significant). Why does this occur? Is it because of outliers or collinearity between variables (even though as I said the variables are significant), or is it something more subtle?


title: Repeating ANOVA analysis after dropping one factor: different means

text: I had to repeat the analyses after dropping out one of the between-subject factors (order of the questionnaires) and the new means I got for other effects were different from those I had before. Can it be like that? Is that because the degrees of freedom changed?


title: Recommendations for visualization type when data has an extremely wide variance

text: I'm working on a web app, and I'm creating some data viz tools for it. For one particular series, I've got an extremely wide variance in data values (0 to millions). We're using a column chart to view the data now, which of course results in some columns that are a pixel high or smaller. We already have some ways to slice the data that helps a bit, but I was wondering if there were different kinds of visualizations out there in common use that deal with this type of situation better. And if so, if there were JS libraries that help implement them.


title: How to quasi match two vectors of strings (in R)?

text: I am not sure how this should be termed, so please correct me if you know a better term. I've got two lists. One of 55 items (e.g: a vector of strings), the other of 92. The item names are similar but not identical. I wish to find the best candidate* _s_ * in the 92 list to the items in the 55 list (I will then go through it and pick the correct fitting). How can it be done? Ideas I had where to:   1. See all the ones that match (using something list ?match)   2. Try a distance matrix between the strings vectors, but I am not sure how to best define it (number of identical letters, what about order of strings?) So what package/functions/field-of-research deals with such a task, and how? Update: Here is an example of the vectors I wish to match               vec55 <- c("Aeropyrum pernix", "Archaeoglobus fulgidus", "Candidatus_Korarchaeum_cryptofilum",      "Candidatus_Methanoregula_boonei_6A8", "Cenarchaeum_symbiosum",      "Desulfurococcus_kamchatkensis", "Ferroplasma acidarmanus", "Haloarcula_marismortui_ATCC_43049",      "Halobacterium sp.", "Halobacterium_salinarum_R1", "Haloferax volcanii",      "Haloquadratum_walsbyi", "Hyperthermus_butylicus", "Ignicoccus_hospitalis_KIN4",      "Metallosphaera_sedula_DSM_5348", "Methanobacterium thermautotrophicus",      "Methanobrevibacter_smithii_ATCC_35061", "Methanococcoides_burtonii_DSM_6242"     )     vec91 <- c("Acidilobus saccharovorans 345-15", "Aciduliprofundum boonei T469",      "Aeropyrum pernix K1", "Archaeoglobus fulgidus DSM 4304", "Archaeoglobus profundus DSM 5631",      "Caldivirga maquilingensis IC-167", "Candidatus Korarchaeum cryptofilum OPF8",      "Candidatus Methanoregula boonei 6A8", "Cenarchaeum symbiosum A",      "Desulfurococcus kamchatkensis 1221n", "Ferroglobus placidus DSM 10642",      "Halalkalicoccus jeotgali B3", "Haloarcula marismortui ATCC 43049",      "Halobacterium salinarum R1", "Halobacterium sp. NRC-1", "Haloferax volcanii DS2",      "Halomicrobium mukohataei DSM 12286", "Haloquadratum walsbyi DSM 16790",      "Halorhabdus utahensis DSM 12940", "Halorubrum lacusprofundi ATCC 49239",      "Haloterrigena turkmenica DSM 5511", "Hyperthermus butylicus DSM 5456",      "Ignicoccus hospitalis KIN4/I", "Ignisphaera aggregans DSM 17230",      "Metallosphaera sedula DSM 5348", "Methanobrevibacter ruminantium M1",      "Methanobrevibacter smithii ATCC 35061", "Methanocaldococcus fervens AG86",      "Methanocaldococcus infernus ME", "Methanocaldococcus jannaschii DSM 2661",      "Methanocaldococcus sp. FS406-22", "Methanocaldococcus vulcanius M7",      "Methanocella paludicola SANAE", "Methanococcoides burtonii DSM 6242",      "Methanococcus aeolicus Nankai-3", "Methanococcus maripaludis C5",      "Methanococcus maripaludis C6", "Methanococcus maripaludis C7",      "Methanococcus maripaludis S2", "Methanococcus vannielii SB",      "Methanococcus voltae A3", "Methanocorpusculum labreanum Z",      "Methanoculleus marisnigri JR1", "Methanohalobium evestigatum Z-7303",      "Methanohalophilus mahii DSM 5219", "Methanoplanus petrolearius DSM 11571",      "Methanopyrus kandleri AV19", "Methanosaeta thermophila PT",      "Methanosarcina acetivorans C2A", "Methanosarcina barkeri str. Fusaro",      "Methanosarcina mazei Go1", "Methanosphaera stadtmanae DSM 3091",      "Methanosphaerula palustris E1-9c", "Methanospirillum hungatei JF-1",      "Methanothermobacter marburgensis str. Marburg", "Methanothermobacter thermautotrophicus str. Delta H",      "Nanoarchaeum equitans Kin4-M", "Natrialba magadii ATCC 43099",      "Natronomonas pharaonis DSM 2160", "Nitrosopumilus maritimus SCM1",      "Picrophilus torridus DSM 9790", "Pyrobaculum aerophilum str. IM2",      "Pyrobaculum arsenaticum DSM 13514", "Pyrobaculum calidifontis JCM 11548",      "Pyrobaculum islandicum DSM 4184", "Pyrococcus abyssi GE5", "Pyrococcus furiosus DSM 3638",      "Pyrococcus horikoshii OT3", "Staphylothermus hellenicus DSM 12710",      "Staphylothermus marinus F1", "Sulfolobus acidocaldarius DSM 639",      "Sulfolobus islandicus L.D.8.5", "Sulfolobus islandicus L.S.2.15",      "Sulfolobus islandicus M.14.25", "Sulfolobus islandicus M.16.27",      "Sulfolobus islandicus M.16.4", "Sulfolobus islandicus Y.G.57.14",      "Sulfolobus islandicus Y.N.15.51", "Sulfolobus solfataricus P2",      "Sulfolobus tokodaii str. 7", "Thermococcus gammatolerans EJ3",      "Thermococcus kodakarensis KOD1", "Thermococcus onnurineus NA1",      "Thermococcus sibiricus MM 739", "Thermofilum pendens Hrk 5",      "Thermoplasma acidophilum DSM 1728", "Thermoplasma volcanium GSS1",      "Thermoproteus neutrophilus V24Sta", "Thermosphaera aggregans DSM 11486",      "Vulcanisaeta distributa DSM 14429", "uncultured methanogenic archaeon RC-I"     )      (Thanks)


title: How do I decide what span to use in LOESS regression in R?

text: I am running LOESS regression models in R, and I want to compare the outputs of 12 different models with varying sample sizes. I can describe the actual models in more details if it helps with answering the question. Here are the sample sizes:               Fastballs vs RHH 2008-09: 2002     Fastballs vs LHH 2008-09: 2209     Fastballs vs RHH 2010: 527      Fastballs vs LHH 2010: 449          Changeups vs RHH 2008-09: 365     Changeups vs LHH 2008-09: 824     Changeups vs RHH 2010: 201     Changeups vs LHH 2010: 330          Curveballs vs RHH 2008-09: 488     Curveballs vs LHH 2008-09: 483     Curveballs vs RHH 2010: 213     Curveballs vs LHH 2010: 162      The LOESS regression model is a surface fit, where the X location and the Y location of each baseball pitch is used to predict sw, swinging strike probability. However, I'd like to compare between all 12 of these models, but setting the same span (i.e. span = 0.5) will bear different results since there is such a wide range of sample sizes. My basic question is how do you determine the span of your model? A higher span smooths out the fit more, while a lower span captures more trends but introduces statistical noise if there is too little data. I use a higher span for smaller sample sizes and a lower span for larger sample sizes. What should I do? What's a good rule of thumb when setting span for LOESS regression models in R? Thanks in advance!


title: Using principal components analysis vs correspondence analysis

text: I am analyzing a data set concerning intertidal communities. The data are percent cover (of seaweed, barnacles, mussels, etc) in quadrats. I am used to thinking about correspondence analysis (CA) in terms of species _counts,_ and principle component analysis (PCA) as something more useful for linear environmental (not species) trends. I haven't really had any luck figuring out if PCA or CA would be a better fit for percent cover (can't find any papers), and I'm not even sure how something that is capped up to 100% would be distributed? I am familiar with the rough guideline that if the length of the first detrended correspondence analysis (DCA) axis is greater than 2, then you can safely assume that CA should be used. The length of DCA axis 1 was 2.17, which I don't find helpful.


title: How to assess change in obesity levels pre and post intervention where obesity is an ordinal variable?

text: I'm looking at changes in child obesity prevalence as well as other weight categories following an intervention. I'm trying to figure out which test to use for changes in a categorical dependent variable (I assigned values of 1-5 for children who were severely obese, obese, overweight, normal, underweight) pre and post intervention.


title: A survey of data-mining software tools

text: Although I was trained as an engineer, I find that I'm becoming more interested in data mining. Right now I'm trying to investigate the field further. In particular, I would like to understand the different categories of software tools that exist and which tools are notable in each category and why. (Note that I didn't say the "best" tools, just the notable ones lest we start a flame war.) Especially make note of the tools that are open-source and freely available - although don't take this to mean that I'm only interested in open-source and free.


title: How to test random effects in a multilevel model in R

text: I have been reading a good book called Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence by Judith Singer and John Willet. The book shows that by modeling in 2 levels, we can model the individual change in level 1 and in level 2 model for systematic interindividual differences in change. The R codes for the examples only show how to use `lme()` to estimate the fixed and random effects. However, the text suggested that we should test the variance components to determine whether the random effects are significant or not. For example, one of the codes does only the following:               library(nlme)          model.a <- lme(alcuse~ 1, alcohol1, random= ~1 |id)     summary(model.a)          Linear mixed-effects model fit by REML      Data: alcohol1             AIC      BIC    logLik       679.0049 689.5087 -336.5025          Random effects:      Formula: ~1 | id             (Intercept)  Residual     StdDev:   0.7570578 0.7494974          Fixed effects: alcuse ~ 1                      Value  Std.Error  DF  t-value p-value     (Intercept) 0.9219549 0.09629638 164 9.574139       0          Standardized Within-Group Residuals:            Min         Q1        Med         Q3        Max      -1.8892070 -0.3079143 -0.3029178  0.6110925  2.8562135           Number of Observations: 246     Number of Groups: 82      But the text lists the following:   * fixed effect: 0.922*** (s = 0.096) -> available in the output    * within person variance: 0.562*** (s = 0.062) -> can be obtained from the output (random effect residual std. dev squared)   * between person variance: 0.564*** (s = 0.119) My work involves a lot of analysis for longitudinal data so I really need to understand this idea. Your help is very much appreciated.


title: Rewriting AR model in State-Space form

text: How can I rewrite an AR(p) model in state-space form? Max(p)=5 and I want to use Kalman Predictor.


title: Measures of publication "importance" in statistics?

text: I am looking through articles citing an article I'm reading - and I wish to judge how much that citing article is "important" or "good" by itself. one way of knowing would be if I had known how "distinguished" that journal was. Which leads me to my question: **What measures are there (and where can I find them) for a journal "importance" or "impact"?** (I know of Impact Factor score. But I wonder if there are other such measures, and which are more applicable to statistical journals)


title: Logistic Regression - Multicollinearity Concerns/Pitfalls

text: In Logistic Regression, is there a need to be as concerned about multicollinearity as you would be in straight up OLS regression? For example, with a logistic regression, where multicollinearity exists, would you need to be cautious (as you would in OLS regression) with taking inference from the Beta coefficients? For OLS regression one "fix" to high multicollinearity is ridge regression, is there something like that for logistic regression? Also, dropping variables, or combining variables. What approaches are reasonable for reducing the effects of multicollinearity in a logistic regression? Are they essentially the same as OLS? (Note: this is not for the purpose of a designed experiment)


title: Efficient Empirical CDF Computation / Storage

text: I'm trying to precompute the distributions of several random variables. In particular, these random variables are the results of functions evaluated at locations in a genome, so there will be on the order of 10^8 or 10^9 values for each. The functions are pretty smooth, so I don't think I'll lose much accuracy by only evaluating at every 2nd/10th/100th? base or so, but regardless there will be a large number of samples. My plan is to precompute quantile tables (maybe percentiles) for each function and reference these in the execution of my main program to avoid having to compute these distribution statistics in every run. But I don't really see how I can easily do this: storing, sorting, and reducing an array of 10^9 floats isn't really feasible, but I can't think of another way that doesn't lose information about the distribution. Is there a way of measuring the quantiles of a sample distribution that doesn't require storing the whole thing in memory?


title: Computing inter-rater reliability in R with variable number of ratings?

text: Wikipedia suggests that one way to look at inter-rater reliability is to use a random effects model to compute intraclass correlation. The example of intraclass correlation talks about looking at $$\frac{\sigma_\alpha^2}{\sigma_\alpha^2+\sigma_\epsilon^2}$$ from a model $$Y_{ij} = \mu + \alpha_i + \epsilon_{ij}$$ "where Yij is the jth observation in the ith group, μ is an unobserved overall mean, αi is an unobserved random effect shared by all values in group i, and εij is an unobserved noise term." This is an attractive model especially because in my data no rater has rated all things (although most have rated 20+), and things are rated a variable number of times (usually 3-4). Question #0: Is "group i" in that example ("group i") a grouping of things being rated? Question #1: If I'm looking for inter-rater-reliability, don't I need a random effects model with two terms, one for the rater, and one for the thing rated? After all, both have possible variation. Question #2: How would I best express this model in R? It looks as if this question has a nice-looking proposal:               lmer(measurement ~ 1 + (1 | subject) + (1 | site), mydata)      I looked at a couple questions, and the syntax of the "random" parameter for lme is opaque to me. I read the help page for lme, but the description for "random" is incomprehensible to me without examples. This question is somewhat similar to a long list of questions, with this the closest. However, most don't address R in detail.


title: Data analysis for mixed within-subjects/between-subjects design?

text: I have run a survey with many respondents. There was a pool of questions, and each respondent was asked 3 of the questions (chosen uniformly at random and independently for each respondent). Now I'm interested in comparing the fraction of respondents who answered question $i$ correctly to the fraction who answered question $j$ correctly. What is the best methodology to do this? Keep in mind that a few of the respondents will have been asked both questions $i$ and $j$ (for them, this was a within-subjects design), while many will have been asked only one or the other (for this subclass of respondents, it is more akin to a between-subjects design). **Option #1.** I could count the total number of respondents who answered question $i$ and the number who answered it correctly to get the fraction who answered question $i$ correctly. Similarly, I could compute the fraction of respondents who answered question $j$ correctly. I could then compare these two numbers. Basically, I'd just be ignoring the fact that some users were asked both questions; don't worry, be happy. This is sort of a hybrid between a within-subjects/between-subjects design. **Option #2.** I could filter out all the respondents who were asked both questions $i$ and $j$. Of the remaining respondents who were asked either question $i$ or question $j$ (but not both), I could then compute the fraction who answered question $i$ correctly and the fraction who answered question $j$ correctly. This is essentially converting this to a between-subjects design. **Option #3.** I could keep only the answers from respondents who saw both questions, and then compare using a standard test for within-subjects designs. The downside is that I will have filtered out most of the responses; probably only a few respondents will have been randomly selected to answer both questions. **Option #4.** I could do something else (you tell me what). Do you have any advice? Or there are any pitfalls with any of these? Should any of these options be ruled out? Perhaps there is a way to combine the best of both worlds, and combine both Option #2 and Option #3 in some clever way?


Where can I find resources to learn about change-point analysis ? Hopefully, someone can advise me a **textbook** to read and it will cover both univariate change-point analysis and _multivariate_ change-point analysis. It will also be very helpful if someone can point me to some resources, which contains a worked example (using change-point analysis) using R. So, I understand the theoretical material and the applied material. Thank you very much.


title: Counterintuitive result for random sample

text: Say that you're making a political research, and you want to know how many people vote to Party A and how many vote to Party B. (There are only two choices.) You make a random survey of 10 people, and find that 1 out of your 10 votes A, while the other 9 vote B. What would you conclude is the average expected ratio of A voters in the general population? Most people's intuition would say `0.1`. But I calculated it and it came out as `0.0909090`, i.e. `1/11`. It's definitely possible I made a mistake in my calculation. Can you tell me whether my result makes sense?


title: Go from cluster to random forest

text: I'm working on improving a random forest model. I've done some clustering of the data using the `pvclust` package in R. I use ward and euclidean distances. My question is how do I go from the results of the clustering analysis to a better random forest. Do I need to run multiple `rf` objects with data split along the edges defined in the clustering results or is it something else. Any suggestions would be appreciated. Most of my work is in R.


title: A question on Poisson distribution

text: Perhaps my question is very simple, however I could not get a reliable solution for my problem. Here is my problem: Consider something as 'event' if this causes me a monetary loss equal to or above $x$ amount (say in USD). The number of such events (in a particular time interval, like say 10 years) follows a Poisson distribution with rate $\alpha$. Now my question is: Under this framework, what will be the distribution of the number of 'another event' which causes me $\$\ x_1$ loss, where $x_1 > x$. I believe that distribution will be Poisson with rate $\alpha\frac{x}{x_1}$? Is it correct? How can I prove this mathematically? Any suggestion will be highly appreciated.


title: Incidence rate ratios (Stata vs R)

text: I asked this question on Stackoverflow: http://stackoverflow.com/questions/8142118/incidence-rate-ratios-in-r-poisson- regression and was advised to post here instead. I have data that looks like this                  sex agecat  cases population          1 male    0-4  12     126526     2 male    5-9  12     128375     3 female  0-4  11     129280     4 male    10-14 4     127910     5 female  0-4  13     127158     6 male    0-4  8      125125      I want to duplicate the output I get in stata with this command               poisson cases i.agecat, exp(pop) irr      which gives output such as:                  cases |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]     ---------+----------------------------------------------------------------       agecat |           2  |   .5125755   .0530442    -6.18   0.000     .4578054    .6669639           3  |    .323456   .0381304    -9.60   0.000     .2665044    .4172274     population | (exposure)      in R with a command such as               glm(cases~agecat, family = poisson(link = "log")      I know I need to exponentiate the coefficients and confidence intervals, but I think I also need some kind of offset so so that the intercept is zero; and adjust for per unit population vs baseline. Can anyone help/advise ? Thanks EDIT: The question on SO has been answered, but I posted more detail here. In particular, I think the issue has to do with adjusting for population size in stata with exp(pop) - and how to replicate this in R.


title: How can I estimate the statistical significance of the optimal Gaussian Mixture model from Mclust?

text: I would like to estimate the statistical significance of a Gaussian Mixture model that is calculated by Mclust. The output of Mclust includes both the optimal Bayesian information criterion (BIC) used to select the optimal model as well as the Log Likelihood used in the calculation. However, as far as I understand these are only meaningful for comparing different models for the same data. Is there a way to estimate the statistical significant of the optimal model?


title: Comparing logistic regression models with same number of parameters

text: What would be an appropriate way to compare two logistic regression models with the same number of parameters (i.e., model 1 is not nested in model 2)? In my case, I cannot compare the models to a superset model with all predictors of both models, because including the differing parameters of both models results in nonidentifiability due to multicolliniearity.


title: Should availability (calculated using mean times to fail/repair) correspond to average proportion of units available over time?

text: If I have a mean time to repair (MTTR) of 3 days and a mean time between failures (MTBF) of 9 days, and so an availability of MTBF/(MTTR + MTBF) = 9/12= 75%, can I expect that, on average, 75% of my machines will be available per day, lets say I have a large number machines? Do any assumptions need to be upheld in order this to occur?


title: How to build event recommender based only on events' descriptions

text: I built a simple app that grabs events using eventbrite and meetup.com APIs and displays it based on your zip code. Now I'm trying to build a simple event recommender that will use events that you marked as interested/not-interested as a training set. Two main questions here:   * What would be the best way to extract features from events' descriptions?   * How big should be the training set? I used LDA to model topics based on descriptions so for every event there's a feature matrix which represents topics. Then I tried different classifiers linear SVC, KNearestNeighbors, RandomForest but nothing gives good results (not better than random classification). Is LDA a good way to extract features for this situation? Also the training set is just 356 events with 254 out of them marked negatively. Am I doing everything completely wrong? ;)


title: What is the expected value of modified Dirichlet distribution? (integration problem)

text: It is easy to produce a random variable with Dirichlet distribution using Gamma variables with the same scale parameter. If: $ X_i \sim \text{Gamma}(\alpha_i, \beta) $ Then: $ \left(\frac{X_1}{\sum_j X_j},\; \ldots\; , \frac{X_n}{\sum_j X_j}\right) \sim \text{Dirichlet}(\alpha_1,\;\ldots\;,\alpha_n) $ _Problem_ What happens if the scale parameters are not equal? $ X_i \sim \text{Gamma}(\alpha_i, \beta_i) $ Then what is the distribution this variable? $ \left(\frac{X_1}{\sum_j X_j},\; \ldots\; , \frac{X_n}{\sum_j X_j}\right) \sim \; ? $ For me it would be sufficient to know the expected value of this distribution.   I need a approximate closed algebraic formula that can be evaluated very very quickly by a computer.   Let's say approximation with accurancy of 0.01 is sufficient.   You can assume that: $ \alpha_i, \beta_i \in \mathbb{N} $ _Note_ In short, the task is to find an approximation of this integral: $ f(\vec{\alpha}, \vec{\beta}) = \int_{\mathbb{R}^n_+} \;\frac{x_1}{\sum_j x_j} \cdot \prod_j \frac{\beta_j^{\alpha_j}}{\Gamma(\alpha_j)} x_j^{\alpha_j - 1} e^{-\beta_j x_j} \;\; dx_1\ldots dx_n$


title: Testing for the existence of dependencies in time series

text: What methods exist to test for the existence of any sort of dependence in a time series? This is in contrast to something like auto-correlation, which tests for a particular type of dependency. Is there something that can say if dependencies exist, in general, without necessarily determining what kind? I have **heard of** and experimented with something called a "differential spectrum test", which states that the histogrammed _changes_ in the time series should be symmetric about zero if the series of changes is independent. (i.e. histogram the set of $x_{t+1} - x_{t}$ and it should be symmetric about zero). If anyone has a reference for this test, I would appreciate it. Update: The comments suggest that failing a white noise test would indicate that some sort of dependencies exist. True? I guess it must depend on the test. Here are some I have found.. any comments on their applicability here? Bartlett and Q Tests Welch's Method Method of Perez et. al. Method of Lobato and Velasco


title: Determine optimal (recurring) frequency of discrete values

text: So I'm trying to determine a semi-efficient way to calculate the optimal (recurring) frequency of a set of data. The data only exists at random time periods, but it is assumed that the data set is complete. **For example:** Say that we had data recorded every time a movie was watched, for which the viewer recorded their liking of the movie. The likability obviously depends on what movie was watched. However, what I am interested in is the frequency at which movies should be watched for optimal likability. Furthermore, even if two movies were watched one after the other, and the viewer loved them, this should not be the optimal frequency unless:   1. it is proven that the viewer could continuously watch movies at that rate.    2. the data is not an outlier, and the results have been reproduced elsewhere in the data set.  How would one go about doing this? Note that in this scenario, the data could be processed as it was forming (we could tell the viewer to try watching movies at x frequency, for optimal results, but there is no guarantee that they will do it). Thanks!


title: Finding local extrema of a density function using splines

text: I am trying to find the local maxima for a probability density function (found using R's `density` method). I cannot do a simple "look around neighbors" method (where one looks around a point to see if it's a local maximum with respect to its neighbors) as there is a large volume of data. Furthermore, it seems more efficient and generic to use something like Spline interpolation and then find the roots of the 1st derivative, as opposed to building a "look around neighbors" with fault tolerance and other parameters. So, my questions:   1. Given a function from `splinefun`, what methods will find the local maxima?    2. Is there an easy / standard way to find derivatives of a function returned using `splinefun` ?   3. Is there a better/standard way to find the local maxima of a probability density function ? For reference, below is a plot of my density function. Other density functions I'm working with are similar in form. I should say that I am new to R, but not new to programming, so there may be a standard library or package for achieving what I need. ![density function](http://i.stack.imgur.com/eF6ot.png) Thanks for your help!!


title: Compare one distribution with multiple other distribution

text: I've a distribution of numeric scores computed using real data ( my reference distribution ). After that I generated 100000 random distributions of numeric scores computed using points picked at random in a space of possible points. So now I want to compare my reference distribtuion with these random distributions and to test whether the reference distribution is significantly different from random distributions. Thanks


title: Variable A predicts B, variable B predicts C

text: I have pre - post data on quality of life, depression, and self-esteem. I hypothesize that changes in self-esteem from pre to post cause changes in depression, which in turn cause changes in quality of life, from pre to post. This looks to me like two separate multiple regressions, but I am unsure whether this could be combined into a single analysis or what it's the same. I don't know what kind of analysis would be best to answer my hypothesis.


title: F-test: Testing unrestricted model vs. restricted model in Gretl and SPSS

text: I have a proposed a model for Corporate Social Responsibility determinants (not important). Among my independent variables I obtain two that are statistically not significant. Following Wooldridge approach, in order to know if I should take out these two variables from the model, I should run a F-test, with the restricted and unrestricted model, right? My question is: can I do this straighaway with a software like Gretl or SPSS? Or should I construct the statistic manually? Thanks for the reply


title: Perform simple regression without raw data

text: I have a dataset that I can collect some quantities from, eg. sum,mean,variance... I want to perform a simple regression on column(x,y). According to Wikipedia, the closed form for $\alpha,\beta$ is \begin{equation} \alpha=\frac{Cov(x,y)}{Var(x)} \quad \mbox{and} \quad \beta=\bar{y}-\beta\bar{x}. \end{equation} However I would also like some more stats from the regression, for example **r, r-square, significance...** Basically I want to get the same result that this apache SimpleRegression class produces by only feeding collected quantities instead of raw data points. Any advice on how to do that?


title: Time series: Stationary R-squared vs R-squared (Forecasting in SPSS)

text: I would like to check how two time series are correlated - how one of them predicts the other one. I have data for 27 years. I conduct my analysis in SPSS (Forecasting command, ARIMA) and I get: ![enter image description here](http://i.stack.imgur.com/lMJpN.png) Could somebody please tell me what the difference between the Stationary R-squared and R-squared is and which of them is better in this case? Judging by the R-squared the model is good I guess, whereas based on Stationary R-squared it's so-so. I've read that stationary R squared is better if there's a trend (do they mean one simple trend?) or seasonal changes. In my case both values have both increased and decreased at different points in time: there were periods in which the values were persistently high/ increased and periods in which they were persistently low/ decreased. So there're complex dynamics, hard to tell whether there's one trend. Does it mean R squared is enought? On the other hand, I read that R-squared is mainly appropriate to analyse stationary data. My data in series are definitety autocorrelated, so if I understand it correctly, it's non-stationary. So do I have to use the stationary R-squared? Thanks in advance! (Sorry, I notice now I didn't change all commas into points in the table. I'm using a non-English version of SPSS. But commas equal points)


title: Variable bounds from infinity to negative infinity - Change Limits

text: I have a variable whose values vary between $-\infty$ and $\infty$.   I'd like to change its limits from $0$ to $1$. I want to give extremely low values when the values are less than $-1$ or greater than $1$ and want to give the highest value to the value $0$. Which distribution would be the best here, and what's the formula to do so?   And if possible, can you guide me to some good links where I can learn more about such stuff?


title: How to compute the marginal probability form conditional probabilities in logspace?

text: Normally the marginal probablity is computed as $p(x) = \sum_y p(x | y) \cdot p(y) $ Now, suppose I have all these probabilities at the right-hand side in logspace (so as logprobabilities). How do I do this?!?


title: Faith in an extrapolated result

text: I would like to be able to predict when I will exhaust a particular resource. My situation is analogous* to a water tank. Each day zero or more rain will fall, filling up the tank. I can not tell in advance how much rain will fall each day, if any at all. There will never be any outflow from the tank. Each day I can measure how much is in the tank and I know that tank's capacity. Since there is a cost to having spare tanks on hand, just in case, and also to filling a tank when no spare is available, I'd like to use my daily measurements to predict when the tank will be full. Experience so far suggests the tank has about four year's capacity. It will take us something like three to six months to get a new tanks set up and working depending on what else is going on. It is something we want to plan for in advance. Of course emergency tanks can be rushed through if a deluge is imminent, but the cost is huge and this should be avoided. I've calculated a linear regression of the daily water levels. This gives me a slope and intercept from which I can calculate the date when the trend line will cross the tank's capacity. This gives me no indication of the faith I can place in this number, however. Is it legitimate to simply use the Confidence Interval formula to similarly calculate when the CI curve crosses the capacity? Is there a better way of achieving this? The data will look like this (invented values): ![Example with trend.](http://i.stack.imgur.com/bW8Dv.png) * * * *The actual case is a data warehouse computer system. The rows have an auto-incrementing column. As a new row is written the value in this column goes up. It never goes down. There is no equivalent of "draining" the water-tank analogy. If I measure the largest value in this column today, and no rows are added overnight, I will get the same vlaue when I ask tomorrow. If rows were written overnight, tomorrow's measurement will be larger than today's. Even if all rows are removed the auto-increment will not reset to zero. The column is defined as an integer which has a maximum value of just over two billion. If we reach this limit we cannot add rows to the table, we have system down time, my boss's head explodes, the oceans turn to blood etc. We don't want that to happen. The disks storing this data cost tens of thousands of dollars, however, so we don't want to have the capital outlay any earlier than we have to, obviously. The definition of "have to" is the thrust of my question.


title: calculate the bias and variance of the MLE in an exponential family

text: I'm trying to calculate the bias and variance of $\hat{\lambda}_{MLE}$ where $X \sim exponential(\lambda)$. Secondly, does it attain the CRLB? Attempt: $$f_x(X) = \lambda e^{-\lambda x}$$ $$L(\lambda; x_1,...,x_n) = \lambda^n e^{-\lambda \sum^n X_i}$$ $$\ell(\lambda; x_1,...x_n) = n ln(x) - \lambda \sum^n X_i$$ $$\hat{\lambda}_{MLE} = \frac{\sum^n X_i}{n} = \lambda^{-1}$$ As $\frac{d}{d\lambda} \ell(\lambda; x_1,...,x_n) = 0$ and $\frac{d^2}{d\lambda^2} \ell(\lambda; x_1,...,x_n) < 0$ when $\hat{\lambda}_{MLE} = \frac{\sum^n X_i}{n}$, I believe that the $\hat{\lambda}_{MLE}$ is unbiased. I am stuck on finding the variance of $\hat{\lambda}_{MLE}$. Since the MLE is unbiased I believe it attains CRLB.


title: Predictive Models for the Soccer World Cup 2014

text: The Goldman Sachs model was already discussed here. I wonder whether there are other statistical forecasts publicly available, or even better, models including raw data. Specifically, I thought that the Goldman Sachs model didn't work so well for the 2010 cup (see the original report on page 7). I also heard that actual or estimated transfer values of players have high predictive power, but they are not included in the Goldman Sachs model.


title: How can I analyze "win rate" over time (i.e. on/off data over time)?

text: I have a series of data points that consist of 1) a time, and 2) a win or a loss. I would like to be able to determine the aggregate win rate for particular time periods, and graph it. For instance, let's say I am dealing with boxers. I would like to be able to say the boxer Bill tends to peak in boxing matches around the 30 mark, with a win rate of 70%, and tends to fall off around the 40 minute mark, with a win rate of 40%. What statistical methods or tools should I use to get insight into this?


title: Reduce Random Forest Model Size

text: I've created a regression model on my data using random forests in R. The output is quite large, I'm wondering if there's any way to reduce this to only the necessary pieces to make a prediction? The training data set contains 20 variables and ~45,000 rows, which is also large. My code is listed below.               data <- readRDS("data.Rds")          require("data.table")     require("doParallel")     require("randomForest")          train <- data[ which(set == "train")]     test <- data[ which(set == "test")]     rm(data)          x <- data.table(train[, 2:21, with=FALSE])     y <- as.vector(as.matrix(train[, 23, with=FALSE]))          cl <- makeCluster(detectCores())     registerDoParallel(cl, cores=4)     time <- system.time({rf.fit <- foreach(ntree=rep(500, 6),                                    .combine=combine,                                    .multicombine=TRUE,                                    .packages="randomForest") %dopar%                         {randomForest(x, y, ntree=ntree)}})     stopCluster(cl)          saveRDS(rf.fit, "rf.fit.Rds")      The output of this is ~230 MB. Once I have the model, is it possible to reduce the size to make it easier to work with? My goals with this are to identify the important variables, and make a prediction on new data.


title: Longitudinal Monitoring: How to test if a parameter changes before the gold standard

text: I have some data that I am really uncertain how to analyse. Each subject who all begin as healthy/normal will undergo a series of imaging examinations during therapy: scan 1 (baseline), scan 2, scan 3, ... scan N. These scans will then be processed and the gold standard parameter, say G, is measured. Now if G changes by 10% this patients therapy will be suspended. What I would like to do is aswell as measuring G for each scan for each subject I would like to measure several other novel parameters (say A,B,C and D) to see if they are predictors that G will drop by 10% in a subsequent scan. Something like A would fall by a certain amount in scan 3 which would predict a significant (10% or greater) fall in G. So far I have split my data into two sets: Those subjects who experience a clinical significant drop in G (abnormal) and those who do not (normal). The change in G can occur at any scan in the series. Any ideas on how best to approach this analysis?


title: Moment Generating Functions and Fourier Transforms?

text: Is a moment-generating function a Fourier transform of a probability density function? In other words, is a moment generating function just the spectral resolution of a probability density distribution of a random variable, i.e. an equivalent way to characterize a function in terms of it's amplitude, phase and frequency instead of in terms of a parameter? If so, can we give a physical interpretation to this beast? I ask because in statistical physics a cumulant generating function, the logarithm of a moment generating function, is an additive quantity that characterizes a physical system. If you think of energy as a random variable, then it's cumulant generating function has a very intuitive interpretation as the spread of energy throughout a system. Is there a similar intuitive interpretation for the moment generating function? I understand the mathematical utility of it, but it's not just a trick concept, surely there's meaning behind it conceptually?


title: Time series for network matrices?

text: Let's say I have two people. Each of whom I create a network matrix for daily for some of their daily habits $(i.e. A, B, C, D)$. I do this for, let's say $60$ days. ![Two People and their daily matrices](http://www.chrisjvargo.com/Matrix%20Time%20Series.png) I didn't bother to change the values in the matrices, but let's assume they vary from day to day. At the end of the $60$ days, I want to see if `Person 1's` matrix of daily habits influenced `Person 2's` matrix. If this wasn't a matrix, I'd do something like an ARIMA time series. Is there a way to process matrices and see if `Person 1` predicts Person 2 or vice versa? Open to any statistic/package that might get the job done.


title: How was this p-value calculated?

text: Sorry to ask such a blunt question, but I'm running out of time and I just can't figure this seemingly simple thing out. It's the p-value for $B_0$, a single variable regression. $H_0: B_0=0$ and $H_1: B0 \neq 0$. $T_{calc} = (0.5681 - 0)/0.9284 = 0.6119$. Then it says $p\text{-value}=2 \cdot (t\ge 0,6119) =2 \cdot 0.281525$. Where is this $0.281525$ coming from? I know it has something to do with the table. $n=8$. Sorry again, I have a test coming up, for some reason I can't figure this out. I spent way too much time on this now, watching youtube videos for the past 1.5 hours trying to figure out where this .2815... came from.


title: Forecasting the business cycle?

text: I am wondering what is the best way to forecast the business cycle based on the past. Currently I feed the seasonally-adjusted GDP index data to a Hodrick–Prescott filter, extract the cyclical component, take the last data point, find the similar points in the past and see what happened afterwards. (See attached picture.) Perhaps there is a better way to do the extrapolation? ![enter image description here](http://i.stack.imgur.com/EpnIZ.jpg)


title: Get the distribution of a dataset

text: this is a conceptual question. In data mining, the problem often arises that scientists/data engineers are using an expert guess for the underlying distribution of their data. Often their assumption is based on previous results of academic research and does not really reflect the real distribution of the dataset. However, I would like to get the real distribution of my dataset? Any recommendations where to start? Does this even make sense to get the real distribution? Can thisv task be automated? I appreciate your answer! btw I am mostly using R as my data analysis tool


title: How to perform pooled cross-sectional time series analysis?

text: For 86 companies and for 103 days, I have collected (i) tweets (variable `hbVol`) about each company and (ii) pageviews for the corporate wikipedia page (`wikiVol`). The dependent variable is each company's stock trading volume (`stockVol0`). My data is structured as follows:               company  date  hbVol    wikiVol   stockVol0  comp1  comp2 ... comp89  marketRet     -------------------------------------------------------------------------------     1        1     200        150     2423325      1      0   ...   0     -2.50     1        2     194        152     2455343      1      0   ...   0     -1.45     .        .      .          .         .         .      .   ...   .     1       103    205        103     2563463      1      0   ...   0      1.90     2        1     752        932     7434124      0      1   ...   0     -2.50     2        2     932        823     7464354      0      1   ...   0     -1.45     .        .      .          .         .         .      .   ...   .     .        .      .          .         .         .      .   ...   .     86      103     3          55      32324       0      0   ...   1      1.90      As I understood, this is called pooled cross-sectional time series data. I have taken the Log-value of all variables to smoothen the big differences between companies. A regression model with both independent variables on the dependent `stockVolo` returns: ![enter image description here](http://i.stack.imgur.com/Ud6hV.png) A Durbin-Watson of 0,276 suggest significant autocorrelation of the residuals. The residuals are, however, bellshaped, as can be seen from the P-P plot below. The partial autocorrelation function shows a significant spike at a lag of 1 to 5 (above upper limit), confirming the conclusions drawn from the Durbin-Watson statistic: ![enter image description here](http://i.stack.imgur.com/kvMce.png) The presence of first-order autocorrelated residuals violates the assumption of uncorrelated residuals that underlies the OLS regression method. Different methods have been developed, however, to handle such series. One method I read about is to include a lagged dependent variable as an independent variable. So I created a lagged `stockVol1` and added it to the model: ![enter image description here](http://i.stack.imgur.com/KQQz7.png) Now, Durbin-Watson is at an accceptable 2,408. But obviously, R-squared is extremely high because of the lagged variable, see also the coefficients below: ![enter image description here](http://i.stack.imgur.com/2vHzf.png) Another method I read about when being confronted with autocorrelation, is autoregression with Prais-Winsten (or Cochrane-Orcutt) method. Once performing this the model reads: ![enter image description here](http://i.stack.imgur.com/zlIUP.png) This is what I don't understand. Two different methods, and I get very different results. Other suggestions for analyzing this data include (i) not including a lagged variable but reformat the dependent variable by differencing (ii) perform AR(1) or ARIMA(1,0,0) models. I haven't calculated those because I am now lost on how to proceed because of the different results of the two tests I did perform. What model should I use to perform a proper regression on my data? I'm very keen on understanding this, but have never had to analyze a timeseries dataset like this before.


title: Analysis to use with categorical predictors and continuous response variables

text: I am trying to figure out what the right statistical test to use is for data where the DV is continuous, and the IV's are categorical. There are a lot of different DV's that I want to run through a model one at a time, so I am having a hard time figuring out how to do exploratory statistics to really figure out the right next steps (i.e. methodology). I had assumed regression would work, but I am worried about the requirement for linearity. Does this requirement still hold when the IV's are categorical vs continuous. Can someone please help me figure out the right test to use.


title: Correlation between predictor variables in an AIC model

text: I'm using multinomial logistic regression analysis to analyse deer behavioural responses to camera traps based on 7 predictor variables. I have 2 models which are very close together in AIC value (season AIC 1005.103; season and percentage cover 1005.023). Season and percentage cover is the best model (but only nominally). Percentage cover is defined as the percentage of the deer covered by topographical features and vegetation relative to the camera. It has an AIC value of 1008.289 when considered alone. However, season and cover appear to be closely correlated (which makes sense). I am wondering if it is redundant to use both variables in the same model or is a certain level of correlation okay? Thanks a mil


title: Get positions/values from heat costs bill

text: for a project I need to extract values from customers yearly heat costs bill. The customer takes a photo of the bill and the program should extract the values heating period of the billing, type of energy source (gas, fuel etc.) and the amount of energy used. I have little knowledge of machine learning (one lecture at university). It's not an easy task, is it? Can it be done with open source software (mostly Python I guess)?


title: Modelling a binary response variable with a binary explanatory variable

text: I'm modeling (using GLM) one binary explanatory variable against a binary response variable. The output looks ok to me, but I'm a biologist with a basic knowledge of stats, I just wanted to ask if there's any rule/caveat/problem for modeling a binary versus a binary?


title: How to define and model consumption bundles?

text: Imagine an a la carte buffet with n different rooms. On entering the buffet you pick a room (let's say American, Mexican or Italian food) where you stay for the duration of your visit. Once in a room, you have a full selection of small dishes to pick from and most people choose multiples of several different dishes, e.g. 1 burrito, 2 salad, 4 tacos, 2 churros. Presumably these are neither independent nor random choices, but you can't get into the heads of your customers. Also, although the rooms are themed (and people are more likely to order according to the theme they've chosen) every dish is available in every room. How would you model the distribution of dishes in a visit to a room, given complete data? How could you quantify the deviation from expectation of a given visitor's selection? Because the dish selection isn't necessarily independent I think I would do the first according to a Dirichlet process, but I'm not sure that's right.


title: Does Sampling Frequency matter for Time Series analysis?

text: I am given two time series of prices between 2009 and 2013. Price series A is weekly data, series B is monthly data. I would like to compare some basic descriptive statistics of these two time series data (such as mean, standard deviation etc). Intuitively, I just wonder if the period sampling difference should be of any concern.... What could be some potential problems, and if there are ways to reduce the impact.... I guess the other way to put it is that, should I just pick the month-end numbers of series A to match the frequency of series B. I am not sure about this approach because the quality of statistics might be compromised for series A. Thank you so much!


title: ANOVA to test overall significance of model in linear regression with only continuous predictors

text: A normal ANOVA model always needs a categorical variable as a predictor(the factor/group) right? Please let me know the reason if you think otherwise. So in case of a linear regression with only continuous predictors, how can the ANOVA table be even used since ANOVA requires at least 1 categorical variable as predictor? Thanks


title: How to analyze this messy design?

text: I need to analyze a data-set, with a very messy design, I am not sure how. I will try to make it simple. A new kind of stitches was invented, and is tested vs. 2 old kind of stitches. I will call this: Treatment, control 1 and control 2. A few patients were selected, needing stitches. Each patient was assigned (by the need, not by random) to one of two stitching techniques. There are others, but only these two were tested. I call this variable "procedure type". Each patient, needed one stitch or more, up to 5 stitches. The problem is, that those patients in need for more than one stitch, received (sometimes) a mix of different stitches. Meaning, that one patient could have treatment and control 1, and other could have 2 treatment and 2 control 2, etc... I hope I am not mixing terms here, but my experimental unit is the patient and the observations unit is the stitch within the patient. I have 40 treatment stitches, and 20 of each control, all coming from 25 patients. The measure being tested, is how many minutes did it take for each stitch until the bleeding stopped. The variable can take values: 0,0.5,1,1.5,2,2.5,... Clinicians claim that the procedure type, i.e. the technique being used is of no clinical importance and should not have any affect. Thus they think the analysis could ignore this factor. I wanted to ask you, how would you analyze this data, with and without taking into account the technique. I think it should be some sort of a mixed model or a generalized estimating equations model, but I am not sure to define this design exactly in terms of "nesting", "blocking" and so on, and I would like to also fit a SAS code to model the differences between the means of the 3 different treatments, and I wouldn't mind verifying it with R. For your convenience, I attach a diagram I made of the design: the blue circles are treatment stitches, the red ones are control 1 stitches, and the green are control 2 stitches. Each yellow rectangle is a patient, and the big brown rectangles are the techniques, which again, I would like to try to fit and to try to ignore the see the outcome of that. ![enter image description here](http://i.stack.imgur.com/l7O51.png)


title: Trying to run statistical tests in R but struggling as I am new to the language

text: Good Day, I believe this issue is more of a lack of understand of R (as I have never used it till recently) than anything else. What I am looking for is references, or documents to help me solve my problem. I have a set of data that I created a pair of decision tress on (pruned and unpruned). I did this using Weka and 10-fold cross-validation. There was no problem here. I got a higher accuracy for the pruned tree and I have individual accuracies for each of the ten folds. I am now working on doing statistical testing to these two trees, I want to do a T-Test and the Wilcoxon test. It was suggested that I use the               experimentalComparison()      function for the cross-validation, but I can already do this in weka. For the actual stat tests I can do the t-test similar to the example here by manually feeding the values from the individual fold-tests similar to what is in the code in that example. The Wilcoxon test is what is throwing me off. I found this example but I was told to use the function               compAnalysis()       in the DMwR package. I looked at the documentation for it and I do not quite understand it. It says to feed in the results from `experiemntalComparison()` for the first value, which means I should be doing my cross-validation with that function (vs. weka). I am not understanding how to get the "trees" I made in weka this function. The documentation refers to some usermade class `cv.rpartXse`. So I am stumped, being new to this language, on how to accomplish these statistical tests via R.


title: If $Y \sim geometric(P)$ and $P \sim \mathcal B(2, 1)$ how to compute $E(Y)$ and marginal pmf of $Y$?

text: $$Y \sim Geometric(P)\\ P \sim \mathcal B(2, 1)$$ I'm trying to compute $E[Y]$ without finding marginal distribution of $Y$. I need some hints here. I also need to find the pmf of $Y$. My approach is as follows: Integrate from 0, 1 over $p: f(p)(1-f(p))^{(y-1)}$ where $f(p)$ is the pdf of $\mathcal B(a, b)$. This leads to some messy algebra. Am I on the right track?


title: Good resources for learning Probabilistic Graphical Models

text: I recently started taking Probabilistic Graphical Models on coursera, and 2 weeks after starting I am starting to believe I am not that great in Probability and as a result of that I am not even able to follow the first topic (Bayesian Network). That being said I want to make an effort to learn this course, so can you suggest me some other resources for PGM or for Probability which can be helpful in understanding this course.


title: Simulation setting for the average treatment effect on the treated (ATT)

text: I have a question about a simulation set up. Assume there are two groups (Z = 0 and Z = 1). The outcome for Z = 0 and Z = 1 are generated by the following equations: $Y_0 = \alpha_0 + \beta_{01}*{x_1} + \beta_{02}*{x_2} + \epsilon_{0}$ $Y_1 = \alpha_1 + \beta_{11}*{x_1} + \beta_{12}*{x_2} + \epsilon_{1}$, where $x_1$ and $\epsilon$ are standard normal $N(0,1)$ and $x_2$ is bernoulli distributed with $\pi = 0.5$. The true values of ($\alpha_0$, $\beta_{01}$, $\beta_{02}$) = (0, 1, 2) and ($\alpha_0$, $\beta_{01}$, $\beta_{02}$) = (1, 2, 3). How can I know the true value of the average treatment effect on the treated (ATT $ = E(Y_1 - Y_0 \mid Z =1)$) based on the above setting? many thanks in advance.


I have been calculating p-values for correlations using bootstrapping. I've just been using the same p-value cutoff regardless of the number of pair-wise calculations. Intuitively, it seems to me that the more correlations that I perform, the more likely I am to obtain a spurious correlation due to random chance. **Does using bootstrapping change how you deal with problems of Type I errors when testing multiple correlations?**


title: k-fold cross-validation strategy for large data set in statistical learning

text: I'm trying to learn the Bayesian network structure from a very large data set, and the R package I used for learning can only handle a very small portion of the data set (~10%) at one time due to the computational limitations. The conventional k-fold cross-validation strategy uses k-1 subsets for training and 1 subset for testing. I want to know if I can use only one random subset for training and another random subset for testing? Is there any better solution?


title: Proper use of weights argument in linear mixed model [lme()]?

text: After attempting to produce a linear mixed model I was left with a great deal of heterogeneity.               lme1 <- lme(Average.payoff ~ Game + Type + Others.Type + Game:Type +        Game:Others.Type + Type:Others.Type, random=~1|Subjects, method="REML",        data=Subjectsm1)      The response term `Average.payoff` is continuous, whilst all explanatory variables are binary. When I look to validate, I can clearly see that the spread of the residuals decreases with the larger fitted values. So, I decided to look at using the weights argument. I came up with the following model which seems to be the best fit (`AIC` and `loglik`).               vf1 <- varComb(varIdent(form=~1|Others.Type), varIdent(form=~1|Game),            varIdent(form=~1|Type), varIdent(form=~1|Male))          lme1 <- lme(Average.payoff ~ Game + Type + Others.Type + Game:Type +        Game:Others.Type + Type:Others.Type, random=~1|Subjects, method="ML",        weights=vf1, data=Subjectsm1)        1. Can I justify using the varIdent function alongside random effects?   2. Can I justify combining varIdent functions so the variance of all explanatory variables are accounted for?   3. If I remove a fixed effect explanatory variable for being insignificant should the same variable be removed from the variance structure?   4. If I get a better fit (using `AIC` and `loglik`) without including the random effects by re-modelling with `gls()`, can I justify using the model without the random effects? The experiment involved $20$ subjects, for each subject I have two samples, one where `Others.Type=0` and one where `Others.Type=1`.


title: Advice on regression modelling

text: I have a large set of data for 37 different clinical units (all oncology) in their respective 37 hospitals. There are two specific outcome variables that I need to analyse: First, drug usage for specific drugs types and classes (aggregated drugs) that are expressed as a rate – DDD (Defined Daily Doses) per 100 patient days. There are individual patient drug use figures for this set. Question1: Which regression approach should I take? From what I can gather I can use a Poisson regression model. IF there is overdispersion in the outcome I could resort to a negative binomial model. Second: I have antibiotic resistance data that is expressed as proportion in the range 0 – 1.These are not available as individual patient data points but aggregated to each of the 37 hospitals. Question 2: Again, which approach? From what I have read I can use a logistic regression model. I have been advised by another statistician to initiall use a logit model and then use a probit model and compare goodness of fit for each model. Does this sound like a reasonable approach? Is there a specific text that you could direct me to in order to upgrade my basic regression modelling skills. I will be using R for the analysis. Thanks in advance.


title: Calculation of an overall top X over several groups

text: I'm not a statistics expert, so I'm hoping someone here can lend me a hand. I've got a bunch of key-value pairs associated with a specific time range. Something like this:               <time>|<key>|<value>     0     |1    |34534     0     |2    |23434     0     |3    |4606     1     |1    |945954     1     |6    |459459     1     |8    |34      There will be 10's of millions of key-value pairs over a variety of time values (24 unique values). I need to be able to efficient calculate the "Top 10" sum of values across all time values grouped by the key value. I don't think that I'll be able to efficiently do this unless I split the problem. So, my theory is that I can split the calculation across each individual value of time. In the example above, that would give me two jobs:               Job1     <time>|<key>|<value>     0     |1    |34534     0     |2    |23434     0     |3    |4606          Job2     <time>|<key>|<value>     1     |1    |945954     1     |6    |459459     1     |8    |34      Now imagine that I have millions of key-value pairs for each time value. I calculate say the Top 1000 sums ordered by the value in a descending fashion. I then aggregate all of the Top 1000 job sums together to calculate the "Top 10" of the full time range. I'm using 1000 as an example here, but what is a more precise value that will guarantee my "Top 10" will be correct? I'm not even quite sure where to start with this. Vinbot